diff --git a/applications/utilities/thermophysical/chemkinToFoam/chemkinToFoam.C b/applications/utilities/thermophysical/chemkinToFoam/chemkinToFoam.C
index cb219fb83..e094ee28a 100644
--- a/applications/utilities/thermophysical/chemkinToFoam/chemkinToFoam.C
+++ b/applications/utilities/thermophysical/chemkinToFoam/chemkinToFoam.C
@@ -67,15 +67,13 @@ int main(int argc, char *argv[])
 
     chemkinReader cr(species, args[1], args[3], args[2], newFormat);
 
-
     OFstream reactionsFile(args[4]);
-    reactionsFile
-        << "elements" << cr.elementNames() << token::END_STATEMENT << nl << nl;
-    reactionsFile
-        << "species" << cr.species() << token::END_STATEMENT << nl << nl;
+    reactionsFile.writeKeyword("elements")
+        << cr.elementNames() << token::END_STATEMENT << nl << nl;
+    reactionsFile.writeKeyword("species")
+        << cr.species() << token::END_STATEMENT << nl << nl;
     cr.reactions().write(reactionsFile);
 
-
     // Temporary hack to splice the specie composition data into the thermo file
     // pending complete integration into the thermodynamics structure
 
@@ -103,6 +101,14 @@ int main(int argc, char *argv[])
 
     thermoDict.write(OFstream(args[5])(), false);
 
+    reactionsFile << nl;
+
+    reactionsFile.writeKeyword("Tlow")
+        << Reaction<gasHThermoPhysics>::TlowDefault
+        << token::END_STATEMENT << nl;
+    reactionsFile.writeKeyword("Thigh")
+        << Reaction<gasHThermoPhysics>::ThighDefault
+        << token::END_STATEMENT << nl << nl;
 
     Info<< "End\n" << endl;
 
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/StandardChemistryModel/StandardChemistryModel.C b/src/thermophysicalModels/chemistryModel/chemistryModel/StandardChemistryModel/StandardChemistryModel.C
index ad2c75577..adc979f62 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/StandardChemistryModel/StandardChemistryModel.C
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/StandardChemistryModel/StandardChemistryModel.C
@@ -56,7 +56,7 @@ Foam::StandardChemistryModel<ReactionThermo, ThermoType>::StandardChemistryModel
         BasicChemistryModel<ReactionThermo>::template lookupOrDefault<scalar>
         (
             "Treact",
-            0.0
+            0
         )
     ),
     RR_(nSpecie_),
@@ -80,7 +80,7 @@ Foam::StandardChemistryModel<ReactionThermo, ThermoType>::StandardChemistryModel
                     IOobject::NO_WRITE
                 ),
                 thermo.p().mesh(),
-                dimensionedScalar("zero", dimMass/dimVolume/dimTime, 0.0)
+                dimensionedScalar("zero", dimMass/dimVolume/dimTime, 0)
             )
         );
     }
@@ -109,8 +109,6 @@ void Foam::StandardChemistryModel<ReactionThermo, ThermoType>::omega
     scalarField& dcdt
 ) const
 {
-    scalar pf, cf, pr, cr;
-    label lRef, rRef;
 
     dcdt = Zero;
 
@@ -118,24 +116,7 @@ void Foam::StandardChemistryModel<ReactionThermo, ThermoType>::omega
     {
         const Reaction<ThermoType>& R = reactions_[i];
 
-        const scalar omegai = omega
-        (
-            R, c, T, p, pf, cf, lRef, pr, cr, rRef
-        );
-
-        forAll(R.lhs(), s)
-        {
-            const label si = R.lhs()[s].index;
-            const scalar sl = R.lhs()[s].stoichCoeff;
-            dcdt[si] -= sl*omegai;
-        }
-
-        forAll(R.rhs(), s)
-        {
-            const label si = R.rhs()[s].index;
-            const scalar sr = R.rhs()[s].stoichCoeff;
-            dcdt[si] += sr*omegai;
-        }
+        R.omega(p, T, c, dcdt);
     }
 }
 
@@ -156,123 +137,11 @@ Foam::scalar Foam::StandardChemistryModel<ReactionThermo, ThermoType>::omegaI
 ) const
 {
     const Reaction<ThermoType>& R = reactions_[index];
-    scalar w = omega(R, c, T, p, pf, cf, lRef, pr, cr, rRef);
+    scalar w = R.omega(p, T, c, pf, cf, lRef, pr, cr, rRef);
     return(w);
 }
 
 
-template<class ReactionThermo, class ThermoType>
-Foam::scalar Foam::StandardChemistryModel<ReactionThermo, ThermoType>::omega
-(
-    const Reaction<ThermoType>& R,
-    const scalarField& c,
-    const scalar T,
-    const scalar p,
-    scalar& pf,
-    scalar& cf,
-    label& lRef,
-    scalar& pr,
-    scalar& cr,
-    label& rRef
-) const
-{
-    const scalar kf = R.kf(p, T, c);
-    const scalar kr = R.kr(kf, p, T, c);
-
-    pf = 1.0;
-    pr = 1.0;
-
-    const label Nl = R.lhs().size();
-    const label Nr = R.rhs().size();
-
-    label slRef = 0;
-    lRef = R.lhs()[slRef].index;
-
-    pf = kf;
-    for (label s = 1; s < Nl; s++)
-    {
-        const label si = R.lhs()[s].index;
-
-        if (c[si] < c[lRef])
-        {
-            const scalar exp = R.lhs()[slRef].exponent;
-            pf *= pow(max(c[lRef], 0.0), exp);
-            lRef = si;
-            slRef = s;
-        }
-        else
-        {
-            const scalar exp = R.lhs()[s].exponent;
-            pf *= pow(max(c[si], 0.0), exp);
-        }
-    }
-    cf = max(c[lRef], 0.0);
-
-    {
-        const scalar exp = R.lhs()[slRef].exponent;
-        if (exp < 1.0)
-        {
-            if (cf > small)
-            {
-                pf *= pow(cf, exp - 1.0);
-            }
-            else
-            {
-                pf = 0.0;
-            }
-        }
-        else
-        {
-            pf *= pow(cf, exp - 1.0);
-        }
-    }
-
-    label srRef = 0;
-    rRef = R.rhs()[srRef].index;
-
-    // Find the matrix element and element position for the rhs
-    pr = kr;
-    for (label s = 1; s < Nr; s++)
-    {
-        const label si = R.rhs()[s].index;
-        if (c[si] < c[rRef])
-        {
-            const scalar exp = R.rhs()[srRef].exponent;
-            pr *= pow(max(c[rRef], 0.0), exp);
-            rRef = si;
-            srRef = s;
-        }
-        else
-        {
-            const scalar exp = R.rhs()[s].exponent;
-            pr *= pow(max(c[si], 0.0), exp);
-        }
-    }
-    cr = max(c[rRef], 0.0);
-
-    {
-        const scalar exp = R.rhs()[srRef].exponent;
-        if (exp < 1.0)
-        {
-            if (cr > small)
-            {
-                pr *= pow(cr, exp - 1.0);
-            }
-            else
-            {
-                pr = 0.0;
-            }
-        }
-        else
-        {
-            pr *= pow(cr, exp - 1.0);
-        }
-    }
-
-    return pf*cf - pr*cr;
-}
-
-
 template<class ReactionThermo, class ThermoType>
 void Foam::StandardChemistryModel<ReactionThermo, ThermoType>::derivatives
 (
@@ -286,29 +155,29 @@ void Foam::StandardChemistryModel<ReactionThermo, ThermoType>::derivatives
 
     forAll(c_, i)
     {
-        c_[i] = max(c[i], 0.0);
+        c_[i] = max(c[i], 0);
     }
 
     omega(c_, T, p, dcdt);
 
     // Constant pressure
     // dT/dt = ...
-    scalar rho = 0.0;
-    scalar cSum = 0.0;
+    scalar rho = 0;
+    scalar cSum = 0;
     for (label i = 0; i < nSpecie_; i++)
     {
         const scalar W = specieThermo_[i].W();
         cSum += c_[i];
         rho += W*c_[i];
     }
-    scalar cp = 0.0;
+    scalar cp = 0;
     for (label i=0; i<nSpecie_; i++)
     {
         cp += c_[i]*specieThermo_[i].cp(p, T);
     }
     cp /= rho;
 
-    scalar dT = 0.0;
+    scalar dT = 0;
     for (label i = 0; i < nSpecie_; i++)
     {
         const scalar hi = specieThermo_[i].ha(p, T);
@@ -319,7 +188,7 @@ void Foam::StandardChemistryModel<ReactionThermo, ThermoType>::derivatives
     dcdt[nSpecie_] = -dT;
 
     // dp/dt = ...
-    dcdt[nSpecie_ + 1] = 0.0;
+    dcdt[nSpecie_ + 1] = 0;
 }
 
 
@@ -329,7 +198,7 @@ void Foam::StandardChemistryModel<ReactionThermo, ThermoType>::jacobian
     const scalar t,
     const scalarField& c,
     scalarField& dcdt,
-    scalarSquareMatrix& dfdc
+    scalarSquareMatrix& J
 ) const
 {
     const scalar T = c[nSpecie_];
@@ -337,131 +206,67 @@ void Foam::StandardChemistryModel<ReactionThermo, ThermoType>::jacobian
 
     forAll(c_, i)
     {
-        c_[i] = max(c[i], 0.0);
+        c_[i] = max(c[i], 0);
     }
 
-    dfdc = Zero;
-
-    // Length of the first argument must be nSpecie_
-    omega(c_, T, p, dcdt);
+    J = Zero;
+    dcdt = Zero;
 
+    // To compute the species derivatives of the temperature term,
+    // the enthalpies of the individual species is needed
+    scalarField hi(nSpecie_);
+    scalarField cpi(nSpecie_);
+    for (label i = 0; i < nSpecie_; i++)
+    {
+        hi[i] = specieThermo_[i].ha(p, T);
+        cpi[i] = specieThermo_[i].cp(p, T);
+    }
+    scalar omegaI = 0;
+    List<label> dummy;
     forAll(reactions_, ri)
     {
         const Reaction<ThermoType>& R = reactions_[ri];
-
-        const scalar kf0 = R.kf(p, T, c_);
-        const scalar kr0 = R.kr(kf0, p, T, c_);
-
-        forAll(R.lhs(), j)
-        {
-            const label sj = R.lhs()[j].index;
-            scalar kf = kf0;
-            forAll(R.lhs(), i)
-            {
-                const label si = R.lhs()[i].index;
-                const scalar el = R.lhs()[i].exponent;
-                if (i == j)
-                {
-                    if (el < 1.0)
-                    {
-                        if (c_[si] > small)
-                        {
-                            kf *= el*pow(c_[si], el - 1.0);
-                        }
-                        else
-                        {
-                            kf = 0.0;
-                        }
-                    }
-                    else
-                    {
-                        kf *= el*pow(c_[si], el - 1.0);
-                    }
-                }
-                else
-                {
-                    kf *= pow(c_[si], el);
-                }
-            }
-
-            forAll(R.lhs(), i)
-            {
-                const label si = R.lhs()[i].index;
-                const scalar sl = R.lhs()[i].stoichCoeff;
-                dfdc(si, sj) -= sl*kf;
-            }
-            forAll(R.rhs(), i)
-            {
-                const label si = R.rhs()[i].index;
-                const scalar sr = R.rhs()[i].stoichCoeff;
-                dfdc(si, sj) += sr*kf;
-            }
-        }
-
-        forAll(R.rhs(), j)
-        {
-            const label sj = R.rhs()[j].index;
-            scalar kr = kr0;
-            forAll(R.rhs(), i)
-            {
-                const label si = R.rhs()[i].index;
-                const scalar er = R.rhs()[i].exponent;
-                if (i == j)
-                {
-                    if (er < 1.0)
-                    {
-                        if (c_[si] > small)
-                        {
-                            kr *= er*pow(c_[si], er - 1.0);
-                        }
-                        else
-                        {
-                            kr = 0.0;
-                        }
-                    }
-                    else
-                    {
-                        kr *= er*pow(c_[si], er - 1.0);
-                    }
-                }
-                else
-                {
-                    kr *= pow(c_[si], er);
-                }
-            }
-
-            forAll(R.lhs(), i)
-            {
-                const label si = R.lhs()[i].index;
-                const scalar sl = R.lhs()[i].stoichCoeff;
-                dfdc(si, sj) += sl*kr;
-            }
-            forAll(R.rhs(), i)
-            {
-                const label si = R.rhs()[i].index;
-                const scalar sr = R.rhs()[i].stoichCoeff;
-                dfdc(si, sj) -= sr*kr;
-            }
-        }
+        scalar kfwd, kbwd;
+        R.dwdc(p, T, c_, J, dcdt, omegaI, kfwd, kbwd, false, dummy);
+        R.dwdT(p, T, c_, omegaI, kfwd, kbwd, J, false, dummy, nSpecie_);
     }
 
-    // Calculate the dcdT elements numerically
-    const scalar delta = 1.0e-3;
-
-    omega(c_, T + delta, p, dcdt_);
+    // The species derivatives of the temperature term are partially computed
+    // while computing dwdc, they are completed hereunder:
+    scalar cpMean = 0;
+    scalar dcpdTMean = 0;
     for (label i=0; i<nSpecie_; i++)
     {
-        dfdc(i, nSpecie_) = dcdt_[i];
+        cpMean += c_[i]*cpi[i]; // J/(m3.K)
+        dcpdTMean += c_[i]*specieThermo_[i].dcpdT(p, T);
     }
-
-    omega(c_, T - delta, p, dcdt_);
+    scalar dTdt = 0.0;
     for (label i=0; i<nSpecie_; i++)
     {
-        dfdc(i, nSpecie_) = 0.5*(dfdc(i, nSpecie_) - dcdt_[i])/delta;
+        dTdt += hi[i]*dcdt[i]; // J/(m3.s)
+    }
+    dTdt /= -cpMean; // K/s
+
+    for (label i = 0; i < nSpecie_; i++)
+    {
+        J(nSpecie_, i) = 0;
+        for (label j = 0; j < nSpecie_; j++)
+        {
+            J(nSpecie_, i) += hi[j]*J(j, i);
+        }
+        J(nSpecie_, i) += cpi[i]*dTdt; // J/(mol.s)
+        J(nSpecie_, i) /= -cpMean;    // K/s / (mol/m3)
     }
 
-    dfdc(nSpecie_, nSpecie_) = 0;
-    dfdc(nSpecie_ + 1, nSpecie_) = 0;
+    // ddT of dTdt
+    J(nSpecie_, nSpecie_) = 0;
+    for (label i = 0; i < nSpecie_; i++)
+    {
+        J(nSpecie_, nSpecie_) += cpi[i]*dcdt[i] + hi[i]*J(i, nSpecie_);
+    }
+    J(nSpecie_, nSpecie_) += dTdt*dcpdTMean;
+    J(nSpecie_, nSpecie_) /= -cpMean;
+    J(nSpecie_, nSpecie_) += dTdt/T;
 }
 
 
@@ -509,7 +314,7 @@ Foam::StandardChemistryModel<ReactionThermo, ThermoType>::tc() const
             const scalar Ti = T[celli];
             const scalar pi = p[celli];
 
-            scalar cSum = 0.0;
+            scalar cSum = 0;
 
             for (label i=0; i<nSpecie_; i++)
             {
@@ -521,7 +326,7 @@ Foam::StandardChemistryModel<ReactionThermo, ThermoType>::tc() const
             {
                 const Reaction<ThermoType>& R = reactions_[i];
 
-                omega(R, c_, Ti, pi, pf, cf, lRef, pr, cr, rRef);
+                R.omega(pi, Ti, c_, pf, cf, lRef, pr, cr, rRef);
 
                 forAll(R.rhs(), s)
                 {
@@ -557,7 +362,7 @@ Foam::StandardChemistryModel<ReactionThermo, ThermoType>::Qdot() const
                 false
             ),
             this->mesh_,
-            dimensionedScalar("zero", dimEnergy/dimVolume/dimTime, 0.0)
+            dimensionedScalar("zero", dimEnergy/dimVolume/dimTime, 0)
         )
     );
 
@@ -600,7 +405,7 @@ Foam::StandardChemistryModel<ReactionThermo, ThermoType>::calculateRR
                 IOobject::NO_WRITE
             ),
             this->mesh(),
-            dimensionedScalar("zero", dimMass/dimVolume/dimTime, 0.0)
+            dimensionedScalar("zero", dimMass/dimVolume/dimTime, 0)
         )
     );
 
@@ -628,7 +433,7 @@ Foam::StandardChemistryModel<ReactionThermo, ThermoType>::calculateRR
         }
 
         const Reaction<ThermoType>& R = reactions_[ri];
-        const scalar omegai = omega(R, c_, Ti, pi, pf, cf, lRef, pr, cr, rRef);
+        const scalar omegai = R.omega(pi, Ti, c_, pf, cf, lRef, pr, cr, rRef);
 
         forAll(R.lhs(), s)
         {
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/StandardChemistryModel/StandardChemistryModel.H b/src/thermophysicalModels/chemistryModel/chemistryModel/StandardChemistryModel/StandardChemistryModel.H
index 1a75627a0..c7f1b5005 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/StandardChemistryModel/StandardChemistryModel.H
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/StandardChemistryModel/StandardChemistryModel.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -163,22 +163,6 @@ public:
             scalarField& dcdt
         ) const;
 
-        //- Return the reaction rate for reaction r and the reference
-        //  species and charateristic times
-        virtual scalar omega
-        (
-            const Reaction<ThermoType>& r,
-            const scalarField& c,
-            const scalar T,
-            const scalar p,
-            scalar& pf,
-            scalar& cf,
-            label& lRef,
-            scalar& pr,
-            scalar& cr,
-            label& rRef
-        ) const;
-
 
         //- Return the reaction rate for iReaction and the reference
         //  species and charateristic times
@@ -254,7 +238,7 @@ public:
                 const scalar t,
                 const scalarField& c,
                 scalarField& dcdt,
-                scalarSquareMatrix& dfdc
+                scalarSquareMatrix& J
             ) const;
 
             virtual void solve
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/TDACChemistryModel.C b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/TDACChemistryModel.C
index 7caf8bf55..777fa225d 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/TDACChemistryModel.C
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/TDACChemistryModel.C
@@ -168,9 +168,9 @@ void Foam::TDACChemistryModel<ReactionThermo, ThermoType>::omega
         {
             const Reaction<ThermoType>& R = this->reactions_[i];
 
-            scalar omegai = omega
+            scalar omegai = R.omega
             (
-                R, c, T, p, pf, cf, lRef, pr, cr, rRef
+                p, T, c, pf, cf, lRef, pr, cr, rRef
             );
 
             forAll(R.lhs(), s)
@@ -232,17 +232,17 @@ Foam::scalar Foam::TDACChemistryModel<ReactionThermo, ThermoType>::omega
         if (c[si] < c[lRef])
         {
             const scalar exp = R.lhs()[slRef].exponent;
-            pf *= pow(max(c[lRef], 0.0), exp);
+            pf *= pow(max(c[lRef], 0), exp);
             lRef = si;
             slRef = s;
         }
         else
         {
             const scalar exp = R.lhs()[s].exponent;
-            pf *= pow(max(c[si], 0.0), exp);
+            pf *= pow(max(c[si], 0), exp);
         }
     }
-    cf = max(c[lRef], 0.0);
+    cf = max(c[lRef], 0);
 
     {
         const scalar exp = R.lhs()[slRef].exponent;
@@ -274,17 +274,17 @@ Foam::scalar Foam::TDACChemistryModel<ReactionThermo, ThermoType>::omega
         if (c[si] < c[rRef])
         {
             const scalar exp = R.rhs()[srRef].exponent;
-            pr *= pow(max(c[rRef], 0.0), exp);
+            pr *= pow(max(c[rRef], 0), exp);
             rRef = si;
             srRef = s;
         }
         else
         {
             const scalar exp = R.rhs()[s].exponent;
-            pr *= pow(max(c[si], 0.0), exp);
+            pr *= pow(max(c[si], 0), exp);
         }
     }
-    cr = max(c[rRef], 0.0);
+    cr = max(c[rRef], 0);
 
     {
         const scalar exp = R.rhs()[srRef].exponent;
@@ -334,14 +334,14 @@ void Foam::TDACChemistryModel<ReactionThermo, ThermoType>::derivatives
         // efficiencies
         for (label i=0; i<NsDAC_; i++)
         {
-            this->c_[simplifiedToCompleteIndex_[i]] = max(c[i], 0.0);
+            this->c_[simplifiedToCompleteIndex_[i]] = max(c[i], 0);
         }
     }
     else
     {
         for (label i=0; i<this->nSpecie(); i++)
         {
-            this->c_[i] = max(c[i], 0.0);
+            this->c_[i] = max(c[i], 0);
         }
     }
 
@@ -398,7 +398,8 @@ void Foam::TDACChemistryModel<ReactionThermo, ThermoType>::jacobian
 (
     const scalar t,
     const scalarField& c,
-    scalarSquareMatrix& dfdc
+    scalarField& dcdt,
+    scalarSquareMatrix& J
 ) const
 {
     const bool reduced = mechRed_->active();
@@ -416,183 +417,119 @@ void Foam::TDACChemistryModel<ReactionThermo, ThermoType>::jacobian
         this->c_ = completeC_;
         for (label i=0; i<NsDAC_; i++)
         {
-            this->c_[simplifiedToCompleteIndex_[i]] = max(c[i], 0.0);
+            this->c_[simplifiedToCompleteIndex_[i]] = max(c[i], 0);
         }
     }
     else
     {
         forAll(this->c_, i)
         {
-            this->c_[i] = max(c[i], 0.0);
+            this->c_[i] = max(c[i], 0);
         }
     }
 
-    dfdc = Zero;
+    J = Zero;
+    dcdt = Zero;
+    scalarField hi(this->c_.size());
+    scalarField cpi(this->c_.size());
+    forAll(hi, i)
+    {
+        hi[i] = this->specieThermo_[i].ha(p, T);
+        cpi[i] = this->specieThermo_[i].cp(p, T);
+    }
+
+    scalar omegaI = 0;
 
     forAll(this->reactions_, ri)
     {
         if (!reactionsDisabled_[ri])
         {
             const Reaction<ThermoType>& R = this->reactions_[ri];
-
-            const scalar kf0 = R.kf(p, T, this->c_);
-            const scalar kr0 = R.kr(kf0, p, T, this->c_);
-
-            forAll(R.lhs(), j)
-            {
-                label sj = R.lhs()[j].index;
-                if (reduced)
-                {
-                    sj = completeToSimplifiedIndex_[sj];
-                }
-                scalar kf = kf0;
-                forAll(R.lhs(), i)
-                {
-                    const label si = R.lhs()[i].index;
-                    const scalar el = R.lhs()[i].exponent;
-                    if (i == j)
-                    {
-                        if (el < 1)
-                        {
-                            if (this->c_[si] > small)
-                            {
-                                kf *= el*pow(this->c_[si], el - 1);
-                            }
-                            else
-                            {
-                                kf = 0;
-                            }
-                        }
-                        else
-                        {
-                            kf *= el*pow(this->c_[si], el - 1);
-                        }
-                    }
-                    else
-                    {
-                        kf *= pow(this->c_[si], el);
-                    }
-                }
-
-                forAll(R.lhs(), i)
-                {
-                    label si = R.lhs()[i].index;
-                    if (reduced)
-                    {
-                        si = completeToSimplifiedIndex_[si];
-                    }
-                    const scalar sl = R.lhs()[i].stoichCoeff;
-                    dfdc(si, sj) -= sl*kf;
-                }
-                forAll(R.rhs(), i)
-                {
-                    label si = R.rhs()[i].index;
-                    if (reduced)
-                    {
-                        si = completeToSimplifiedIndex_[si];
-                    }
-                    const scalar sr = R.rhs()[i].stoichCoeff;
-                    dfdc(si, sj) += sr*kf;
-                }
-            }
-
-            forAll(R.rhs(), j)
-            {
-                label sj = R.rhs()[j].index;
-                if (reduced)
-                {
-                    sj = completeToSimplifiedIndex_[sj];
-                }
-                scalar kr = kr0;
-                forAll(R.rhs(), i)
-                {
-                    const label si = R.rhs()[i].index;
-                    const scalar er = R.rhs()[i].exponent;
-                    if (i == j)
-                    {
-                        if (er < 1)
-                        {
-                            if (this->c_[si] > small)
-                            {
-                                kr *= er*pow(this->c_[si], er - 1);
-                            }
-                            else
-                            {
-                                kr = 0;
-                            }
-                        }
-                        else
-                        {
-                            kr *= er*pow(this->c_[si], er - 1);
-                        }
-                    }
-                    else
-                    {
-                        kr *= pow(this->c_[si], er);
-                    }
-                }
-
-                forAll(R.lhs(), i)
-                {
-                    label si = R.lhs()[i].index;
-                    if (reduced)
-                    {
-                        si = completeToSimplifiedIndex_[si];
-                    }
-                    const scalar sl = R.lhs()[i].stoichCoeff;
-                    dfdc(si, sj) += sl*kr;
-                }
-                forAll(R.rhs(), i)
-                {
-                    label si = R.rhs()[i].index;
-                    if (reduced)
-                    {
-                        si = completeToSimplifiedIndex_[si];
-                    }
-                    const scalar sr = R.rhs()[i].stoichCoeff;
-                    dfdc(si, sj) -= sr*kr;
-                }
-            }
+            scalar kfwd, kbwd;
+            R.dwdc
+            (
+                p,
+                T,
+                this->c_,
+                J,
+                dcdt,
+                omegaI,
+                kfwd,
+                kbwd,
+                reduced,
+                completeToSimplifiedIndex_
+            );
+            R.dwdT
+            (
+                p,
+                T,
+                this->c_,
+                omegaI,
+                kfwd,
+                kbwd,
+                J,
+                reduced,
+                completeToSimplifiedIndex_,
+                this->nSpecie_
+            );
         }
     }
 
-    // Calculate the dcdT elements numerically
-    const scalar delta = 1e-3;
-
-    omega(this->c_, T + delta, p, this->dcdt_);
-    for (label i=0; i<this->nSpecie_; i++)
+    // The species derivatives of the temperature term are partially computed
+    // while computing dwdc, they are completed hereunder:
+    scalar cpMean = 0;
+    scalar dcpdTMean = 0;
+    forAll(this->c_, i)
     {
-        dfdc(i, this->nSpecie_) = this->dcdt_[i];
+        cpMean += this->c_[i]*cpi[i]; // J/(m3.K)
+        // Already multiplied by rho
+        dcpdTMean += this->c_[i]*this->specieThermo_[i].dcpdT(p, T);
     }
 
-    omega(this->c_, T - delta, p, this->dcdt_);
-    for (label i=0; i<this->nSpecie_; i++)
+    scalar dTdt = 0;
+    forAll(hi, i)
     {
-        dfdc(i, this->nSpecie_) =
-            0.5*(dfdc(i, this->nSpecie_) - this->dcdt_[i])/delta;
+        if (reduced)
+        {
+            const label si = completeToSimplifiedIndex_[i];
+            if (si != -1)
+            {
+                dTdt += hi[i]*dcdt[si]; // J/(m3.s)
+            }
+        }
+        else
+        {
+            dTdt += hi[i]*dcdt[i]; // J/(m3.s)
+        }
+    }
+    dTdt /= -cpMean; // K/s
+    dcdt[this->nSpecie_] = dTdt;
+
+    for (label i = 0; i < this->nSpecie_; i++)
+    {
+        J(this->nSpecie_, i) = 0;
+        for (label j = 0; j < this->nSpecie_; j++)
+        {
+            const label sj = reduced ? simplifiedToCompleteIndex_[j] : j;
+            J(this->nSpecie_, i) += hi[sj]*J(j, i);
+        }
+        const label si = reduced ? simplifiedToCompleteIndex_[i] : i;
+        J(this->nSpecie_, i) += cpi[si]*dTdt; // J/(mol.s)
+        J(this->nSpecie_, i) /= -cpMean;    // K/s / (mol/m3)
     }
 
-    dfdc(this->nSpecie_, this->nSpecie_) = 0;
-    dfdc(this->nSpecie_ + 1, this->nSpecie_) = 0;
-}
-
-
-template<class ReactionThermo, class ThermoType>
-void Foam::TDACChemistryModel<ReactionThermo, ThermoType>::jacobian
-(
-    const scalar t,
-    const scalarField& c,
-    scalarField& dcdt,
-    scalarSquareMatrix& dfdc
-) const
-{
-    jacobian(t, c, dfdc);
-
-    const scalar T = c[this->nSpecie_];
-    const scalar p = c[this->nSpecie_ + 1];
-
-    // Note: Uses the c_ field initialized by the call to jacobian above
-    omega(this->c_, T, p, dcdt);
+    // ddT of dTdt
+    J(this->nSpecie_, this->nSpecie_) = 0;
+    for (label i = 0; i < this->nSpecie_; i++)
+    {
+        const label si = reduced ? simplifiedToCompleteIndex_[i] : i;
+        J(this->nSpecie_, this->nSpecie_) +=
+            cpi[si]*dcdt[i]
+          + hi[si]*J(i, this->nSpecie_);
+    }
+    J(this->nSpecie_, this->nSpecie_) += dTdt*dcpdTMean;
+    J(this->nSpecie_, this->nSpecie_) /= -cpMean;
+    J(this->nSpecie_, this->nSpecie_) += dTdt/T;
 }
 
 
@@ -915,7 +852,7 @@ setTabulationResultsAdd
     const label celli
 )
 {
-    tabulationResults_[celli] = 0.0;
+    tabulationResults_[celli] = 0;
 }
 
 
@@ -923,7 +860,7 @@ template<class ReactionThermo, class ThermoType>
 void Foam::TDACChemistryModel<ReactionThermo, ThermoType>::
 setTabulationResultsGrow(const label celli)
 {
-    tabulationResults_[celli] = 1.0;
+    tabulationResults_[celli] = 1;
 }
 
 
@@ -934,7 +871,7 @@ setTabulationResultsRetrieve
     const label celli
 )
 {
-    tabulationResults_[celli] = 2.0;
+    tabulationResults_[celli] = 2;
 }
 
 
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/TDACChemistryModel.H b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/TDACChemistryModel.H
index 87cc80eb6..f2cc41a68 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/TDACChemistryModel.H
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/TDACChemistryModel.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2016-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2016-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -221,20 +221,12 @@ public:
                 scalarField& dcdt
             ) const;
 
-            //- Pure jacobian function for tabulation
-            void jacobian
-            (
-                const scalar t,
-                const scalarField& c,
-                scalarSquareMatrix& dfdc
-            ) const;
-
             virtual void jacobian
             (
                 const scalar t,
                 const scalarField& c,
                 scalarField& dcdt,
-                scalarSquareMatrix& dfdc
+                scalarSquareMatrix& J
             ) const;
 
             virtual void solve
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/ISAT/ISAT.C b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/ISAT/ISAT.C
index f8d957707..9b3d86233 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/ISAT/ISAT.C
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/ISAT/ISAT.C
@@ -65,7 +65,7 @@ Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::ISAT
     (
         this->coeffsDict_.lookupOrDefault
         (
-            "minBalanceThreshold",0.1*chemisTree_.maxNLeafs()
+            "minBalanceThreshold", 0.1*chemisTree_.maxNLeafs()
         )
     ),
     MRURetrieve_(this->coeffsDict_.lookupOrDefault("MRURetrieve", false)),
@@ -169,17 +169,8 @@ void Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::addToMRU
         {
             if (MRUList_.size() == maxMRUSize_)
             {
-                if (iter() == MRUList_.last())
-                {
-                    MRUList_.remove(iter);
-                    MRUList_.insert(phi0);
-                }
-                else
-                {
-                    FatalErrorInFunction
-                        << "wrong MRUList construction"
-                        << exit(FatalError);
-                }
+                MRUList_.remove(iter);
+                MRUList_.insert(phi0);
             }
             else
             {
@@ -238,13 +229,13 @@ void Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::calcNewC
 
                 // As we use an approximation of A, Rphiq should be checked for
                 // negative values
-                Rphiq[i] = max(0.0,Rphiq[i]);
+                Rphiq[i] = max(0, Rphiq[i]);
             }
             // The species is not active A(i, j) = I(i, j)
             else
             {
                 Rphiq[i] += dphi[i];
-                Rphiq[i] = max(0.0,Rphiq[i]);
+                Rphiq[i] = max(0, Rphiq[i]);
             }
         }
         else // Mechanism reduction is not active
@@ -255,7 +246,7 @@ void Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::calcNewC
             }
             // As we use a first order gradient matrix, Rphiq should be checked
             // for negative values
-            Rphiq[i] = max(0.0,Rphiq[i]);
+            Rphiq[i] = max(0, Rphiq[i]);
         }
     }
 }
@@ -322,6 +313,9 @@ Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::cleanAndBalance()
         }
         x = xtmp;
     }
+
+    MRUList_.clear();
+
     // Check if the tree should be balanced according to criterion:
     //  -the depth of the tree bigger than a*log2(size), log2(size) being the
     //      ideal depth (e.g. 4 leafs can be stored in a tree of depth 2)
@@ -333,7 +327,6 @@ Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::cleanAndBalance()
     )
     {
         chemisTree_.balance();
-        MRUList_.clear();
         treeModified = true;
     }
 
@@ -365,7 +358,7 @@ void Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::computeA
         Rcq[i] = rhoi*Rphiq[s2c]/this->chemistry_.specieThermo()[s2c].W();
     }
     Rcq[speciesNumber] = Rphiq[Rphiq.size() - nAdditionalEqns_];
-    Rcq[speciesNumber+1] = Rphiq[Rphiq.size() - nAdditionalEqns_ + 1];
+    Rcq[speciesNumber + 1] = Rphiq[Rphiq.size() - nAdditionalEqns_ + 1];
     if (this->variableTimeStep())
     {
         Rcq[speciesNumber + 2] = Rphiq[Rphiq.size() - nAdditionalEqns_ + 2];
@@ -373,15 +366,15 @@ void Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::computeA
 
     // Aaa is computed implicitly,
     // A is given by A = C(psi0, t0+dt), where C is obtained through solving
-    // d/dt C(psi0,t) = J(psi(t))C(psi0,t)
+    // d/dt C(psi0, t) = J(psi(t))C(psi0, t)
     // If we solve it implicitly:
-    // (C(psi0, t0+dt) - C(psi0,t0))/dt = J(psi(t0+dt))C(psi0,t0+dt)
+    // (C(psi0, t0+dt) - C(psi0, t0))/dt = J(psi(t0+dt))C(psi0, t0+dt)
     // The Jacobian is thus computed according to the mapping
     // C(psi0,t0+dt)*(I-dt*J(psi(t0+dt))) = C(psi0, t0)
     // A = C(psi0,t0)/(I-dt*J(psi(t0+dt)))
     // where C(psi0,t0) = I
-
-    this->chemistry_.jacobian(runTime_.value(), Rcq, A);
+    scalarField dcdt(speciesNumber + 2, Zero);
+    this->chemistry_.jacobian(runTime_.value(), Rcq, dcdt, A);
 
     // The jacobian is computed according to the molar concentration
     // the following conversion allows the code to use A with mass fraction
@@ -410,21 +403,48 @@ void Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::computeA
         // Columns for pressure and temperature
         A(i, speciesNumber) *=
             -dt*this->chemistry_.specieThermo()[si].W()/rhoi;
-        A(i, speciesNumber+1) *=
+        A(i, speciesNumber + 1) *=
             -dt*this->chemistry_.specieThermo()[si].W()/rhoi;
     }
 
-    // For temperature and pressure, only unity on the diagonal
-    A(speciesNumber, speciesNumber) = 1;
-    A(speciesNumber + 1, speciesNumber + 1) = 1;
+    // For the temperature and pressure lines, ddc(dTdt)
+    // should be converted in ddY(dTdt)
+    for (label i=0; i<speciesNumber; i++)
+    {
+        label si = i;
+        if (mechRedActive)
+        {
+            si = this->chemistry_.simplifiedToCompleteIndex()[i];
+        }
+
+        A(speciesNumber, i) *=
+            -dt*rhoi/this->chemistry_.specieThermo()[si].W();
+        A(speciesNumber + 1, i) *=
+            -dt*rhoi/this->chemistry_.specieThermo()[si].W();
+    }
+
+    A(speciesNumber, speciesNumber) = -dt*A(speciesNumber, speciesNumber) + 1;
+
+    A(speciesNumber + 1, speciesNumber + 1) =
+        -dt*A(speciesNumber + 1, speciesNumber + 1) + 1;
+
     if (this->variableTimeStep())
     {
-        A[speciesNumber + 2][speciesNumber + 2] = 1;
+        A(speciesNumber + 2, speciesNumber + 2) = 1;
     }
 
     // Inverse of (I-dt*J(psi(t0+dt)))
     LUscalarMatrix LUA(A);
     LUA.inv(A);
+
+    // After inversion, lines of p and T are set to 0 except diagonal.  This
+    // avoid skewness of the ellipsoid of accuracy and potential issues in the
+    // binary tree.
+    for (label i=0; i<speciesNumber; i++)
+    {
+        A(speciesNumber, i) = 0;
+        A(speciesNumber + 1, i) = 0;
+    }
 }
 
 
@@ -498,7 +518,7 @@ bool Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::retrieve
         }
         lastSearch_->lastTimeUsed() = this->chemistry_.timeSteps();
         addToMRU(phi0);
-        calcNewC(phi0,phiq, Rphiq);
+        calcNewC(phi0, phiq, Rphiq);
         nRetrieved_++;
         return true;
     }
@@ -525,11 +545,12 @@ Foam::label Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::add
     // option is on, the code first tries to grow the point hold by lastSearch_
     if (lastSearch_ && growPoints_)
     {
-        if (grow(lastSearch_,phiq, Rphiq))
+        if (grow(lastSearch_, phiq, Rphiq))
         {
             nGrowth_++;
             growthOrAddFlag = 0;
-            // the structure of the tree is not modified, return false
+            addToMRU(lastSearch_);
+            //the structure of the tree is not modified, return false
             return growthOrAddFlag;
         }
     }
@@ -605,7 +626,10 @@ Foam::label Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::add
         scaleFactor_.size(),
         lastSearch_ // lastSearch_ may be nullptr (handled by binaryTree)
     );
-
+    if (lastSearch_ != nullptr)
+    {
+        addToMRU(lastSearch_);
+    }
     nAdd_++;
 
     return growthOrAddFlag;
diff --git a/src/thermophysicalModels/radiation/derivedFvPatchFields/radiationCoupledBase/radiationCoupledBase.C b/src/thermophysicalModels/radiation/derivedFvPatchFields/radiationCoupledBase/radiationCoupledBase.C
index bc41bce27..3b410cd38 100644
--- a/src/thermophysicalModels/radiation/derivedFvPatchFields/radiationCoupledBase/radiationCoupledBase.C
+++ b/src/thermophysicalModels/radiation/derivedFvPatchFields/radiationCoupledBase/radiationCoupledBase.C
@@ -28,7 +28,7 @@ License
 #include "mappedPatchBase.H"
 #include "fvPatchFieldMapper.H"
 #include "radiationModel.H"
-#include "noRadiation.H"
+#include "opaqueSolid.H"
 #include "absorptionEmissionModel.H"
 
 // * * * * * * * * * * * * * Static Member Data  * * * * * * * * * * * * * * //
@@ -151,29 +151,20 @@ Foam::scalarField Foam::radiationCoupledBase::emissivity() const
 
             const polyMesh& nbrMesh = mpp.sampleMesh();
 
-            const radiation::radiationModel& radiation =
-                nbrMesh.lookupObject<radiation::radiationModel>
+            const radiation::opaqueSolid& radiation =
+                nbrMesh.lookupObject<radiation::opaqueSolid>
                 (
                     "radiationProperties"
                 );
 
-            if (isType<radiation::noRadiation>(radiation))
-            {
-                FatalErrorInFunction
-                    << "No radiation model defined for region "
-                    << nbrMesh.name() << nl
-                    << "    required for option 'solidRadiation' "
-                       "of patchField type " << type()
-                    << " on patch " << patch_.name()
-                    << abort(FatalError);
-            }
-
             const fvMesh& nbrFvMesh = refCast<const fvMesh>(nbrMesh);
 
             const fvPatch& nbrPatch =
                 nbrFvMesh.boundary()[mpp.samplePolyPatch().index()];
 
-
+            // NOTE: for an opaqueSolid the absorptionEmission model returns the
+            // emissivity of the surface rather than the emission coefficient
+            // and the input specification MUST correspond to this.
             scalarField emissivity
             (
                 radiation.absorptionEmission().e()().boundaryField()
diff --git a/src/thermophysicalModels/reactionThermo/chemistryReaders/chemkinReader/chemkinLexer.L b/src/thermophysicalModels/reactionThermo/chemistryReaders/chemkinReader/chemkinLexer.L
index 552d759d7..73cc72a70 100644
--- a/src/thermophysicalModels/reactionThermo/chemistryReaders/chemkinReader/chemkinLexer.L
+++ b/src/thermophysicalModels/reactionThermo/chemistryReaders/chemkinReader/chemkinLexer.L
@@ -651,6 +651,18 @@ bool finishReaction = false;
     }
 
 <readThermoSpecieName>{end} {
+        Reaction<gasHThermoPhysics>::TlowDefault = max
+        (
+            Reaction<gasHThermoPhysics>::TlowDefault,
+            currentLowT
+        );
+
+        Reaction<gasHThermoPhysics>::ThighDefault = min
+        (
+            Reaction<gasHThermoPhysics>::ThighDefault,
+            currentHighT
+        );
+
         BEGIN(INITIAL);
     }
 
diff --git a/src/thermophysicalModels/reactionThermo/chemistryReaders/chemkinReader/chemkinReader.C b/src/thermophysicalModels/reactionThermo/chemistryReaders/chemkinReader/chemkinReader.C
index 029956ba8..710192740 100644
--- a/src/thermophysicalModels/reactionThermo/chemistryReaders/chemkinReader/chemkinReader.C
+++ b/src/thermophysicalModels/reactionThermo/chemistryReaders/chemkinReader/chemkinReader.C
@@ -26,6 +26,7 @@ License
 #include "chemkinReader.H"
 #include <fstream>
 #include "atomicWeights.H"
+#include "ReactionProxy.H"
 #include "IrreversibleReaction.H"
 #include "ReversibleReaction.H"
 #include "NonEquilibriumReversibleReaction.H"
@@ -182,7 +183,7 @@ void Foam::chemkinReader::addReactionType
                 new IrreversibleReaction
                 <Reaction, gasHThermoPhysics, ReactionRateType>
                 (
-                    Reaction<gasHThermoPhysics>
+                    ReactionProxy<gasHThermoPhysics>
                     (
                         speciesTable_,
                         lhs.shrink(),
@@ -202,7 +203,7 @@ void Foam::chemkinReader::addReactionType
                 new ReversibleReaction
                 <Reaction, gasHThermoPhysics, ReactionRateType>
                 (
-                    Reaction<gasHThermoPhysics>
+                    ReactionProxy<gasHThermoPhysics>
                     (
                         speciesTable_,
                         lhs.shrink(),
@@ -490,7 +491,7 @@ void Foam::chemkinReader::addReaction
                     new NonEquilibriumReversibleReaction
                         <Reaction, gasHThermoPhysics, ArrheniusReactionRate>
                     (
-                        Reaction<gasHThermoPhysics>
+                        ReactionProxy<gasHThermoPhysics>
                         (
                             speciesTable_,
                             lhs.shrink(),
@@ -546,7 +547,7 @@ void Foam::chemkinReader::addReaction
                         thirdBodyArrheniusReactionRate
                     >
                     (
-                        Reaction<gasHThermoPhysics>
+                        ReactionProxy<gasHThermoPhysics>
                         (
                             speciesTable_,
                             lhs.shrink(),
@@ -646,9 +647,13 @@ void Foam::chemkinReader::addReaction
                 reactions_.append
                 (
                     new NonEquilibriumReversibleReaction
-                        <Reaction, gasHThermoPhysics, LandauTellerReactionRate>
+                    <
+                        Reaction,
+                        gasHThermoPhysics,
+                        LandauTellerReactionRate
+                    >
                     (
-                        Reaction<gasHThermoPhysics>
+                        ReactionProxy<gasHThermoPhysics>
                         (
                             speciesTable_,
                             lhs.shrink(),
@@ -780,6 +785,9 @@ void Foam::chemkinReader::read
     const fileName& transportFileName
 )
 {
+    Reaction<gasHThermoPhysics>::TlowDefault = 0;
+    Reaction<gasHThermoPhysics>::ThighDefault = great;
+
     transportDict_.read(IFstream(transportFileName)());
 
     if (thermoFileName != fileName::null)
diff --git a/src/thermophysicalModels/solidSpecie/reaction/Reactions/solidReaction/solidReaction.H b/src/thermophysicalModels/solidSpecie/reaction/Reactions/solidReaction/solidReaction.H
index 84e4f61f2..a98126696 100644
--- a/src/thermophysicalModels/solidSpecie/reaction/Reactions/solidReaction/solidReaction.H
+++ b/src/thermophysicalModels/solidSpecie/reaction/Reactions/solidReaction/solidReaction.H
@@ -59,7 +59,7 @@ inline Ostream& operator<<(Ostream&, const solidReaction<ReactionThermo>&);
 template<class ReactionThermo>
 class solidReaction
 :
-    public  Reaction<ReactionThermo>
+    public Reaction<ReactionThermo>
 {
 
 private:
@@ -129,28 +129,6 @@ public:
         );
 
 
-        //- Construct and return a clone
-        virtual autoPtr<Reaction<ReactionThermo>> clone() const
-        {
-            return autoPtr<Reaction<ReactionThermo>>
-            (
-                new solidReaction<ReactionThermo>(*this)
-            );
-        }
-
-        //- Construct and return a clone with new speciesTable
-        virtual autoPtr<Reaction<ReactionThermo>> clone
-        (
-            const speciesTable& species
-        ) const
-        {
-            return autoPtr<Reaction<ReactionThermo>>
-            (
-                new solidReaction<ReactionThermo>(*this, species)
-            );
-        }
-
-
     //- Destructor
     virtual ~solidReaction()
     {}
@@ -160,12 +138,13 @@ public:
 
         // Access
 
-
-            // - Access to gas components of the reaction
+            //- Access to the gas components of the left hand side
             virtual const List<specieCoeffs>& grhs() const;
+
+            //- Access to the gas components of the right hand side
             virtual const List<specieCoeffs>& glhs() const;
 
-            // - Access to gas specie list
+            //- Access to gas specie list
             virtual const speciesTable& gasSpecies() const;
 
 
diff --git a/src/thermophysicalModels/solidSpecie/reaction/reactionRate/solidArrheniusReactionRate/solidArrheniusReactionRate.H b/src/thermophysicalModels/solidSpecie/reaction/reactionRate/solidArrheniusReactionRate/solidArrheniusReactionRate.H
index c5642e627..86ff4639d 100644
--- a/src/thermophysicalModels/solidSpecie/reaction/reactionRate/solidArrheniusReactionRate/solidArrheniusReactionRate.H
+++ b/src/thermophysicalModels/solidSpecie/reaction/reactionRate/solidArrheniusReactionRate/solidArrheniusReactionRate.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -104,6 +104,37 @@ public:
             const scalarField& c
         ) const;
 
+        inline scalar ddT
+        (
+            const scalar p,
+            const scalar T,
+            const scalarField& c
+        ) const;
+
+        //- Third-body efficiencies (beta = 1-alpha)
+        //  non-empty only for third-body reactions
+        //  with enhanced molecularity (alpha != 1)
+        inline const List<Tuple2<label, scalar>>& beta() const;
+
+        //- Species concentration derivative of the pressure dependent term
+        //  By default this value is 1 as it multiplies the third-body term
+        inline void dcidc
+        (
+            const scalar p,
+            const scalar T,
+            const scalarField& c,
+            scalarField& dcidc
+        ) const;
+
+        //- Temperature derivative of the pressure dependent term
+        //  By default this value is 0 since ddT of molecularity is approx.0
+        inline scalar dcidT
+        (
+            const scalar p,
+            const scalar T,
+            const scalarField& c
+        ) const;
+
 
          //- Write to stream
         inline void write(Ostream& os) const;
diff --git a/src/thermophysicalModels/solidSpecie/reaction/reactionRate/solidArrheniusReactionRate/solidArrheniusReactionRateI.H b/src/thermophysicalModels/solidSpecie/reaction/reactionRate/solidArrheniusReactionRate/solidArrheniusReactionRateI.H
index cd64c2f36..bf32298db 100644
--- a/src/thermophysicalModels/solidSpecie/reaction/reactionRate/solidArrheniusReactionRate/solidArrheniusReactionRateI.H
+++ b/src/thermophysicalModels/solidSpecie/reaction/reactionRate/solidArrheniusReactionRate/solidArrheniusReactionRateI.H
@@ -60,18 +60,60 @@ inline Foam::scalar Foam::solidArrheniusReactionRate::operator()
     const scalarField&
 ) const
 {
-    scalar ak = A_;
-
     if (T < Tcrit_)
     {
-        ak *= 0.0;
+        return 0;
     }
     else
     {
-        ak *= exp(-Ta_/T);
+        return A_*exp(-Ta_/T);
     }
+}
 
-    return ak;
+
+inline Foam::scalar Foam::solidArrheniusReactionRate::ddT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField&
+) const
+{
+    if (T < Tcrit_)
+    {
+        return 0;
+    }
+    else
+    {
+        return A_*exp(-Ta_/T)*Ta_/sqr(T);
+    }
+}
+
+
+inline const Foam::List<Foam::Tuple2<Foam::label, Foam::scalar>>&
+Foam::solidArrheniusReactionRate::beta() const
+{
+    return NullObjectRef<List<Tuple2<label, scalar>>>();
+}
+
+
+inline void Foam::solidArrheniusReactionRate::dcidc
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c,
+    scalarField& dcidc
+) const
+{}
+
+
+inline Foam::scalar Foam::solidArrheniusReactionRate::dcidT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c
+) const
+{
+    return 0;
 }
 
 
diff --git a/src/thermophysicalModels/specie/reaction/Reactions/IrreversibleReaction/IrreversibleReaction.C b/src/thermophysicalModels/specie/reaction/Reactions/IrreversibleReaction/IrreversibleReaction.C
index 32a301fab..25b843724 100644
--- a/src/thermophysicalModels/specie/reaction/Reactions/IrreversibleReaction/IrreversibleReaction.C
+++ b/src/thermophysicalModels/specie/reaction/Reactions/IrreversibleReaction/IrreversibleReaction.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -106,6 +106,160 @@ Foam::scalar Foam::IrreversibleReaction
 }
 
 
+template
+<
+    template<class> class ReactionType,
+    class ReactionThermo,
+    class ReactionRate
+>
+Foam::scalar Foam::IrreversibleReaction
+<
+    ReactionType,
+    ReactionThermo,
+    ReactionRate
+>::kr
+(
+    const scalar kfwd,
+    const scalar p,
+    const scalar T,
+    const scalarField& c
+) const
+{
+    return 0;
+}
+
+
+template
+<
+    template<class> class ReactionType,
+    class ReactionThermo,
+    class ReactionRate
+>
+Foam::scalar Foam::IrreversibleReaction
+<
+    ReactionType,
+    ReactionThermo,
+    ReactionRate
+>::kr
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c
+) const
+{
+    return 0;
+}
+
+
+template
+<
+    template<class> class ReactionType,
+    class ReactionThermo,
+    class ReactionRate
+>
+Foam::scalar Foam::IrreversibleReaction
+<
+    ReactionType,
+    ReactionThermo,
+    ReactionRate
+>::dkfdT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c
+) const
+{
+    return k_.ddT(p, T, c);
+}
+
+
+template
+<
+    template<class> class ReactionType,
+    class ReactionThermo,
+    class ReactionRate
+>
+Foam::scalar Foam::IrreversibleReaction
+<
+    ReactionType,
+    ReactionThermo,
+    ReactionRate
+>::dkrdT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c,
+    const scalar dkfdT,
+    const scalar kr
+) const
+{
+    return 0;
+}
+
+
+template
+<
+    template<class> class ReactionType,
+    class ReactionThermo,
+    class ReactionRate
+>
+const Foam::List<Foam::Tuple2<Foam::label, Foam::scalar>>&
+Foam::IrreversibleReaction
+<
+    ReactionType,
+    ReactionThermo,
+    ReactionRate
+>::beta() const
+{
+    return k_.beta();
+}
+
+
+template
+<
+    template<class> class ReactionType,
+    class ReactionThermo,
+    class ReactionRate
+>
+void Foam::IrreversibleReaction
+<
+    ReactionType,
+    ReactionThermo,
+    ReactionRate
+>::dcidc
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c,
+    scalarField& dcidc
+) const
+{
+    k_.dcidc(p, T, c, dcidc);
+}
+
+
+template
+<
+    template<class> class ReactionType,
+    class ReactionThermo,
+    class ReactionRate
+>
+Foam::scalar Foam::IrreversibleReaction
+<
+    ReactionType,
+    ReactionThermo,
+    ReactionRate
+>::dcidT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c
+) const
+{
+    return k_.dcidT(p, T, c);
+}
+
+
 template
 <
     template<class> class ReactionType,
diff --git a/src/thermophysicalModels/specie/reaction/Reactions/IrreversibleReaction/IrreversibleReaction.H b/src/thermophysicalModels/specie/reaction/Reactions/IrreversibleReaction/IrreversibleReaction.H
index 3c3aa7ead..12771b748 100644
--- a/src/thermophysicalModels/specie/reaction/Reactions/IrreversibleReaction/IrreversibleReaction.H
+++ b/src/thermophysicalModels/specie/reaction/Reactions/IrreversibleReaction/IrreversibleReaction.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -25,8 +25,7 @@ Class
     Foam::IrreversibleReaction
 
 Description
-    Simple extension of Reaction to handle reversible reactions using
-    equilibrium thermodynamics.
+    Simple extension of Reaction to handle irreversible reactions
 
 SourceFiles
     IrreversibleReaction.C
@@ -164,6 +163,71 @@ public:
                 const scalarField& c
             ) const;
 
+            //- Reverse rate constant from the given forward rate constant
+            //  Returns 0
+            virtual scalar kr
+            (
+                const scalar kfwd,
+                const scalar p,
+                const scalar T,
+                const scalarField& c
+            ) const;
+
+            //- Reverse rate constant
+            //  Returns 0
+            virtual scalar kr
+            (
+                const scalar p,
+                const scalar T,
+                const scalarField& c
+            ) const;
+
+
+        // IrreversibleReaction Jacobian functions
+
+            //- Temperature derivative of forward rate
+            virtual scalar dkfdT
+            (
+                const scalar p,
+                const scalar T,
+                const scalarField& c
+            ) const;
+
+            //- Temperature derivative of reverse rate
+            //  Returns 0
+            virtual scalar dkrdT
+            (
+                const scalar p,
+                const scalar T,
+                const scalarField& c,
+                const scalar dkfdT,
+                const scalar kr
+            ) const;
+
+            //- Third-body efficiencies (beta = 1-alpha)
+            //  non-empty only for third-body reactions
+            //  with enhanced molecularity (alpha != 1)
+            virtual const List<Tuple2<label, scalar>>& beta() const;
+
+            //- Species concentration derivative of the pressure dependent term
+            //  By default this value is 1 as it multiplies the third-body term
+            virtual void dcidc
+            (
+                const scalar p,
+                const scalar T,
+                const scalarField& c,
+                scalarField& dcidc
+            ) const;
+
+            //- Temperature derivative of the pressure dependent term
+            //  By default this value is 0 since ddT of molecularity is approx.0
+            virtual scalar dcidT
+            (
+                const scalar p,
+                const scalar T,
+                const scalarField& c
+            ) const;
+
 
         //- Write
         virtual void write(Ostream&) const;
diff --git a/src/thermophysicalModels/specie/reaction/Reactions/NonEquilibriumReversibleReaction/NonEquilibriumReversibleReaction.C b/src/thermophysicalModels/specie/reaction/Reactions/NonEquilibriumReversibleReaction/NonEquilibriumReversibleReaction.C
index 255cf2b93..d3b43b8c6 100644
--- a/src/thermophysicalModels/specie/reaction/Reactions/NonEquilibriumReversibleReaction/NonEquilibriumReversibleReaction.C
+++ b/src/thermophysicalModels/specie/reaction/Reactions/NonEquilibriumReversibleReaction/NonEquilibriumReversibleReaction.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -178,6 +178,117 @@ Foam::NonEquilibriumReversibleReaction
 }
 
 
+template
+<
+    template<class> class ReactionType,
+    class ReactionThermo,
+    class ReactionRate
+>
+Foam::scalar
+Foam::NonEquilibriumReversibleReaction
+<
+    ReactionType,
+    ReactionThermo,
+    ReactionRate
+>::dkfdT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c
+) const
+{
+    return fk_.ddT(p, T, c);
+}
+
+
+template
+<
+    template<class> class ReactionType,
+    class ReactionThermo,
+    class ReactionRate
+>
+Foam::scalar
+Foam::NonEquilibriumReversibleReaction
+<
+    ReactionType,
+    ReactionThermo,
+    ReactionRate
+>::dkrdT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c,
+    const scalar dkfdT,
+    const scalar kr
+) const
+{
+    return rk_.ddT(p, T, c);
+}
+
+
+template
+<
+    template<class> class ReactionType,
+    class ReactionThermo,
+    class ReactionRate
+>
+const Foam::List<Foam::Tuple2<Foam::label, Foam::scalar>>&
+Foam::NonEquilibriumReversibleReaction
+<
+    ReactionType,
+    ReactionThermo,
+    ReactionRate
+>::beta() const
+{
+    return rk_.beta();
+}
+
+
+template
+<
+    template<class> class ReactionType,
+    class ReactionThermo,
+    class ReactionRate
+>
+void Foam::NonEquilibriumReversibleReaction
+<
+    ReactionType,
+    ReactionThermo,
+    ReactionRate
+>::dcidc
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c,
+    scalarField& dcidc
+) const
+{
+    fk_.dcidc(p, T, c, dcidc);
+}
+
+
+template
+<
+    template<class> class ReactionType,
+    class ReactionThermo,
+    class ReactionRate
+>
+Foam::scalar Foam::NonEquilibriumReversibleReaction
+<
+    ReactionType,
+    ReactionThermo,
+    ReactionRate
+>::dcidT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c
+) const
+{
+    return fk_.dcidT(p, T, c);
+}
+
+
 template
 <
     template<class> class ReactionType,
diff --git a/src/thermophysicalModels/specie/reaction/Reactions/NonEquilibriumReversibleReaction/NonEquilibriumReversibleReaction.H b/src/thermophysicalModels/specie/reaction/Reactions/NonEquilibriumReversibleReaction/NonEquilibriumReversibleReaction.H
index 8e76691d3..939c0e680 100644
--- a/src/thermophysicalModels/specie/reaction/Reactions/NonEquilibriumReversibleReaction/NonEquilibriumReversibleReaction.H
+++ b/src/thermophysicalModels/specie/reaction/Reactions/NonEquilibriumReversibleReaction/NonEquilibriumReversibleReaction.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -167,6 +167,51 @@ public:
             ) const;
 
 
+        // ReversibleReaction Jacobian functions
+
+            //- Temperature derivative of forward rate
+            virtual scalar dkfdT
+            (
+                const scalar p,
+                const scalar T,
+                const scalarField& c
+            ) const;
+
+            //- Temperature derivative of backward rate
+            virtual scalar dkrdT
+            (
+                const scalar p,
+                const scalar T,
+                const scalarField& c,
+                const scalar dkfdT,
+                const scalar kr
+            ) const;
+
+            //- Third-body efficiencies (beta = 1-alpha)
+            //  non-empty only for third-body reactions
+            //  with enhanced molecularity (alpha != 1)
+            virtual const List<Tuple2<label, scalar>>& beta() const;
+
+            //- Species concentration derivative of the pressure dependent term
+            //  By default this value is 1 as it multiplies the third-body term
+            virtual void dcidc
+            (
+                const scalar p,
+                const scalar T,
+                const scalarField& c,
+                scalarField& dcidc
+            ) const;
+
+            //- Temperature derivative of the pressure dependent term
+            //  By default this value is 0 since ddT of molecularity is approx.0
+            virtual scalar dcidT
+            (
+                const scalar p,
+                const scalar T,
+                const scalarField& c
+            ) const;
+
+
         //- Write
         virtual void write(Ostream&) const;
 };
diff --git a/src/thermophysicalModels/specie/reaction/Reactions/Reaction/Reaction.C b/src/thermophysicalModels/specie/reaction/Reactions/Reaction/Reaction.C
index dc2392eca..f17036f37 100644
--- a/src/thermophysicalModels/specie/reaction/Reactions/Reaction/Reaction.C
+++ b/src/thermophysicalModels/specie/reaction/Reactions/Reaction/Reaction.C
@@ -29,7 +29,13 @@ License
 // * * * * * * * * * * * * * * * * Static Data * * * * * * * * * * * * * * * //
 
 template<class ReactionThermo>
-Foam::label Foam::Reaction<ReactionThermo>::nUnNamedReactions = 0;
+Foam::label Foam::Reaction<ReactionThermo>::nUnNamedReactions(0);
+
+template<class ReactionThermo>
+Foam::scalar Foam::Reaction<ReactionThermo>::TlowDefault(0);
+
+template<class ReactionThermo>
+Foam::scalar Foam::Reaction<ReactionThermo>::ThighDefault(great);
 
 // * * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * //
 
@@ -160,6 +166,8 @@ Foam::Reaction<ReactionThermo>::Reaction
     ReactionThermo::thermoType(*thermoDatabase[species[0]]),
     name_("un-named-reaction-" + Foam::name(getNewReactionID())),
     species_(species),
+    Tlow_(TlowDefault),
+    Thigh_(ThighDefault),
     lhs_(lhs),
     rhs_(rhs)
 {
@@ -177,6 +185,8 @@ Foam::Reaction<ReactionThermo>::Reaction
     ReactionThermo::thermoType(r),
     name_(r.name() + "Copy"),
     species_(species),
+    Tlow_(r.Tlow()),
+    Thigh_(r.Thigh()),
     lhs_(r.lhs_),
     rhs_(r.rhs_)
 {}
@@ -197,7 +207,7 @@ Foam::Reaction<ReactionThermo>::specieCoeffs::specieCoeffs
     }
     else
     {
-        stoichCoeff = 1.0;
+        stoichCoeff = 1;
     }
 
     exponent = stoichCoeff;
@@ -330,7 +340,9 @@ Foam::Reaction<ReactionThermo>::Reaction
 :
     ReactionThermo::thermoType(*thermoDatabase[species[0]]),
     name_(dict.dictName()),
-    species_(species)
+    species_(species),
+    Tlow_(dict.lookupOrDefault<scalar>("Tlow", TlowDefault)),
+    Thigh_(dict.lookupOrDefault<scalar>("Thigh", ThighDefault))
 {
     setLRhs
     (
@@ -388,39 +400,402 @@ void Foam::Reaction<ReactionThermo>::write(Ostream& os) const
 
 
 template<class ReactionThermo>
-Foam::scalar Foam::Reaction<ReactionThermo>::kf
+void Foam::Reaction<ReactionThermo>::ddot
 (
     const scalar p,
     const scalar T,
-    const scalarField& c
+    const scalarField& c,
+    scalarField& d
 ) const
 {
-    return 0.0;
 }
 
 
 template<class ReactionThermo>
-Foam::scalar Foam::Reaction<ReactionThermo>::kr
+void Foam::Reaction<ReactionThermo>::fdot
 (
+    const scalar p,
+    const scalar T,
+    const scalarField& c,
+    scalarField& f
+) const
+{
+}
+
+
+template<class ReactionThermo>
+void Foam::Reaction<ReactionThermo>::omega
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c,
+    scalarField& dcdt
+) const
+{
+    scalar pf, cf, pr, cr;
+    label lRef, rRef;
+
+    scalar omegaI = omega
+    (
+        p, T, c, pf, cf, lRef, pr, cr, rRef
+    );
+
+    forAll(lhs_, i)
+    {
+        const label si = lhs_[i].index;
+        const scalar sl = lhs_[i].stoichCoeff;
+        dcdt[si] -= sl*omegaI;
+    }
+    forAll(rhs_, i)
+    {
+        const label si = rhs_[i].index;
+        const scalar sr = rhs_[i].stoichCoeff;
+        dcdt[si] += sr*omegaI;
+    }
+}
+
+
+template<class ReactionThermo>
+Foam::scalar Foam::Reaction<ReactionThermo>::omega
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c,
+    scalar& pf,
+    scalar& cf,
+    label& lRef,
+    scalar& pr,
+    scalar& cr,
+    label& rRef
+) const
+{
+
+    scalar clippedT = min(max(T, this->Tlow()), this->Thigh());
+
+    const scalar kf = this->kf(p, clippedT, c);
+    const scalar kr = this->kr(kf, p, clippedT, c);
+
+    pf = 1;
+    pr = 1;
+
+    const label Nl = lhs_.size();
+    const label Nr = rhs_.size();
+
+    label slRef = 0;
+    lRef = lhs_[slRef].index;
+
+    pf = kf;
+    for (label s = 1; s < Nl; s++)
+    {
+        const label si = lhs_[s].index;
+
+        if (c[si] < c[lRef])
+        {
+            const scalar exp = lhs_[slRef].exponent;
+            pf *= pow(max(c[lRef], 0), exp);
+            lRef = si;
+            slRef = s;
+        }
+        else
+        {
+            const scalar exp = lhs_[s].exponent;
+            pf *= pow(max(c[si], 0), exp);
+        }
+    }
+    cf = max(c[lRef], 0);
+
+    {
+        const scalar exp = lhs_[slRef].exponent;
+        if (exp < 1)
+        {
+            if (cf > small)
+            {
+                pf *= pow(cf, exp - 1);
+            }
+            else
+            {
+                pf = 0;
+            }
+        }
+        else
+        {
+            pf *= pow(cf, exp - 1);
+        }
+    }
+
+    label srRef = 0;
+    rRef = rhs_[srRef].index;
+
+    // Find the matrix element and element position for the rhs
+    pr = kr;
+    for (label s = 1; s < Nr; s++)
+    {
+        const label si = rhs_[s].index;
+        if (c[si] < c[rRef])
+        {
+            const scalar exp = rhs_[srRef].exponent;
+            pr *= pow(max(c[rRef], 0), exp);
+            rRef = si;
+            srRef = s;
+        }
+        else
+        {
+            const scalar exp = rhs_[s].exponent;
+            pr *= pow(max(c[si], 0), exp);
+        }
+    }
+    cr = max(c[rRef], 0);
+
+    {
+        const scalar exp = rhs_[srRef].exponent;
+        if (exp < 1)
+        {
+            if (cr > small)
+            {
+                pr *= pow(cr, exp - 1);
+            }
+            else
+            {
+                pr = 0;
+            }
+        }
+        else
+        {
+            pr *= pow(cr, exp - 1);
+        }
+    }
+
+    return pf*cf - pr*cr;
+}
+
+
+template<class ReactionThermo>
+void Foam::Reaction<ReactionThermo>::dwdc
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c,
+    scalarSquareMatrix& J,
+    scalarField& dcdt,
+    scalar& omegaI,
+    scalar& kfwd,
+    scalar& kbwd,
+    const bool reduced,
+    const List<label>& c2s
+) const
+{
+
+    scalar pf, cf, pr, cr;
+    label lRef, rRef;
+
+    omegaI = omega(p, T, c, pf, cf, lRef, pr, cr, rRef);
+
+    forAll(lhs_, i)
+    {
+        const label si = reduced ? c2s[lhs_[i].index] : lhs_[i].index;
+        const scalar sl = lhs_[i].stoichCoeff;
+        dcdt[si] -= sl*omegaI;
+    }
+    forAll(rhs_, i)
+    {
+        const label si = reduced ? c2s[rhs_[i].index] : rhs_[i].index;
+        const scalar sr = rhs_[i].stoichCoeff;
+        dcdt[si] += sr*omegaI;
+    }
+
+    kfwd = this->kf(p, T, c);
+    kbwd = this->kr(kfwd, p, T, c);
+
+    forAll(lhs_, j)
+    {
+        const label sj = reduced ? c2s[lhs_[j].index] : lhs_[j].index;
+        scalar kf = kfwd;
+        forAll(lhs_, i)
+        {
+            const label si = lhs_[i].index;
+            const scalar el = lhs_[i].exponent;
+            if (i == j)
+            {
+                if (el < 1)
+                {
+                    if (c[si] > SMALL)
+                    {
+                        kf *= el*pow(c[si] + VSMALL, el - 1);
+                    }
+                    else
+                    {
+                        kf = 0;
+                    }
+                }
+                else
+                {
+                    kf *= el*pow(c[si], el - 1);
+                }
+            }
+            else
+            {
+                kf *= pow(c[si], el);
+            }
+        }
+
+        forAll(lhs_, i)
+        {
+            const label si = reduced ? c2s[lhs_[i].index] : lhs_[i].index;
+            const scalar sl = lhs_[i].stoichCoeff;
+            J(si, sj) -= sl*kf;
+        }
+        forAll(rhs_, i)
+        {
+            const label si = reduced ? c2s[rhs_[i].index] : rhs_[i].index;
+            const scalar sr = rhs_[i].stoichCoeff;
+            J(si, sj) += sr*kf;
+        }
+    }
+
+    forAll(rhs_, j)
+    {
+        const label sj = reduced ? c2s[rhs_[j].index] : rhs_[j].index;
+        scalar kr = kbwd;
+        forAll(rhs_, i)
+        {
+            const label si = rhs_[i].index;
+            const scalar er = rhs_[i].exponent;
+            if (i == j)
+            {
+                if (er < 1)
+                {
+                    if (c[si] > SMALL)
+                    {
+                        kr *= er*pow(c[si] + VSMALL, er - 1);
+                    }
+                    else
+                    {
+                        kr = 0;
+                    }
+                }
+                else
+                {
+                    kr *= er*pow(c[si], er - 1);
+                }
+            }
+            else
+            {
+                kr *= pow(c[si], er);
+            }
+        }
+
+        forAll(lhs_, i)
+        {
+            const label si = reduced ? c2s[lhs_[i].index] : lhs_[i].index;
+            const scalar sl = lhs_[i].stoichCoeff;
+            J(si, sj) += sl*kr;
+        }
+        forAll(rhs_, i)
+        {
+            const label si = reduced ? c2s[rhs_[i].index] : rhs_[i].index;
+            const scalar sr = rhs_[i].stoichCoeff;
+            J(si, sj) -= sr*kr;
+        }
+    }
+
+    // When third-body species are involved, additional terms are added
+    // beta function returns an empty list when third-body are not involved
+    if (this->beta().size() > 0)
+    {
+        scalarField dcidc(this->beta().size());
+        this->dcidc(p, T, c, dcidc);
+        forAll(this->beta(), j)
+        {
+            label sj = this-> beta()[j].first();
+            sj = reduced ? c2s[sj] : sj;
+            if (sj != -1)
+            {
+                forAll(lhs_, i)
+                {
+                    const label si =
+                        reduced ? c2s[lhs_[i].index] : lhs_[i].index;
+                    const scalar sl = lhs_[i].stoichCoeff;
+                    J(si, sj) -= sl*dcidc[j]*omegaI;
+                }
+                forAll(rhs_, i)
+                {
+                    const label si =
+                        reduced ? c2s[rhs_[i].index] : rhs_[i].index;
+                    const scalar sr = rhs_[i].stoichCoeff;
+                    J(si, sj) += sr*dcidc[j]*omegaI;
+                }
+            }
+        }
+    }
+}
+
+
+template<class ReactionThermo>
+void Foam::Reaction<ReactionThermo>::dwdT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c,
+    const scalar omegaI,
     const scalar kfwd,
-    const scalar p,
-    const scalar T,
-    const scalarField& c
+    const scalar kbwd,
+    scalarSquareMatrix& J,
+    const bool reduced,
+    const List<label>& c2s,
+    const label indexT
 ) const
 {
-    return 0.0;
-}
+    scalar kf = kfwd;
+    scalar kr = kbwd;
 
+    scalar dkfdT = this->dkfdT(p, T, c);
+    scalar dkrdT = this->dkrdT(p, T, c, dkfdT, kr);
 
-template<class ReactionThermo>
-Foam::scalar Foam::Reaction<ReactionThermo>::kr
-(
-    const scalar p,
-    const scalar T,
-    const scalarField& c
-) const
-{
-    return 0.0;
+    scalar sumExp = 0.0;
+    forAll(lhs_, i)
+    {
+        const label si = lhs_[i].index;
+        const scalar el = lhs_[i].exponent;
+        const scalar cExp = pow(c[si], el);
+        dkfdT *= cExp;
+        kf *= cExp;
+        sumExp += el;
+    }
+    kf *= -sumExp/T;
+
+    sumExp = 0.0;
+    forAll(rhs_, i)
+    {
+        const label si = rhs_[i].index;
+        const scalar er = rhs_[i].exponent;
+        const scalar cExp = pow(c[si], er);
+        dkrdT *= cExp;
+        kr *= cExp;
+        sumExp += er;
+    }
+    kr *= -sumExp/T;
+
+    // dqidT includes the third-body (or pressure dependent) effect
+    scalar dqidT = dkfdT - dkrdT + kf - kr;
+
+    // For reactions including third-body efficiencies or pressure dependent
+    // reaction, an additional term is needed
+    scalar dcidT = this->dcidT(p, T, c);
+    dcidT *= omegaI;
+
+    // J(i, indexT) = sum_reactions nu_i dqdT
+    forAll(lhs_, i)
+    {
+        const label si = reduced ? c2s[lhs_[i].index] : lhs_[i].index;
+        const scalar sl = lhs_[i].stoichCoeff;
+        J(si, indexT) -= sl*(dqidT + dcidT);
+    }
+    forAll(rhs_, i)
+    {
+        const label si = reduced ? c2s[rhs_[i].index] : rhs_[i].index;
+        const scalar sr = rhs_[i].stoichCoeff;
+        J(si, indexT) += sr*(dqidT + dcidT);
+    }
 }
 
 
diff --git a/src/thermophysicalModels/specie/reaction/Reactions/Reaction/Reaction.H b/src/thermophysicalModels/specie/reaction/Reactions/Reaction/Reaction.H
index 13f6d9a53..dce5e7688 100644
--- a/src/thermophysicalModels/specie/reaction/Reactions/Reaction/Reaction.H
+++ b/src/thermophysicalModels/specie/reaction/Reactions/Reaction/Reaction.H
@@ -40,6 +40,8 @@ SourceFiles
 #include "speciesTable.H"
 #include "HashPtrTable.H"
 #include "scalarField.H"
+#include "simpleMatrix.H"
+#include "Tuple2.H"
 #include "typeInfo.H"
 #include "runTimeSelectionTables.H"
 
@@ -84,6 +86,8 @@ public:
         //- Number of un-named reactions
         static label nUnNamedReactions;
 
+        //- Default temperature limits of applicability of reaction rates
+        static scalar TlowDefault, ThighDefault;
 
     // Public data types
 
@@ -134,6 +138,9 @@ private:
         //- List of specie names present in reaction system
         const speciesTable& species_;
 
+        //- Temperature limits of applicability of reaction rates
+        scalar Tlow_, Thigh_;
+
         //- Specie info for the left-hand-side of the reaction
         List<specieCoeffs> lhs_;
 
@@ -201,25 +208,13 @@ public:
         );
 
         //- Construct and return a clone
-        virtual autoPtr<Reaction<ReactionThermo>> clone() const //
-        {
-            return autoPtr<Reaction<ReactionThermo>>
-            (
-                new Reaction<ReactionThermo>(*this)
-            );
-        }
+        virtual autoPtr<Reaction<ReactionThermo>> clone() const = 0;
 
         //- Construct and return a clone with new speciesTable
         virtual autoPtr<Reaction<ReactionThermo>> clone
         (
             const speciesTable& species
-        ) const
-        {
-            return autoPtr<Reaction<ReactionThermo>>
-            (
-                new Reaction<ReactionThermo>(*this, species)
-            );
-        }
+        ) const = 0;
 
 
     // Selectors
@@ -242,14 +237,22 @@ public:
 
         // Access
 
+            inline word& name();
             inline const word& name() const;
 
-            // - Access to basic components of the reaction
+            inline scalar Tlow() const;
+            inline scalar Thigh() const;
+
+            //- Access to the components of the left hand side
             inline const List<specieCoeffs>& lhs() const;
+
+            //- Access to the components of the right hand side
             inline const List<specieCoeffs>& rhs() const;
 
-            // - Access to gas components of the reaction
+            //- Access to the gas components of the left hand side
             virtual const List<specieCoeffs>& grhs() const;
+
+            //- Access to the gas components of the right hand side
             virtual const List<specieCoeffs>& glhs() const;
 
             //- Access to specie list
@@ -268,6 +271,49 @@ public:
             );
 
 
+        // Reaction rate coefficients
+
+            //- Forward reaction rate
+            void ddot
+            (
+                const scalar p,
+                const scalar T,
+                const scalarField& c,
+                scalarField& d
+            ) const;
+
+            //- Backward reaction rate
+            void fdot
+            (
+                const scalar p,
+                const scalar T,
+                const scalarField& c,
+                scalarField& f
+            ) const;
+
+            //- Net reaction rate for individual species
+            void omega
+            (
+                const scalar p,
+                const scalar T,
+                const scalarField& c,
+                scalarField& dcdt
+            ) const;
+
+            //- Net reaction rate
+            scalar omega
+            (
+                const scalar p,
+                const scalar T,
+                const scalarField& c,
+                scalar& pf,
+                scalar& cf,
+                label& lRef,
+                scalar& pr,
+                scalar& cr,
+                label& rRef
+            ) const;
+
         // Reaction rate coefficients
 
             //- Forward rate constant
@@ -276,7 +322,7 @@ public:
                 const scalar p,
                 const scalar T,
                 const scalarField& c
-            ) const;
+            ) const = 0;
 
             //- Reverse rate constant from the given forward rate constant
             virtual scalar kr
@@ -285,18 +331,91 @@ public:
                 const scalar p,
                 const scalar T,
                 const scalarField& c
-            ) const;
+            ) const = 0;
 
-            //- Reverse rate constant.
-            //  Note this evaluates the forward rate constant and divides by the
-            //  equilibrium constant
+            //- Reverse rate constant
             virtual scalar kr
             (
                 const scalar p,
                 const scalar T,
                 const scalarField& c
+            ) const = 0;
+
+
+        // Jacobian coefficients
+
+            //- Derivative of the net reaction rate for each species involved
+            //  w.r.t. the species concentration
+            void dwdc
+            (
+                const scalar p,
+                const scalar T,
+                const scalarField& c,
+                scalarSquareMatrix& J,
+                scalarField& dcdt,
+                scalar& omegaI,
+                scalar& kfwd,
+                scalar& kbwd,
+                const bool reduced,
+                const List<label>& c2s
             ) const;
 
+            //- Derivative of the net reaction rate for each species involved
+            //  w.r.t. the temperature
+            void dwdT
+            (
+                const scalar p,
+                const scalar T,
+                const scalarField& c,
+                const scalar omegaI,
+                const scalar kfwd,
+                const scalar kbwd,
+                scalarSquareMatrix& J,
+                const bool reduced,
+                const List<label>& c2s,
+                const label indexT
+            ) const;
+
+            //- Temperature derivative of forward rate
+            virtual scalar dkfdT
+            (
+                const scalar p,
+                const scalar T,
+                const scalarField& c
+            ) const = 0;
+
+            //- Temperature derivative of reverse rate
+            virtual scalar dkrdT
+            (
+                const scalar p,
+                const scalar T,
+                const scalarField& c,
+                const scalar dkfdT,
+                const scalar kr
+            ) const = 0;
+
+            //- Third-body efficiencies (beta = 1-alpha)
+            //  non-empty only for third-body reactions
+            //  with enhanced molecularity (alpha != 1)
+            virtual const List<Tuple2<label, scalar>>& beta() const = 0;
+
+            //- Species concentration derivative of the pressure dependent term
+            virtual void dcidc
+            (
+                const scalar p,
+                const scalar T,
+                const scalarField& c,
+                scalarField& dcidc
+            ) const = 0;
+
+            //- Temperature derivative of the pressure dependent term
+            virtual scalar dcidT
+            (
+                const scalar p,
+                const scalar T,
+                const scalarField& c
+            ) const = 0;
+
 
         //- Write
         virtual void write(Ostream&) const;
diff --git a/src/thermophysicalModels/specie/reaction/Reactions/Reaction/ReactionI.H b/src/thermophysicalModels/specie/reaction/Reactions/Reaction/ReactionI.H
index eba046f64..b52565c63 100644
--- a/src/thermophysicalModels/specie/reaction/Reactions/Reaction/ReactionI.H
+++ b/src/thermophysicalModels/specie/reaction/Reactions/Reaction/ReactionI.H
@@ -32,6 +32,13 @@ namespace Foam
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+template<class ReactionThermo>
+inline word& Reaction<ReactionThermo>::name()
+{
+    return name_;
+}
+
+
 template<class ReactionThermo>
 inline const word& Reaction<ReactionThermo>::name() const
 {
@@ -39,6 +46,20 @@ inline const word& Reaction<ReactionThermo>::name() const
 }
 
 
+template<class ReactionThermo>
+inline scalar Reaction<ReactionThermo>::Tlow() const
+{
+    return Tlow_;
+}
+
+
+template<class ReactionThermo>
+inline scalar Reaction<ReactionThermo>::Thigh() const
+{
+    return Thigh_;
+}
+
+
 template<class ReactionThermo>
 inline const List<typename Reaction<ReactionThermo>::specieCoeffs>&
 Reaction<ReactionThermo>::lhs() const
diff --git a/src/thermophysicalModels/specie/reaction/Reactions/ReactionList/ReactionList.C b/src/thermophysicalModels/specie/reaction/Reactions/ReactionList/ReactionList.C
index 6d3f2539f..cfae96819 100644
--- a/src/thermophysicalModels/specie/reaction/Reactions/ReactionList/ReactionList.C
+++ b/src/thermophysicalModels/specie/reaction/Reactions/ReactionList/ReactionList.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -84,6 +84,13 @@ bool Foam::ReactionList<ThermoType>::readReactionDict()
 {
     const dictionary& reactions(dict_.subDict("reactions"));
 
+    // Set general temperature limits from the dictionary
+    Reaction<ThermoType>::TlowDefault =
+        dict_.lookupOrDefault<scalar>("Tlow", 0);
+
+    Reaction<ThermoType>::ThighDefault =
+        dict_.lookupOrDefault<scalar>("Thigh", great);
+
     forAllConstIter(dictionary, reactions, iter)
     {
         const word reactionName = iter().keyword();
diff --git a/src/thermophysicalModels/specie/reaction/Reactions/ReactionProxy/ReactionProxy.C b/src/thermophysicalModels/specie/reaction/Reactions/ReactionProxy/ReactionProxy.C
new file mode 100644
index 000000000..992521fd5
--- /dev/null
+++ b/src/thermophysicalModels/specie/reaction/Reactions/ReactionProxy/ReactionProxy.C
@@ -0,0 +1,207 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2018 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "ReactionProxy.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class ReactionThermo>
+Foam::ReactionProxy<ReactionThermo>::ReactionProxy
+(
+    const speciesTable& species,
+    const List<typename Reaction<ReactionThermo>::specieCoeffs>& lhs,
+    const List<typename Reaction<ReactionThermo>::specieCoeffs>& rhs,
+    const HashPtrTable<ReactionThermo>& thermoDatabase
+)
+:
+    Reaction<ReactionThermo>
+    (
+        species,
+        lhs,
+        rhs,
+        thermoDatabase
+    )
+{}
+
+
+template<class ReactionThermo>
+Foam::ReactionProxy<ReactionThermo>::ReactionProxy
+(
+    const Reaction<ReactionThermo>& r,
+    const speciesTable& species
+)
+:
+    Reaction<ReactionThermo>
+    (
+        r,
+        species
+    )
+{}
+
+
+template<class ReactionThermo>
+Foam::ReactionProxy<ReactionThermo>::ReactionProxy
+(
+    const speciesTable& species,
+    const HashPtrTable<ReactionThermo>& thermoDatabase,
+    const dictionary& dict
+)
+:
+    Reaction<ReactionThermo>
+    (
+        species,
+        thermoDatabase,
+        dict
+    )
+{}
+
+
+template<class ReactionThermo>
+Foam::autoPtr<Foam::Reaction<ReactionThermo>>
+Foam::ReactionProxy<ReactionThermo>::clone() const
+{
+    NotImplemented;
+    return autoPtr<Foam::Reaction<ReactionThermo>>();
+}
+
+
+template<class ReactionThermo>
+Foam::autoPtr<Foam::Reaction<ReactionThermo>>
+Foam::ReactionProxy<ReactionThermo>::clone
+(
+    const speciesTable& species
+) const
+{
+    NotImplemented;
+    return autoPtr<Reaction<ReactionThermo>>();
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class ReactionThermo>
+Foam::scalar Foam::ReactionProxy<ReactionThermo>::kf
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c
+) const
+{
+    NotImplemented;
+    return 0;
+}
+
+
+template<class ReactionThermo>
+Foam::scalar Foam::ReactionProxy<ReactionThermo>::kr
+(
+    const scalar kfwd,
+    const scalar p,
+    const scalar T,
+    const scalarField& c
+) const
+{
+    NotImplemented;
+    return 0;
+}
+
+
+template<class ReactionThermo>
+Foam::scalar Foam::ReactionProxy<ReactionThermo>::kr
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c
+) const
+{
+    NotImplemented;
+    return 0;
+}
+
+
+template<class ReactionThermo>
+Foam::scalar Foam::ReactionProxy<ReactionThermo>::dkfdT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c
+) const
+{
+    NotImplemented;
+    return 0;
+}
+
+
+template<class ReactionThermo>
+Foam::scalar Foam::ReactionProxy<ReactionThermo>::dkrdT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c,
+    const scalar dkfdT,
+    const scalar kr
+) const
+{
+    NotImplemented;
+    return 0;
+}
+
+
+template<class ReactionThermo>
+const Foam::List<Foam::Tuple2<Foam::label, Foam::scalar>>&
+Foam::ReactionProxy<ReactionThermo>::beta() const
+{
+    NotImplemented;
+    return NullObjectRef<List<Tuple2<label, scalar>>>();
+}
+
+
+template<class ReactionThermo>
+void Foam::ReactionProxy<ReactionThermo>::dcidc
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c,
+    scalarField& dcidc
+) const
+{
+    NotImplemented;
+}
+
+
+template<class ReactionThermo>
+Foam::scalar Foam::ReactionProxy<ReactionThermo>::dcidT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c
+) const
+{
+    NotImplemented;
+    return 0;
+}
+
+
+// ************************************************************************* //
diff --git a/src/thermophysicalModels/specie/reaction/Reactions/ReactionProxy/ReactionProxy.H b/src/thermophysicalModels/specie/reaction/Reactions/ReactionProxy/ReactionProxy.H
new file mode 100644
index 000000000..2349dd27c
--- /dev/null
+++ b/src/thermophysicalModels/specie/reaction/Reactions/ReactionProxy/ReactionProxy.H
@@ -0,0 +1,189 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2018 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::ReactionProxy
+
+Description
+    Proxy version of Reaction which provides dummy implementations of the
+    abstract virtual functions.
+
+    Used for read-construction and format conversion.
+
+SourceFiles
+    ReactionProxy.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef ReactionProxy_H
+#define ReactionProxy_H
+
+#include "Reaction.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                        Class ReactionProxy Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class ReactionThermo>
+class ReactionProxy
+:
+    public Reaction<ReactionThermo>
+{
+
+public:
+
+    // Constructors
+
+        //- Construct from components
+        ReactionProxy
+        (
+            const speciesTable& species,
+            const List<typename Reaction<ReactionThermo>::specieCoeffs>& lhs,
+            const List<typename Reaction<ReactionThermo>::specieCoeffs>& rhs,
+            const HashPtrTable<ReactionThermo>& thermoDatabase
+        );
+
+        //- Construct as copy given new speciesTable
+        ReactionProxy
+        (
+            const Reaction<ReactionThermo>&,
+            const speciesTable& species
+        );
+
+        //- Construct from dictionary
+        ReactionProxy
+        (
+            const speciesTable& species,
+            const HashPtrTable<ReactionThermo>& thermoDatabase,
+            const dictionary& dict
+        );
+
+        //- Construct and return a clone
+        virtual autoPtr<Reaction<ReactionThermo>> clone() const;
+
+        //- Construct and return a clone with new speciesTable
+        virtual autoPtr<Reaction<ReactionThermo>> clone
+        (
+            const speciesTable& species
+        ) const;
+
+
+    //- Destructor
+    virtual ~ReactionProxy()
+    {}
+
+
+    // Member Functions
+
+        // Reaction rate coefficients
+
+            //- Forward rate constant
+            virtual scalar kf
+            (
+                const scalar p,
+                const scalar T,
+                const scalarField& c
+            ) const;
+
+            //- Reverse rate constant from the given forward rate constant
+            virtual scalar kr
+            (
+                const scalar kfwd,
+                const scalar p,
+                const scalar T,
+                const scalarField& c
+            ) const;
+
+            //- Reverse rate constant
+            virtual scalar kr
+            (
+                const scalar p,
+                const scalar T,
+                const scalarField& c
+            ) const;
+
+
+        // Jacobian coefficients
+
+            //- Temperature derivative of forward rate
+            virtual scalar dkfdT
+            (
+                const scalar p,
+                const scalar T,
+                const scalarField& c
+            ) const;
+
+            //- Temperature derivative of reverse rate
+            virtual scalar dkrdT
+            (
+                const scalar p,
+                const scalar T,
+                const scalarField& c,
+                const scalar dkfdT,
+                const scalar kr
+            ) const;
+
+            //- Third-body efficiencies (beta = 1-alpha)
+            //  non-empty only for third-body reactions
+            //  with enhanced molecularity (alpha != 1)
+            virtual const List<Tuple2<label, scalar>>& beta() const;
+
+            //- Species concentration derivative of the pressure dependent term
+            virtual void dcidc
+            (
+                const scalar p,
+                const scalar T,
+                const scalarField& c,
+                scalarField& dcidc
+            ) const;
+
+            //- Temperature derivative of the pressure dependent term
+            virtual scalar dcidT
+            (
+                const scalar p,
+                const scalar T,
+                const scalarField& c
+            ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "ReactionProxy.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/thermophysicalModels/specie/reaction/Reactions/ReversibleReaction/ReversibleReaction.C b/src/thermophysicalModels/specie/reaction/Reactions/ReversibleReaction/ReversibleReaction.C
index 2ccc2fd88..58964a9ba 100644
--- a/src/thermophysicalModels/specie/reaction/Reactions/ReversibleReaction/ReversibleReaction.C
+++ b/src/thermophysicalModels/specie/reaction/Reactions/ReversibleReaction/ReversibleReaction.C
@@ -151,6 +151,117 @@ Foam::scalar Foam::ReversibleReaction
 }
 
 
+template
+<
+    template<class> class ReactionType,
+    class ReactionThermo,
+    class ReactionRate
+>
+Foam::scalar Foam::ReversibleReaction
+<
+    ReactionType,
+    ReactionThermo,
+    ReactionRate
+>::dkfdT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c
+) const
+{
+    return k_.ddT(p, T, c);
+}
+
+
+template
+<
+    template<class> class ReactionType,
+    class ReactionThermo,
+    class ReactionRate
+>
+Foam::scalar Foam::ReversibleReaction
+<
+    ReactionType,
+    ReactionThermo,
+    ReactionRate
+>::dkrdT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c,
+    const scalar dkfdT,
+    const scalar kr
+) const
+{
+    scalar Kc = max(this->Kc(p, T), rootSmall);
+
+    return dkfdT/Kc - kr*this->dKcdTbyKc(p, T);
+}
+
+
+template
+<
+    template<class> class ReactionType,
+    class ReactionThermo,
+    class ReactionRate
+>
+const Foam::List<Foam::Tuple2<Foam::label, Foam::scalar>>&
+Foam::ReversibleReaction
+<
+    ReactionType,
+    ReactionThermo,
+    ReactionRate
+>::beta() const
+{
+    return k_.beta();
+}
+
+
+template
+<
+    template<class> class ReactionType,
+    class ReactionThermo,
+    class ReactionRate
+>
+void Foam::ReversibleReaction
+<
+    ReactionType,
+    ReactionThermo,
+    ReactionRate
+>::dcidc
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c,
+    scalarField& dcidc
+) const
+{
+    k_.dcidc(p, T, c, dcidc);
+}
+
+
+template
+<
+    template<class> class ReactionType,
+    class ReactionThermo,
+    class ReactionRate
+>
+Foam::scalar Foam::ReversibleReaction
+<
+    ReactionType,
+    ReactionThermo,
+    ReactionRate
+>::dcidT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c
+) const
+{
+    return k_.dcidT(p, T, c);
+}
+
+
 template
 <
     template<class> class ReactionType,
diff --git a/src/thermophysicalModels/specie/reaction/Reactions/ReversibleReaction/ReversibleReaction.H b/src/thermophysicalModels/specie/reaction/Reactions/ReversibleReaction/ReversibleReaction.H
index 798b75759..9cb1c4d64 100644
--- a/src/thermophysicalModels/specie/reaction/Reactions/ReversibleReaction/ReversibleReaction.H
+++ b/src/thermophysicalModels/specie/reaction/Reactions/ReversibleReaction/ReversibleReaction.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -180,6 +180,51 @@ public:
             ) const;
 
 
+        // ReversibleReaction Jacobian functions
+
+            //- Temperature derivative of forward rate
+            virtual scalar dkfdT
+            (
+                const scalar p,
+                const scalar T,
+                const scalarField& c
+            ) const;
+
+            //- Temperature derivative of backward rate
+            virtual scalar dkrdT
+            (
+                const scalar p,
+                const scalar T,
+                const scalarField& c,
+                const scalar dkfdT,
+                const scalar kr
+            ) const;
+
+            //- Third-body efficiencies (beta = 1-alpha)
+            //  non-empty only for third-body reactions
+            //  with enhanced molecularity (alpha != 1)
+            virtual const List<Tuple2<label, scalar>>& beta() const;
+
+            //- Species concentration derivative of the pressure dependent term
+            //  By default this value is 1 as it multiplies the third-body term
+            virtual void dcidc
+            (
+                const scalar p,
+                const scalar T,
+                const scalarField& c,
+                scalarField& dcidc
+            ) const;
+
+            //- Temperature derivative of the pressure dependent term
+            //  By default this value is 0 since ddT of molecularity is approx.0
+            virtual scalar dcidT
+            (
+                const scalar p,
+                const scalar T,
+                const scalarField& c
+            ) const;
+
+
         //- Write
         virtual void write(Ostream&) const;
 };
diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/ArrheniusReactionRate/ArrheniusReactionRate.H b/src/thermophysicalModels/specie/reaction/reactionRate/ArrheniusReactionRate/ArrheniusReactionRate.H
index fbe89154b..604161773 100644
--- a/src/thermophysicalModels/specie/reaction/reactionRate/ArrheniusReactionRate/ArrheniusReactionRate.H
+++ b/src/thermophysicalModels/specie/reaction/reactionRate/ArrheniusReactionRate/ArrheniusReactionRate.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -100,6 +100,37 @@ public:
             const scalarField& c
         ) const;
 
+        inline scalar ddT
+        (
+            const scalar p,
+            const scalar T,
+            const scalarField& c
+        ) const;
+
+        //- Third-body efficiencies (beta = 1-alpha)
+        //  non-empty only for third-body reactions
+        //  with enhanced molecularity (alpha != 1)
+        inline const List<Tuple2<label, scalar>>& beta() const;
+
+        //- Species concentration derivative of the pressure dependent term
+        //  By default this value is 1 as it multiplies the third-body term
+        inline void dcidc
+        (
+            const scalar p,
+            const scalar T,
+            const scalarField& c,
+            scalarField& dcidc
+        ) const;
+
+        //- Temperature derivative of the pressure dependent term
+        //  By default this value is 0 since ddT of molecularity is approx.0
+        inline scalar dcidT
+        (
+            const scalar p,
+            const scalar T,
+            const scalarField& c
+        ) const;
+
         //- Write to stream
         inline void write(Ostream& os) const;
 
diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/ArrheniusReactionRate/ArrheniusReactionRateI.H b/src/thermophysicalModels/specie/reaction/reactionRate/ArrheniusReactionRate/ArrheniusReactionRateI.H
index a856eb911..61bb2e1b6 100644
--- a/src/thermophysicalModels/specie/reaction/reactionRate/ArrheniusReactionRate/ArrheniusReactionRateI.H
+++ b/src/thermophysicalModels/specie/reaction/reactionRate/ArrheniusReactionRate/ArrheniusReactionRateI.H
@@ -75,6 +75,57 @@ inline Foam::scalar Foam::ArrheniusReactionRate::operator()
 }
 
 
+inline Foam::scalar Foam::ArrheniusReactionRate::ddT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField&
+) const
+{
+    scalar ak = A_;
+
+    if (mag(beta_) > vSmall)
+    {
+        ak *= pow(T, beta_);
+    }
+
+    if (mag(Ta_) > vSmall)
+    {
+        ak *= exp(-Ta_/T);
+    }
+
+    return ak*(beta_+Ta_/T)/T;
+}
+
+
+inline const Foam::List<Foam::Tuple2<Foam::label, Foam::scalar>>&
+Foam::ArrheniusReactionRate::beta() const
+{
+    return NullObjectRef<List<Tuple2<label, scalar>>>();
+}
+
+
+inline void Foam::ArrheniusReactionRate::dcidc
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c,
+    scalarField& dcidc
+) const
+{}
+
+
+inline Foam::scalar Foam::ArrheniusReactionRate::dcidT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c
+) const
+{
+    return 0;
+}
+
+
 inline void Foam::ArrheniusReactionRate::write(Ostream& os) const
 {
     os.writeKeyword("A") << A_ << token::END_STATEMENT << nl;
diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/ChemicallyActivatedReactionRate/ChemicallyActivatedReactionRate.H b/src/thermophysicalModels/specie/reaction/reactionRate/ChemicallyActivatedReactionRate/ChemicallyActivatedReactionRate.H
index 87eab0d40..e70a3977b 100644
--- a/src/thermophysicalModels/specie/reaction/reactionRate/ChemicallyActivatedReactionRate/ChemicallyActivatedReactionRate.H
+++ b/src/thermophysicalModels/specie/reaction/reactionRate/ChemicallyActivatedReactionRate/ChemicallyActivatedReactionRate.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -69,6 +69,7 @@ class ChemicallyActivatedReactionRate
         ReactionRate kInf_;
         ChemicallyActivationFunction F_;
         thirdBodyEfficiencies thirdBodyEfficiencies_;
+        List<Tuple2<label, scalar>> beta_;
 
 
 public:
@@ -109,6 +110,40 @@ public:
             const scalarField& c
         ) const;
 
+        inline scalar ddT
+        (
+            const scalar p,
+            const scalar T,
+            const scalarField& c
+        ) const;
+
+        //- Third-body efficiencies (beta = 1-alpha)
+        //  non-empty only for third-body reactions
+        //  with enhanced molecularity (alpha != 1)
+        inline const List<Tuple2<label, scalar>>& beta() const
+        {
+            return beta_;
+        }
+
+        //- Species concentration derivative of the pressure dependent term
+        //  By default this value is 1 as it multiplies the third-body term
+        inline void dcidc
+        (
+            const scalar p,
+            const scalar T,
+            const scalarField& c,
+            scalarField& dcidc
+        ) const;
+
+        //- Temperature derivative of the pressure dependent term
+        //  By default this value is 0 since ddT of molecularity is approx.0
+        inline scalar dcidT
+        (
+            const scalar p,
+            const scalar T,
+            const scalarField& c
+        ) const;
+
         //- Write to stream
         inline void write(Ostream& os) const;
 
diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/ChemicallyActivatedReactionRate/ChemicallyActivatedReactionRateI.H b/src/thermophysicalModels/specie/reaction/reactionRate/ChemicallyActivatedReactionRate/ChemicallyActivatedReactionRateI.H
index 3e54f9a6a..d85a8c5bb 100644
--- a/src/thermophysicalModels/specie/reaction/reactionRate/ChemicallyActivatedReactionRate/ChemicallyActivatedReactionRateI.H
+++ b/src/thermophysicalModels/specie/reaction/reactionRate/ChemicallyActivatedReactionRate/ChemicallyActivatedReactionRateI.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -42,7 +42,12 @@ inline Foam::ChemicallyActivatedReactionRate
     kInf_(kInf),
     F_(F),
     thirdBodyEfficiencies_(tbes)
-{}
+{
+    forAll(tbes, i)
+    {
+        beta_.append(Tuple2<label, scalar>(i, tbes[i]));
+    }
+}
 
 
 template<class ReactionRate, class ChemicallyActivationFunction>
@@ -60,7 +65,19 @@ inline Foam::ChemicallyActivatedReactionRate
     kInf_(species, dict),
     F_(dict),
     thirdBodyEfficiencies_(species, dict)
-{}
+{
+    forAll(thirdBodyEfficiencies_, i)
+    {
+        beta_.append
+        (
+            Tuple2<label, scalar>
+            (
+                i,
+                thirdBodyEfficiencies_[i]
+            )
+        );
+    }
+}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
@@ -77,15 +94,108 @@ inline Foam::scalar Foam::ChemicallyActivatedReactionRate
     const scalarField& c
 ) const
 {
-    scalar k0 = k0_(p, T, c);
-    scalar kInf = kInf_(p, T, c);
-
-    scalar Pr = k0*thirdBodyEfficiencies_.M(c)/kInf;
+    const scalar k0 = k0_(p, T, c);
+    const scalar kInf = kInf_(p, T, c);
+    const scalar Pr = k0*thirdBodyEfficiencies_.M(c)/kInf;
 
     return k0*(1/(1 + Pr))*F_(T, Pr);
 }
 
 
+template<class ReactionRate, class ChemicallyActivationFunction>
+inline Foam::scalar Foam::ChemicallyActivatedReactionRate
+<
+    ReactionRate,
+    ChemicallyActivationFunction
+>::ddT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c
+) const
+{
+    const scalar k0 = k0_(p, T, c);
+    const scalar kInf = kInf_(p, T, c);
+    const scalar Pr = k0*thirdBodyEfficiencies_.M(c)/kInf;
+
+    return (1/(1 + Pr))*F_(T, Pr)*k0_.ddT(p, T, c);
+}
+
+
+template<class ReactionRate, class ChemicallyActivationFunction>
+inline void Foam::ChemicallyActivatedReactionRate
+<
+    ReactionRate,
+    ChemicallyActivationFunction
+>::dcidc
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c,
+    scalarField& dcidc
+) const
+{
+    const scalar M = thirdBodyEfficiencies_.M(c);
+
+    if (M > small)
+    {
+        const scalar k0 = k0_(p, T, c);
+        const scalar kInf = kInf_(p, T, c);
+
+        const scalar Pr = k0*M/kInf;
+        const scalar F = F_(T, Pr);
+
+        forAll(beta_, i)
+        {
+            const scalar dPrdci = -beta_[i].second()*k0/kInf;
+            const scalar dFdci = F_.ddc(Pr, F, dPrdci, T);
+            dcidc[i] = (-dPrdci/(1 + Pr) + dFdci/F);
+        }
+    }
+    else
+    {
+        forAll(beta_, i)
+        {
+            dcidc[i] = 0;
+        }
+    }
+}
+
+
+template<class ReactionRate, class ChemicallyActivationFunction>
+inline Foam::scalar Foam::ChemicallyActivatedReactionRate
+<
+    ReactionRate,
+    ChemicallyActivationFunction
+>::dcidT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c
+) const
+{
+    const scalar M = thirdBodyEfficiencies_.M(c);
+
+    if (M > small)
+    {
+        const scalar k0 = k0_(p, T, c);
+        const scalar kInf = kInf_(p, T, c);
+
+        const scalar Pr = k0*thirdBodyEfficiencies_.M(c)/kInf;
+        const scalar F = F_(T, Pr);
+        const scalar dPrdT =
+            Pr*(k0_.ddT(p, T, c)/k0 - kInf_.ddT(p, T, c)/kInf - 1/T);
+        const scalar dFdT = F_.ddT(Pr, F, dPrdT, T);
+
+        return (-dPrdT/(1 + Pr) + dFdT/F);
+    }
+    else
+    {
+        return 0;
+    }
+}
+
+
 template<class ReactionRate, class ChemicallyActivationFunction>
 inline void Foam::ChemicallyActivatedReactionRate
 <
diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/FallOffReactionRate/FallOffReactionRate.H b/src/thermophysicalModels/specie/reaction/reactionRate/FallOffReactionRate/FallOffReactionRate.H
index e89aa8949..eb2498f60 100644
--- a/src/thermophysicalModels/specie/reaction/reactionRate/FallOffReactionRate/FallOffReactionRate.H
+++ b/src/thermophysicalModels/specie/reaction/reactionRate/FallOffReactionRate/FallOffReactionRate.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -68,6 +68,7 @@ class FallOffReactionRate
         ReactionRate kInf_;
         FallOffFunction F_;
         thirdBodyEfficiencies thirdBodyEfficiencies_;
+        List<Tuple2<label, scalar>> beta_;
 
 
 public:
@@ -106,6 +107,38 @@ public:
             const scalarField& c
         ) const;
 
+        inline scalar ddT
+        (
+            const scalar p,
+            const scalar T,
+            const scalarField& c
+        ) const;
+
+        //- Third-body efficiencies (beta = 1-alpha)
+        //  non-empty only for third-body reactions
+        //  with enhanced molecularity (alpha != 1)
+        inline const List<Tuple2<label, scalar>>& beta() const
+        {
+            return beta_;
+        }
+
+        //- Species concentration derivative of the pressure dependent term
+        inline void dcidc
+        (
+            const scalar p,
+            const scalar T,
+            const scalarField& c,
+            scalarField& dcidc
+        ) const;
+
+        //- Temperature derivative of the pressure dependent term
+        inline scalar dcidT
+        (
+            const scalar p,
+            const scalar T,
+            const scalarField& c
+        ) const;
+
         //- Write to stream
         inline void write(Ostream& os) const;
 
diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/FallOffReactionRate/FallOffReactionRateI.H b/src/thermophysicalModels/specie/reaction/reactionRate/FallOffReactionRate/FallOffReactionRateI.H
index fff348d48..228776e61 100644
--- a/src/thermophysicalModels/specie/reaction/reactionRate/FallOffReactionRate/FallOffReactionRateI.H
+++ b/src/thermophysicalModels/specie/reaction/reactionRate/FallOffReactionRate/FallOffReactionRateI.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -39,7 +39,12 @@ FallOffReactionRate
     kInf_(kInf),
     F_(F),
     thirdBodyEfficiencies_(tbes)
-{}
+{
+    forAll(tbes, i)
+    {
+        beta_.append(Tuple2<label, scalar>(i, tbes[i]));
+    }
+}
 
 
 template<class ReactionRate, class FallOffFunction>
@@ -54,7 +59,19 @@ FallOffReactionRate
     kInf_(species, dict.subDict("kInf")),
     F_(dict.subDict("F")),
     thirdBodyEfficiencies_(species, dict.subDict("thirdBodyEfficiencies"))
-{}
+{
+    forAll(thirdBodyEfficiencies_, i)
+    {
+        beta_.append
+        (
+            Tuple2<label, scalar>
+            (
+                i,
+                thirdBodyEfficiencies_[i]
+            )
+        );
+    }
+}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
@@ -68,15 +85,97 @@ Foam::FallOffReactionRate<ReactionRate, FallOffFunction>::operator()
     const scalarField& c
 ) const
 {
-    scalar k0 = k0_(p, T, c);
-    scalar kInf = kInf_(p, T, c);
-
-    scalar Pr = k0*thirdBodyEfficiencies_.M(c)/kInf;
+    const scalar k0 = k0_(p, T, c);
+    const scalar kInf = kInf_(p, T, c);
+    const scalar Pr = k0*thirdBodyEfficiencies_.M(c)/kInf;
 
     return kInf*(Pr/(1 + Pr))*F_(T, Pr);
 }
 
 
+template<class ReactionRate, class FallOffFunction>
+inline Foam::scalar
+Foam::FallOffReactionRate<ReactionRate, FallOffFunction>::ddT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c
+) const
+{
+    const scalar k0 = k0_(p, T, c);
+    const scalar kInf = kInf_(p, T, c);
+    const scalar Pr = k0*thirdBodyEfficiencies_.M(c)/kInf;
+
+    return (Pr/(1 + Pr))*F_(T, Pr)*kInf_.ddT(p, T, c);
+}
+
+
+template<class ReactionRate, class FallOffFunction>
+inline void Foam::FallOffReactionRate<ReactionRate, FallOffFunction>::dcidc
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c,
+    scalarField& dcidc
+) const
+{
+    const scalar M = thirdBodyEfficiencies_.M(c);
+
+    if (M > small)
+    {
+        const scalar k0 = k0_(p, T, c);
+        const scalar kInf = kInf_(p, T, c);
+        const scalar Pr = k0*M/kInf;
+        const scalar F = F_(T, Pr);
+
+        forAll(beta_, i)
+        {
+            const scalar dPrdci = -beta_[i].second()*k0/kInf;
+            const scalar dFdci = F_.ddc(Pr, F, dPrdci, T);
+            dcidc[i] = (dPrdci/(Pr*(1 + Pr)) + dFdci/F);
+        }
+    }
+    else
+    {
+        forAll(beta_, i)
+        {
+            dcidc[i] = 0;
+        }
+    }
+}
+
+
+template<class ReactionRate, class FallOffFunction>
+inline Foam::scalar
+Foam::FallOffReactionRate<ReactionRate, FallOffFunction>::dcidT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c
+) const
+{
+    const scalar M = thirdBodyEfficiencies_.M(c);
+
+    if (M > small)
+    {
+        const scalar k0 = k0_(p, T, c);
+        const scalar kInf = kInf_(p, T, c);
+
+        const scalar Pr = k0*thirdBodyEfficiencies_.M(c)/kInf;
+        const scalar F = F_(T, Pr);
+        const scalar dPrdT =
+            Pr*(k0_.ddT(p, T, c)/k0 - kInf_.ddT(p, T, c)/kInf - 1/T);
+        const scalar dFdT = F_.ddT(Pr, F, dPrdT, T);
+
+        return (dPrdT/(Pr*(1 + Pr)) + dFdT/F);
+    }
+    else
+    {
+        return 0;
+    }
+}
+
+
 template<class ReactionRate, class FallOffFunction>
 inline void Foam::FallOffReactionRate<ReactionRate, FallOffFunction>::write
 (
diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/JanevReactionRate/JanevReactionRate.H b/src/thermophysicalModels/specie/reaction/reactionRate/JanevReactionRate/JanevReactionRate.H
index ab37e054f..a1b813ab4 100644
--- a/src/thermophysicalModels/specie/reaction/reactionRate/JanevReactionRate/JanevReactionRate.H
+++ b/src/thermophysicalModels/specie/reaction/reactionRate/JanevReactionRate/JanevReactionRate.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -103,6 +103,35 @@ public:
             const scalarField& c
         ) const;
 
+        inline scalar ddT
+        (
+            const scalar p,
+            const scalar T,
+            const scalarField& c
+        ) const;
+
+        //- Third-body efficiencies (beta = 1-alpha)
+        //  non-empty only for third-body reactions
+        //  with enhanced molecularity (alpha != 1)
+        inline const List<Tuple2<label, scalar>>& beta() const;
+
+        //- Species concentration derivative of the pressure dependent term
+        inline void dcidc
+        (
+            const scalar p,
+            const scalar T,
+            const scalarField& c,
+            scalarField& dcidc
+        ) const;
+
+        //- Temperature derivative of the pressure dependent term
+        inline scalar dcidT
+        (
+            const scalar p,
+            const scalar T,
+            const scalarField& c
+        ) const;
+
         //- Write to stream
         inline void write(Ostream& os) const;
 
diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/JanevReactionRate/JanevReactionRateI.H b/src/thermophysicalModels/specie/reaction/reactionRate/JanevReactionRate/JanevReactionRateI.H
index 636c3b0c2..cdfb3dc73 100644
--- a/src/thermophysicalModels/specie/reaction/reactionRate/JanevReactionRate/JanevReactionRateI.H
+++ b/src/thermophysicalModels/specie/reaction/reactionRate/JanevReactionRate/JanevReactionRateI.H
@@ -69,7 +69,7 @@ inline Foam::scalar Foam::JanevReactionRate::operator()
         lta *= pow(T, beta_);
     }
 
-    scalar expArg = 0.0;
+    scalar expArg = 0;
 
     if (mag(Ta_) > vSmall)
     {
@@ -89,6 +89,75 @@ inline Foam::scalar Foam::JanevReactionRate::operator()
 }
 
 
+inline Foam::scalar Foam::JanevReactionRate::ddT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField&
+) const
+{
+    scalar lta = A_;
+
+    if (mag(beta_) > vSmall)
+    {
+        lta *= pow(T, beta_);
+    }
+
+    scalar expArg = 0;
+
+    if (mag(Ta_) > vSmall)
+    {
+        expArg -= Ta_/T;
+    }
+
+    scalar lnT = log(T);
+
+    for (int n=0; n<nb_; n++)
+    {
+        expArg += b_[n]*pow(lnT, n);
+    }
+
+    scalar deriv = b_[1];
+
+    for (int n=2; n<nb_; n++)
+    {
+        deriv += n*b_[n]*pow(lnT, n-1);
+    }
+
+    lta *= exp(expArg);
+
+    return lta*(beta_+Ta_/T+deriv)/T;
+}
+
+
+inline const Foam::List<Foam::Tuple2<Foam::label, Foam::scalar>>&
+Foam::JanevReactionRate::beta() const
+{
+    return NullObjectRef<List<Tuple2<label, scalar>>>();
+}
+
+
+inline void Foam::JanevReactionRate::dcidc
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c,
+    scalarField& dcidc
+) const
+{}
+
+
+inline Foam::scalar Foam::JanevReactionRate::dcidT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c
+) const
+{
+    return 0;
+}
+
+
 inline void Foam::JanevReactionRate::write(Ostream& os) const
 {
     os.writeKeyword("A") << A_ << nl;
diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/LandauTellerReactionRate/LandauTellerReactionRate.H b/src/thermophysicalModels/specie/reaction/reactionRate/LandauTellerReactionRate/LandauTellerReactionRate.H
index c353475be..8b1c3572e 100644
--- a/src/thermophysicalModels/specie/reaction/reactionRate/LandauTellerReactionRate/LandauTellerReactionRate.H
+++ b/src/thermophysicalModels/specie/reaction/reactionRate/LandauTellerReactionRate/LandauTellerReactionRate.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -102,6 +102,35 @@ public:
             const scalarField& c
         ) const;
 
+        inline scalar ddT
+        (
+            const scalar p,
+            const scalar T,
+            const scalarField& c
+        ) const;
+
+        //- Third-body efficiencies (beta = 1-alpha)
+        //  non-empty only for third-body reactions
+        //  with enhanced molecularity (alpha != 1)
+        inline const List<Tuple2<label, scalar>>& beta() const;
+
+        //- Species concentration derivative of the pressure dependent term
+        inline void dcidc
+        (
+            const scalar p,
+            const scalar T,
+            const scalarField& c,
+            scalarField& dcidc
+        ) const;
+
+        //- Temperature derivative of the pressure dependent term
+        inline scalar dcidT
+        (
+            const scalar p,
+            const scalar T,
+            const scalarField& c
+        ) const;
+
         //- Write to stream
         inline void write(Ostream& os) const;
 
diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/LandauTellerReactionRate/LandauTellerReactionRateI.H b/src/thermophysicalModels/specie/reaction/reactionRate/LandauTellerReactionRate/LandauTellerReactionRateI.H
index 4d8e92b55..913dbd215 100644
--- a/src/thermophysicalModels/specie/reaction/reactionRate/LandauTellerReactionRate/LandauTellerReactionRateI.H
+++ b/src/thermophysicalModels/specie/reaction/reactionRate/LandauTellerReactionRate/LandauTellerReactionRateI.H
@@ -72,7 +72,7 @@ inline Foam::scalar Foam::LandauTellerReactionRate::operator()
         lta *= pow(T, beta_);
     }
 
-    scalar expArg = 0.0;
+    scalar expArg = 0;
 
     if (mag(Ta_) > vSmall)
     {
@@ -98,6 +98,81 @@ inline Foam::scalar Foam::LandauTellerReactionRate::operator()
 }
 
 
+inline Foam::scalar Foam::LandauTellerReactionRate::ddT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField&
+) const
+{
+    scalar lta = A_;
+
+    if (mag(beta_) > vSmall)
+    {
+        lta *= pow(T, beta_);
+    }
+
+    scalar expArg = 0;
+    scalar deriv = 0;
+
+    if (mag(Ta_) > vSmall)
+    {
+        scalar TaT = Ta_/T;
+        expArg -= TaT;
+        deriv += TaT;
+    }
+
+    if (mag(B_) > vSmall)
+    {
+        scalar BT = B_/cbrt(T);
+        expArg += BT;
+        deriv -= BT/3;
+    }
+
+    if (mag(C_) > vSmall)
+    {
+        scalar CT = C_/pow(T, 2.0/3.0);
+        expArg += CT;
+        deriv -= 2*CT/3;
+    }
+
+    if (mag(expArg) > vSmall)
+    {
+        lta *= exp(expArg);
+    }
+
+    return lta*(beta_+deriv)/T;
+}
+
+
+inline const Foam::List<Foam::Tuple2<Foam::label, Foam::scalar>>&
+Foam::LandauTellerReactionRate::beta() const
+{
+    return NullObjectRef<List<Tuple2<label, scalar>>>();
+}
+
+
+inline void Foam::LandauTellerReactionRate::dcidc
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c,
+    scalarField& dcidc
+) const
+{}
+
+
+inline Foam::scalar Foam::LandauTellerReactionRate::dcidT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c
+) const
+{
+    return 0;
+}
+
+
 inline void Foam::LandauTellerReactionRate::write(Ostream& os) const
 {
     os.writeKeyword("A") << A_ << token::END_STATEMENT << nl;
diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/LangmuirHinshelwood/LangmuirHinshelwoodReactionRate.H b/src/thermophysicalModels/specie/reaction/reactionRate/LangmuirHinshelwood/LangmuirHinshelwoodReactionRate.H
index 1d1c5dd33..4b38fc897 100644
--- a/src/thermophysicalModels/specie/reaction/reactionRate/LangmuirHinshelwood/LangmuirHinshelwoodReactionRate.H
+++ b/src/thermophysicalModels/specie/reaction/reactionRate/LangmuirHinshelwood/LangmuirHinshelwoodReactionRate.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -104,6 +104,35 @@ public:
             const scalarField& c
         ) const;
 
+        inline scalar ddT
+        (
+            const scalar p,
+            const scalar T,
+            const scalarField& c
+        ) const;
+
+        //- Third-body efficiencies (beta = 1-alpha)
+        //  non-empty only for third-body reactions
+        //  with enhanced molecularity (alpha != 1)
+        inline const List<Tuple2<label, scalar>>& beta() const;
+
+        //- Species concentration derivative of the pressure dependent term
+        inline void dcidc
+        (
+            const scalar p,
+            const scalar T,
+            const scalarField& c,
+            scalarField& dcidc
+        ) const;
+
+        //- Temperature derivative of the pressure dependent term
+        inline scalar dcidT
+        (
+            const scalar p,
+            const scalar T,
+            const scalarField& c
+        ) const;
+
         //- Write to stream
         inline void write(Ostream& os) const;
 
diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/LangmuirHinshelwood/LangmuirHinshelwoodReactionRateI.H b/src/thermophysicalModels/specie/reaction/reactionRate/LangmuirHinshelwood/LangmuirHinshelwoodReactionRateI.H
index d0b13587b..7429b2d73 100644
--- a/src/thermophysicalModels/specie/reaction/reactionRate/LangmuirHinshelwood/LangmuirHinshelwoodReactionRateI.H
+++ b/src/thermophysicalModels/specie/reaction/reactionRate/LangmuirHinshelwood/LangmuirHinshelwoodReactionRateI.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -89,6 +89,77 @@ inline Foam::scalar Foam::LangmuirHinshelwoodReactionRate::operator()
 }
 
 
+inline Foam::scalar Foam::LangmuirHinshelwoodReactionRate::ddT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c
+) const
+{
+    scalar den =
+    (
+        T
+       *sqr(1 + A_[1]*exp(-Ta_[1]/T)*c[co_] + A_[2]*exp(-Ta_[2]/T)*c[c3h6_])
+       *(1 + A_[3]*exp(-Ta_[3]/T)*sqr(c[co_])*sqr(c[c3h6_]))
+       *(1 + A_[4]*exp(-Ta_[4]/T)*pow(c[no_], 0.7))
+    );
+    scalar rate = A_[0]*exp(-Ta_[0]/T)/ den;
+
+    scalar derivDen =
+    (
+        sqr(1 + A_[1]*exp(-Ta_[1]/T)*c[co_] + A_[2]*exp(-Ta_[2]/T)*c[c3h6_])
+       *(1 + A_[3]*exp(-Ta_[3]/T)*sqr(c[co_])*sqr(c[c3h6_]))
+       *(1 + A_[4]*exp(-Ta_[4]/T)*pow(c[no_], 0.7))
+      + 2*T*(1 + A_[1]*exp(-Ta_[1]/T)*c[co_] + A_[2]*exp(-Ta_[2]/T)*c[c3h6_])
+       *
+        (
+            A_[1]*exp(-Ta_[1]/T)*c[co_]*Ta_[1]/sqr(T)
+          + A_[2]*exp(-Ta_[2]/T)*c[c3h6_]*Ta_[2]/sqr(T)
+        )
+       *(1 + A_[3]*exp(-Ta_[3]/T)*sqr(c[co_])*sqr(c[c3h6_]))
+       *(1 + A_[4]*exp(-Ta_[4]/T)*pow(c[no_], 0.7))
+      + T
+       *sqr(1 + A_[1]*exp(-Ta_[1]/T)*c[co_] + A_[2]*exp(-Ta_[2]/T)*c[c3h6_])
+       *(A_[3]*exp(-Ta_[3]/T)*Ta_[3]*sqr(c[co_])*sqr(c[c3h6_])/sqr(T))
+       *(1 + A_[4]*exp(-Ta_[4]/T)*pow(c[no_], 0.7))
+      + T
+       *sqr(1 + A_[1]*exp(-Ta_[1]/T)*c[co_] + A_[2]*exp(-Ta_[2]/T)*c[c3h6_])
+       *(1 + A_[3]*exp(-Ta_[3]/T)*sqr(c[co_])*sqr(c[c3h6_]))
+       *(A_[4]*exp(-Ta_[4]/T)*Ta_[4]*pow(c[no_], 0.7))/sqr(T)
+    );
+
+    return rate*(Ta_[0]/sqr(T) - derivDen/den);
+}
+
+
+inline const Foam::List<Foam::Tuple2<Foam::label, Foam::scalar>>&
+Foam::LangmuirHinshelwoodReactionRate::beta() const
+{
+    return NullObjectRef<List<Tuple2<label, scalar>>>();
+}
+
+
+inline void Foam::LangmuirHinshelwoodReactionRate::dcidc
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c,
+    scalarField& dcidc
+) const
+{}
+
+
+inline Foam::scalar Foam::LangmuirHinshelwoodReactionRate::dcidT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c
+) const
+{
+    return 0;
+}
+
+
 inline void Foam::LangmuirHinshelwoodReactionRate::write(Ostream& os) const
 {
     FixedList<Tuple2<scalar, scalar>, n_> coeffs;
diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/fallOffFunctions/LindemannFallOffFunction/LindemannFallOffFunction.H b/src/thermophysicalModels/specie/reaction/reactionRate/fallOffFunctions/LindemannFallOffFunction/LindemannFallOffFunction.H
index 76c2f3acc..23409455b 100644
--- a/src/thermophysicalModels/specie/reaction/reactionRate/fallOffFunctions/LindemannFallOffFunction/LindemannFallOffFunction.H
+++ b/src/thermophysicalModels/specie/reaction/reactionRate/fallOffFunctions/LindemannFallOffFunction/LindemannFallOffFunction.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -80,6 +80,22 @@ public:
             const scalar Pr
         ) const;
 
+        inline scalar ddT
+        (
+            const scalar Pr,
+            const scalar F,
+            const scalar dPrdT,
+            const scalar T
+        ) const;
+
+        inline scalar ddc
+        (
+            const scalar Pr,
+            const scalar F,
+            const scalar dPrdc,
+            const scalar T
+        ) const;
+
         //- Write to stream
         inline void write(Ostream& os) const;
 
diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/fallOffFunctions/LindemannFallOffFunction/LindemannFallOffFunctionI.H b/src/thermophysicalModels/specie/reaction/reactionRate/fallOffFunctions/LindemannFallOffFunction/LindemannFallOffFunctionI.H
index 8421898b6..3ab557fed 100644
--- a/src/thermophysicalModels/specie/reaction/reactionRate/fallOffFunctions/LindemannFallOffFunction/LindemannFallOffFunctionI.H
+++ b/src/thermophysicalModels/specie/reaction/reactionRate/fallOffFunctions/LindemannFallOffFunction/LindemannFallOffFunctionI.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -48,6 +48,30 @@ inline Foam::scalar Foam::LindemannFallOffFunction::operator()
 }
 
 
+inline Foam::scalar Foam::LindemannFallOffFunction::ddT
+(
+    const scalar,
+    const scalar,
+    const scalar,
+    const scalar
+) const
+{
+    return 0;
+}
+
+
+inline Foam::scalar Foam::LindemannFallOffFunction::ddc
+(
+    const scalar,
+    const scalar,
+    const scalar,
+    const scalar
+) const
+{
+    return 0;
+}
+
+
 inline void Foam::LindemannFallOffFunction::write(Ostream& os) const
 {}
 
diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/fallOffFunctions/SRIFallOffFunction/SRIFallOffFunction.H b/src/thermophysicalModels/specie/reaction/reactionRate/fallOffFunctions/SRIFallOffFunction/SRIFallOffFunction.H
index 2a73db4be..7e099ea0d 100644
--- a/src/thermophysicalModels/specie/reaction/reactionRate/fallOffFunctions/SRIFallOffFunction/SRIFallOffFunction.H
+++ b/src/thermophysicalModels/specie/reaction/reactionRate/fallOffFunctions/SRIFallOffFunction/SRIFallOffFunction.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -91,6 +91,22 @@ public:
             const scalar Pr
         ) const;
 
+        inline scalar ddT
+        (
+            const scalar Pr,
+            const scalar F,
+            const scalar dPrdT,
+            const scalar T
+        ) const;
+
+        inline scalar ddc
+        (
+            const scalar Pr,
+            const scalar F,
+            const scalar dPrdc,
+            const scalar T
+        ) const;
+
         //- Write to stream
         inline void write(Ostream& os) const;
 
diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/fallOffFunctions/SRIFallOffFunction/SRIFallOffFunctionI.H b/src/thermophysicalModels/specie/reaction/reactionRate/fallOffFunctions/SRIFallOffFunction/SRIFallOffFunctionI.H
index 82fdf78c8..58ac126ff 100644
--- a/src/thermophysicalModels/specie/reaction/reactionRate/fallOffFunctions/SRIFallOffFunction/SRIFallOffFunctionI.H
+++ b/src/thermophysicalModels/specie/reaction/reactionRate/fallOffFunctions/SRIFallOffFunction/SRIFallOffFunctionI.H
@@ -60,11 +60,48 @@ inline Foam::scalar Foam::SRIFallOffFunction::operator()
     const scalar Pr
 ) const
 {
-    scalar X = 1.0/(1.0 + sqr(log10(max(Pr, small))));
+    scalar X = 1.0/(1 + sqr(log10(max(Pr, small))));
     return d_*pow(a_*exp(-b_/T) + exp(-T/c_), X)*pow(T, e_);
 }
 
 
+inline Foam::scalar Foam::SRIFallOffFunction::ddT
+(
+    const scalar Pr,
+    const scalar F,
+    const scalar dPrdT,
+    const scalar T
+) const
+{
+    scalar X = 1.0/(1 + sqr(log10(max(Pr, small))));
+    scalar dXdPr = -X*X*2*log10(Pr)/Pr/log(10.0);
+    return
+        F
+       *(
+            e_/T
+          + X
+           *(a_*b_*exp(-b_/T)/sqr(T) - exp(-T/c_)/c_)
+           /(a_*exp(-b_/T) + exp(-T/c_))
+          + dXdPr*dPrdT*log(a_*exp(-b_/T) + exp(-T/c_))
+        );
+}
+
+
+inline Foam::scalar Foam::SRIFallOffFunction::ddc
+(
+    const scalar Pr,
+    const scalar F,
+    const scalar dPrdc,
+    const scalar T
+) const
+{
+    scalar X = 1.0/(1 + sqr(log10(max(Pr, small))));
+    scalar dXdPr = -X*X*2*log10(Pr)/Pr/log(10.0);
+    scalar dXdc = dXdPr*dPrdc;
+    return F*dXdc*log(a_*exp(-b_/T) + exp(-T/c_));
+}
+
+
 inline void Foam::SRIFallOffFunction::write(Ostream& os) const
 {
     os.writeKeyword("a") << a_ << token::END_STATEMENT << nl;
diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/fallOffFunctions/TroeFallOffFunction/TroeFallOffFunction.H b/src/thermophysicalModels/specie/reaction/reactionRate/fallOffFunctions/TroeFallOffFunction/TroeFallOffFunction.H
index 18b39038f..4a9025fa4 100644
--- a/src/thermophysicalModels/specie/reaction/reactionRate/fallOffFunctions/TroeFallOffFunction/TroeFallOffFunction.H
+++ b/src/thermophysicalModels/specie/reaction/reactionRate/fallOffFunctions/TroeFallOffFunction/TroeFallOffFunction.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -91,6 +91,22 @@ public:
             const scalar Pr
         ) const;
 
+        inline scalar ddT
+        (
+            const scalar Pr,
+            const scalar F,
+            const scalar dPrdT,
+            const scalar T
+        ) const;
+
+        inline scalar ddc
+        (
+            const scalar Pr,
+            const scalar F,
+            const scalar dPrdc,
+            const scalar T
+        ) const;
+
         //- Write to stream
         inline void write(Ostream& os) const;
 
diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/fallOffFunctions/TroeFallOffFunction/TroeFallOffFunctionI.H b/src/thermophysicalModels/specie/reaction/reactionRate/fallOffFunctions/TroeFallOffFunction/TroeFallOffFunctionI.H
index 3f009702c..79cc84ed0 100644
--- a/src/thermophysicalModels/specie/reaction/reactionRate/fallOffFunctions/TroeFallOffFunction/TroeFallOffFunctionI.H
+++ b/src/thermophysicalModels/specie/reaction/reactionRate/fallOffFunctions/TroeFallOffFunction/TroeFallOffFunctionI.H
@@ -71,7 +71,102 @@ inline Foam::scalar Foam::TroeFallOffFunction::operator()
     scalar n = 0.75 - 1.27*logFcent;
 
     scalar logPr = log10(max(Pr, small));
-    return pow(10.0, logFcent/(1.0 + sqr((logPr + c)/(n - d*(logPr + c)))));
+    return pow(10, logFcent/(1 + sqr((logPr + c)/(n - d*(logPr + c)))));
+}
+
+
+inline Foam::scalar Foam::TroeFallOffFunction::ddT
+(
+    const scalar Pr,
+    const scalar F,
+    const scalar dPrdT,
+    const scalar T
+) const
+{
+    scalar logPr = log10(max(Pr, small));
+    scalar logTen = log(10.0);
+    scalar Fcent =
+    (
+        max
+        (
+            (1 - alpha_)*exp(-T/Tsss_) + alpha_*exp(-T/Ts_) + exp(-Tss_/T),
+            small
+        )
+    );
+    scalar logFcent = log10(Fcent);
+
+    scalar dFcentdT =
+    (
+        (alpha_ - 1)*exp(-T/Tsss_)/Tsss_
+      - alpha_*exp(-T/Ts_)/Ts_
+      + Tss_*exp(-Tss_/T)/sqr(T)
+    );
+
+    scalar d = 0.14;
+    scalar dlogFcentdT = dFcentdT/Fcent/logTen;
+    scalar c = -0.4 - 0.67*logFcent;
+    scalar dcdT = -0.67*dlogFcentdT;
+    scalar n = 0.75 - 1.27*logFcent;
+    scalar dndT = -1.27*dlogFcentdT;
+
+    scalar dlogPrdT = dPrdT/Pr/logTen;
+
+    scalar dParentdT =
+        2.0*(logPr + c)/sqr(n - d*(logPr + c))
+       *(
+            (dlogPrdT + dcdT)
+          - (logPr + c)*(dndT - d*(dlogPrdT + dcdT))/(n - d*(logPr + c))
+        );
+
+    return
+    (
+        F*logTen
+       *(
+            dlogFcentdT/(1.0 + sqr((logPr + c)/(n - d*(logPr + c))))
+          - logFcent*dParentdT/sqr(1.0 + sqr((logPr + c)/(n - d*(logPr + c))))
+        )
+    );
+}
+
+
+inline Foam::scalar Foam::TroeFallOffFunction::ddc
+(
+    const scalar Pr,
+    const scalar F,
+    const scalar dPrdc,
+    const scalar T
+) const
+{
+    scalar logPr = log10(max(Pr, small));
+    scalar logTen = log(10.0);
+    scalar logFcent = log10
+    (
+        max
+        (
+            (1 - alpha_)*exp(-T/Tsss_) + alpha_*exp(-T/Ts_) + exp(-Tss_/T),
+            small
+        )
+    );
+
+    scalar dlogPrdc = dPrdc/Pr/logTen;
+    scalar d = 0.14;
+    scalar c = -0.4 - 0.67*logFcent;
+    scalar n = 0.75 - 1.27*logFcent;
+
+    scalar dParentdc =
+        2.0*(logPr + c)/sqr(n - d*(logPr + c))
+       *(
+            (dlogPrdc)
+          - (logPr + c)*(-d*(dlogPrdc))/(n - d*(logPr + c))
+        );
+
+    return
+    (
+        F*logTen
+       *(
+          - logFcent*dParentdc/sqr(1.0 + sqr((logPr + c)/(n - d*(logPr + c))))
+        )
+    );
 }
 
 
diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/infiniteReactionRate/infiniteReactionRate.H b/src/thermophysicalModels/specie/reaction/reactionRate/infiniteReactionRate/infiniteReactionRate.H
index ec8abe2cd..0dc9f13b2 100644
--- a/src/thermophysicalModels/specie/reaction/reactionRate/infiniteReactionRate/infiniteReactionRate.H
+++ b/src/thermophysicalModels/specie/reaction/reactionRate/infiniteReactionRate/infiniteReactionRate.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -87,6 +87,35 @@ public:
             const scalarField& c
         ) const;
 
+        inline scalar ddT
+        (
+            const scalar p,
+            const scalar T,
+            const scalarField& c
+        ) const;
+
+        //- Third-body efficiencies (beta = 1-alpha)
+        //  non-empty only for third-body reactions
+        //  with enhanced molecularity (alpha != 1)
+        inline const List<Tuple2<label, scalar>>& beta() const;
+
+        //- Species concentration derivative of the pressure dependent term
+        inline void dcidc
+        (
+            const scalar p,
+            const scalar T,
+            const scalarField& c,
+            scalarField& dcidc
+        ) const;
+
+        //- Temperature derivative of the pressure dependent term
+        inline scalar dcidT
+        (
+            const scalar p,
+            const scalar T,
+            const scalarField& c
+        ) const;
+
         //- Write to stream
         inline void write(Ostream& os) const;
 
diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/infiniteReactionRate/infiniteReactionRateI.H b/src/thermophysicalModels/specie/reaction/reactionRate/infiniteReactionRate/infiniteReactionRateI.H
index 918603451..137cf0ec9 100644
--- a/src/thermophysicalModels/specie/reaction/reactionRate/infiniteReactionRate/infiniteReactionRateI.H
+++ b/src/thermophysicalModels/specie/reaction/reactionRate/infiniteReactionRate/infiniteReactionRateI.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -53,6 +53,44 @@ inline Foam::scalar Foam::infiniteReactionRate::operator()
     return (1);
 }
 
+inline Foam::scalar Foam::infiniteReactionRate::ddT
+(
+    const scalar p,
+    const scalar,
+    const scalarField&
+) const
+{
+    return (0);
+}
+
+
+inline const Foam::List<Foam::Tuple2<Foam::label, Foam::scalar>>&
+Foam::infiniteReactionRate::beta() const
+{
+    return NullObjectRef<List<Tuple2<label, scalar>>>();
+}
+
+
+inline void Foam::infiniteReactionRate::dcidc
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c,
+    scalarField& dcidc
+) const
+{}
+
+
+inline Foam::scalar Foam::infiniteReactionRate::dcidT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c
+) const
+{
+    return 0;
+}
+
 
 inline Foam::Ostream& Foam::operator<<
 (
diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/powerSeries/powerSeriesReactionRate.H b/src/thermophysicalModels/specie/reaction/reactionRate/powerSeries/powerSeriesReactionRate.H
index f12d3ba0d..4187365aa 100644
--- a/src/thermophysicalModels/specie/reaction/reactionRate/powerSeries/powerSeriesReactionRate.H
+++ b/src/thermophysicalModels/specie/reaction/reactionRate/powerSeries/powerSeriesReactionRate.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -103,6 +103,35 @@ public:
             const scalarField& c
         ) const;
 
+        inline scalar ddT
+        (
+            const scalar p,
+            const scalar T,
+            const scalarField& c
+        ) const;
+
+        //- Third-body efficiencies (beta = 1-alpha)
+        //  non-empty only for third-body reactions
+        //  with enhanced molecularity (alpha != 1)
+        inline const List<Tuple2<label, scalar>>& beta() const;
+
+        //- Species concentration derivative of the pressure dependent term
+        inline void dcidc
+        (
+            const scalar p,
+            const scalar T,
+            const scalarField& c,
+            scalarField& dcidc
+        ) const;
+
+        //- Temperature derivative of the pressure dependent term
+        inline scalar dcidT
+        (
+            const scalar p,
+            const scalar T,
+            const scalarField& c
+        ) const;
+
         //- Write to stream
         inline void write(Ostream& os) const;
 
diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/powerSeries/powerSeriesReactionRateI.H b/src/thermophysicalModels/specie/reaction/reactionRate/powerSeries/powerSeriesReactionRateI.H
index f727c8866..446b29886 100644
--- a/src/thermophysicalModels/specie/reaction/reactionRate/powerSeries/powerSeriesReactionRateI.H
+++ b/src/thermophysicalModels/specie/reaction/reactionRate/powerSeries/powerSeriesReactionRateI.H
@@ -69,7 +69,7 @@ inline Foam::scalar Foam::powerSeriesReactionRate::operator()
         lta *= pow(T, beta_);
     }
 
-    scalar expArg = 0.0;
+    scalar expArg = 0;
 
     forAll(coeffs_, n)
     {
@@ -82,6 +82,64 @@ inline Foam::scalar Foam::powerSeriesReactionRate::operator()
 }
 
 
+inline Foam::scalar Foam::powerSeriesReactionRate::ddT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField&
+) const
+{
+    scalar lta = A_;
+
+    if (mag(beta_) > vSmall)
+    {
+        lta *= pow(T, beta_);
+    }
+
+    scalar expArg = 0;
+    scalar deriv = 0;
+
+    forAll(coeffs_, n)
+    {
+        scalar cT = coeffs_[n]/pow(T, n + 1);
+        expArg += cT;
+        deriv -= (n + 1)*cT;
+    }
+
+    lta *= exp(expArg);
+
+    return lta*(beta_+deriv)/T;
+}
+
+
+inline const Foam::List<Foam::Tuple2<Foam::label, Foam::scalar>>&
+Foam::powerSeriesReactionRate::beta() const
+{
+    return NullObjectRef<List<Tuple2<label, scalar>>>();
+}
+
+
+inline void Foam::powerSeriesReactionRate::dcidc
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c,
+    scalarField& dcidc
+) const
+{}
+
+
+inline Foam::scalar Foam::powerSeriesReactionRate::dcidT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c
+) const
+{
+    return 0;
+}
+
+
 inline void Foam::powerSeriesReactionRate::write(Ostream& os) const
 {
     os.writeKeyword("A") << A_ << token::END_STATEMENT << nl;
diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/thirdBodyArrheniusReactionRate/thirdBodyArrheniusReactionRate.H b/src/thermophysicalModels/specie/reaction/reactionRate/thirdBodyArrheniusReactionRate/thirdBodyArrheniusReactionRate.H
index 3c7375d73..9e7bae406 100644
--- a/src/thermophysicalModels/specie/reaction/reactionRate/thirdBodyArrheniusReactionRate/thirdBodyArrheniusReactionRate.H
+++ b/src/thermophysicalModels/specie/reaction/reactionRate/thirdBodyArrheniusReactionRate/thirdBodyArrheniusReactionRate.H
@@ -61,6 +61,7 @@ class thirdBodyArrheniusReactionRate
     // Private data
 
         thirdBodyEfficiencies thirdBodyEfficiencies_;
+        List<Tuple2<label, scalar>> beta_;
 
 
 public:
@@ -99,6 +100,38 @@ public:
             const scalarField& c
         ) const;
 
+        inline scalar ddT
+        (
+            const scalar p,
+            const scalar T,
+            const scalarField& c
+        ) const;
+
+        //- Third-body efficiencies (beta = 1-alpha)
+        //  non-empty only for third-body reactions
+        //  with enhanced molecularity (alpha != 1)
+        inline const List<Tuple2<label, scalar>>& beta() const
+        {
+            return beta_;
+        }
+
+        //- Species concentration derivative of the pressure dependent term
+        inline void dcidc
+        (
+            const scalar p,
+            const scalar T,
+            const scalarField& c,
+            scalarField& dcidc
+        ) const;
+
+        //- Temperature derivative of the pressure dependent term
+        inline scalar dcidT
+        (
+            const scalar p,
+            const scalar T,
+            const scalarField& c
+        ) const;
+
         //- Write to stream
         inline void write(Ostream& os) const;
 
diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/thirdBodyArrheniusReactionRate/thirdBodyArrheniusReactionRateI.H b/src/thermophysicalModels/specie/reaction/reactionRate/thirdBodyArrheniusReactionRate/thirdBodyArrheniusReactionRateI.H
index d9b39b61a..7061aeea4 100644
--- a/src/thermophysicalModels/specie/reaction/reactionRate/thirdBodyArrheniusReactionRate/thirdBodyArrheniusReactionRateI.H
+++ b/src/thermophysicalModels/specie/reaction/reactionRate/thirdBodyArrheniusReactionRate/thirdBodyArrheniusReactionRateI.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -35,7 +35,12 @@ inline Foam::thirdBodyArrheniusReactionRate::thirdBodyArrheniusReactionRate
 :
     ArrheniusReactionRate(A, beta, Ta),
     thirdBodyEfficiencies_(tbes)
-{}
+{
+    forAll(tbes, i)
+    {
+        beta_.append(Tuple2<label, scalar>(i, tbes[i]));
+    }
+}
 
 
 inline Foam::thirdBodyArrheniusReactionRate::thirdBodyArrheniusReactionRate
@@ -50,7 +55,19 @@ inline Foam::thirdBodyArrheniusReactionRate::thirdBodyArrheniusReactionRate
         dict
     ),
     thirdBodyEfficiencies_(species, dict)
-{}
+{
+    forAll(thirdBodyEfficiencies_, i)
+    {
+        beta_.append
+        (
+            Tuple2<label, scalar>
+            (
+                i,
+                thirdBodyEfficiencies_[i]
+            )
+        );
+    }
+}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
@@ -68,6 +85,46 @@ inline Foam::scalar Foam::thirdBodyArrheniusReactionRate::operator()
 }
 
 
+inline Foam::scalar Foam::thirdBodyArrheniusReactionRate::ddT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c
+) const
+{
+    return
+        thirdBodyEfficiencies_.M(c)
+       *ArrheniusReactionRate::ddT(p, T, c);
+}
+
+
+inline void Foam::thirdBodyArrheniusReactionRate::dcidc
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c,
+    scalarField& dcidc
+) const
+{
+    scalar M = thirdBodyEfficiencies_.M(c);
+    forAll(beta_, i)
+    {
+        dcidc[i] = beta_[i].second()/max(M, small);
+    }
+}
+
+
+inline Foam::scalar Foam::thirdBodyArrheniusReactionRate::dcidT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c
+) const
+{
+    return -1.0/T;
+}
+
+
 inline void Foam::thirdBodyArrheniusReactionRate::write(Ostream& os) const
 {
     ArrheniusReactionRate::write(os);
diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/thirdBodyEfficiencies/thirdBodyEfficienciesI.H b/src/thermophysicalModels/specie/reaction/reactionRate/thirdBodyEfficiencies/thirdBodyEfficienciesI.H
index 97a0e8eee..a1aa8265b 100644
--- a/src/thermophysicalModels/specie/reaction/reactionRate/thirdBodyEfficiencies/thirdBodyEfficienciesI.H
+++ b/src/thermophysicalModels/specie/reaction/reactionRate/thirdBodyEfficiencies/thirdBodyEfficienciesI.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -83,7 +83,7 @@ inline Foam::thirdBodyEfficiencies::thirdBodyEfficiencies
 
 inline Foam::scalar Foam::thirdBodyEfficiencies::M(const scalarList& c) const
 {
-    scalar M = 0.0;
+    scalar M = 0;
     forAll(*this, i)
     {
         M += operator[](i)*c[i];
diff --git a/src/thermophysicalModels/specie/thermo/eConst/eConstThermo.H b/src/thermophysicalModels/specie/thermo/eConst/eConstThermo.H
index 2a9c599ff..b212d93ec 100644
--- a/src/thermophysicalModels/specie/thermo/eConst/eConstThermo.H
+++ b/src/thermophysicalModels/specie/thermo/eConst/eConstThermo.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -155,6 +155,16 @@ public:
             inline scalar S(const scalar p, const scalar T) const;
 
 
+        // Derivative term used for Jacobian
+
+            //- Derivative of Gibbs free energy w.r.t. temperature
+            inline scalar dGdT(const scalar p, const scalar T) const;
+
+            //- Temperature derivative of heat capacity at constant pressure
+            inline scalar dCpdT(const scalar p, const scalar T) const;
+
+
+
         // I-O
 
             //- Write to Ostream
diff --git a/src/thermophysicalModels/specie/thermo/eConst/eConstThermoI.H b/src/thermophysicalModels/specie/thermo/eConst/eConstThermoI.H
index f913d48f9..c2522b34f 100644
--- a/src/thermophysicalModels/specie/thermo/eConst/eConstThermoI.H
+++ b/src/thermophysicalModels/specie/thermo/eConst/eConstThermoI.H
@@ -139,6 +139,28 @@ inline Foam::scalar Foam::eConstThermo<EquationOfState>::S
 }
 
 
+template<class EquationOfState>
+inline Foam::scalar Foam::eConstThermo<EquationOfState>::dGdT
+(
+    const scalar p,
+    const scalar T
+) const
+{
+    return 0;
+}
+
+
+template<class EquationOfState>
+inline Foam::scalar Foam::eConstThermo<EquationOfState>::dCpdT
+(
+    const scalar p,
+    const scalar T
+) const
+{
+    return 0;
+}
+
+
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
 template<class EquationOfState>
diff --git a/src/thermophysicalModels/specie/thermo/hConst/hConstThermo.H b/src/thermophysicalModels/specie/thermo/hConst/hConstThermo.H
index dadbb0712..5a5c7626e 100644
--- a/src/thermophysicalModels/specie/thermo/hConst/hConstThermo.H
+++ b/src/thermophysicalModels/specie/thermo/hConst/hConstThermo.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -148,6 +148,15 @@ public:
             inline scalar S(const scalar p, const scalar T) const;
 
 
+        // Derivative term used for Jacobian
+
+            //- Derivative of Gibbs free energy w.r.t. temperature
+            inline scalar dGdT(const scalar p, const scalar T) const;
+
+            //- Temperature derivative of heat capacity at constant pressure
+            inline scalar dCpdT(const scalar p, const scalar T) const;
+
+
         // I-O
 
             //- Write to Ostream
diff --git a/src/thermophysicalModels/specie/thermo/hConst/hConstThermoI.H b/src/thermophysicalModels/specie/thermo/hConst/hConstThermoI.H
index cf75c27e2..5c48335b3 100644
--- a/src/thermophysicalModels/specie/thermo/hConst/hConstThermoI.H
+++ b/src/thermophysicalModels/specie/thermo/hConst/hConstThermoI.H
@@ -136,6 +136,26 @@ inline Foam::scalar Foam::hConstThermo<EquationOfState>::S
 }
 
 
+template<class EquationOfState>
+inline Foam::scalar Foam::hConstThermo<EquationOfState>::dGdT
+(
+    const scalar p, const scalar T
+) const
+{
+    return 0;
+}
+
+
+template<class EquationOfState>
+inline Foam::scalar Foam::hConstThermo<EquationOfState>::dCpdT
+(
+    const scalar p, const scalar T
+) const
+{
+    return 0;
+}
+
+
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
 template<class EquationOfState>
diff --git a/src/thermophysicalModels/specie/thermo/hPolynomial/hPolynomialThermo.H b/src/thermophysicalModels/specie/thermo/hPolynomial/hPolynomialThermo.H
index f9b73bd93..c9b5adb60 100644
--- a/src/thermophysicalModels/specie/thermo/hPolynomial/hPolynomialThermo.H
+++ b/src/thermophysicalModels/specie/thermo/hPolynomial/hPolynomialThermo.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -198,6 +198,15 @@ public:
             inline scalar S(const scalar p, const scalar T) const;
 
 
+        // Derivative term used for Jacobian
+
+            //- Derivative of Gibbs free energy w.r.t. temperature
+            inline scalar dGdT(const scalar p, const scalar T) const;
+
+            //- Temperature derivative of heat capacity at constant pressure
+            inline scalar dCpdT(const scalar p, const scalar T) const;
+
+
         // I-O
 
             //- Write to Ostream
diff --git a/src/thermophysicalModels/specie/thermo/hPolynomial/hPolynomialThermoI.H b/src/thermophysicalModels/specie/thermo/hPolynomial/hPolynomialThermoI.H
index 370b3502e..84a96ed19 100644
--- a/src/thermophysicalModels/specie/thermo/hPolynomial/hPolynomialThermoI.H
+++ b/src/thermophysicalModels/specie/thermo/hPolynomial/hPolynomialThermoI.H
@@ -126,6 +126,36 @@ inline Foam::scalar Foam::hPolynomialThermo<EquationOfState, PolySize>::S
 }
 
 
+template<class EquationOfState, int PolySize>
+inline Foam::scalar Foam::hPolynomialThermo<EquationOfState, PolySize>::dGdT
+(
+    const scalar p,
+    const scalar T
+) const
+{
+    return
+    (
+        hCoeffs_.derivative(T)
+      - T*sCoeffs_.derivative(T)
+      - sCoeffs_.value(T)
+    );
+}
+
+
+template<class EquationOfState, int PolySize>
+inline Foam::scalar Foam::hPolynomialThermo<EquationOfState, PolySize>::dCpdT
+(
+    const scalar p,
+    const scalar T
+) const
+{
+    return
+    (
+        CpCoeffs_.derivative(T)
+    );
+}
+
+
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
 template<class EquationOfState, int PolySize>
diff --git a/src/thermophysicalModels/specie/thermo/hPower/hPowerThermo.H b/src/thermophysicalModels/specie/thermo/hPower/hPowerThermo.H
index 7ca3b2bbd..cf19e37e4 100644
--- a/src/thermophysicalModels/specie/thermo/hPower/hPowerThermo.H
+++ b/src/thermophysicalModels/specie/thermo/hPower/hPowerThermo.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -169,6 +169,15 @@ public:
             inline scalar S(const scalar p, const scalar T) const;
 
 
+        // Derivative term used for Jacobian
+
+            //- Derivative of Gibbs free energy w.r.t. temperature
+            inline scalar dGdT(const scalar p, const scalar T) const;
+
+            //- Temperature derivative of heat capacity at constant pressure
+            inline scalar dCpdT(const scalar p, const scalar T) const;
+
+
     // Member operators
 
         inline void operator+=(const hPowerThermo&);
diff --git a/src/thermophysicalModels/specie/thermo/hPower/hPowerThermoI.H b/src/thermophysicalModels/specie/thermo/hPower/hPowerThermoI.H
index 2b2f61016..115a9af7f 100644
--- a/src/thermophysicalModels/specie/thermo/hPower/hPowerThermoI.H
+++ b/src/thermophysicalModels/specie/thermo/hPower/hPowerThermoI.H
@@ -164,6 +164,30 @@ inline Foam::scalar Foam::hPowerThermo<EquationOfState>::S
 }
 
 
+template<class EquationOfState>
+inline Foam::scalar Foam::hPowerThermo<EquationOfState>::dGdT
+(
+    const scalar p, const scalar T
+) const
+{
+    // To be implemented
+    NotImplemented;
+    return 0;
+}
+
+
+template<class EquationOfState>
+inline Foam::scalar Foam::hPowerThermo<EquationOfState>::dCpdT
+(
+    const scalar p, const scalar T
+) const
+{
+    // To be implemented
+    NotImplemented;
+    return 0;
+}
+
+
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
 template<class EquationOfState>
diff --git a/src/thermophysicalModels/specie/thermo/hRefConst/hRefConstThermo.H b/src/thermophysicalModels/specie/thermo/hRefConst/hRefConstThermo.H
index 5af026540..a139fa564 100644
--- a/src/thermophysicalModels/specie/thermo/hRefConst/hRefConstThermo.H
+++ b/src/thermophysicalModels/specie/thermo/hRefConst/hRefConstThermo.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2015-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2015-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -152,6 +152,15 @@ public:
             inline scalar S(const scalar p, const scalar T) const;
 
 
+        // Derivative term used for Jacobian
+
+            //- Derivative of Gibbs free energy w.r.t. temperature
+            inline scalar dGdT(const scalar p, const scalar T) const;
+
+            //- Temperature derivative of heat capacity at constant pressure
+            inline scalar dCpdT(const scalar p, const scalar T) const;
+
+
         // I-O
 
             //- Write to Ostream
diff --git a/src/thermophysicalModels/specie/thermo/hRefConst/hRefConstThermoI.H b/src/thermophysicalModels/specie/thermo/hRefConst/hRefConstThermoI.H
index e5ebfd954..15e7988ca 100644
--- a/src/thermophysicalModels/specie/thermo/hRefConst/hRefConstThermoI.H
+++ b/src/thermophysicalModels/specie/thermo/hRefConst/hRefConstThermoI.H
@@ -142,6 +142,26 @@ inline Foam::scalar Foam::hRefConstThermo<EquationOfState>::S
 }
 
 
+template<class EquationOfState>
+inline Foam::scalar Foam::hRefConstThermo<EquationOfState>::dGdT
+(
+    const scalar p, const scalar T
+) const
+{
+    return 0;
+}
+
+
+template<class EquationOfState>
+inline Foam::scalar Foam::hRefConstThermo<EquationOfState>::dCpdT
+(
+    const scalar p, const scalar T
+) const
+{
+    return 0;
+}
+
+
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
 template<class EquationOfState>
diff --git a/src/thermophysicalModels/specie/thermo/janaf/janafThermo.H b/src/thermophysicalModels/specie/thermo/janaf/janafThermo.H
index f16b76ebf..8b487c38c 100644
--- a/src/thermophysicalModels/specie/thermo/janaf/janafThermo.H
+++ b/src/thermophysicalModels/specie/thermo/janaf/janafThermo.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -187,6 +187,15 @@ public:
             inline scalar S(const scalar p, const scalar T) const;
 
 
+        // Derivative term used for Jacobian
+
+            //- Derivative of Gibbs free energy w.r.t. temperature
+            inline scalar dGdT(const scalar p, const scalar T) const;
+
+            //- Temperature derivative of heat capacity at constant pressure
+            inline scalar dCpdT(const scalar p, const scalar T) const;
+
+
         // I-O
 
             //- Write to Ostream
diff --git a/src/thermophysicalModels/specie/thermo/janaf/janafThermoI.H b/src/thermophysicalModels/specie/thermo/janaf/janafThermoI.H
index ad8b1a9fa..7391c4b30 100644
--- a/src/thermophysicalModels/specie/thermo/janaf/janafThermoI.H
+++ b/src/thermophysicalModels/specie/thermo/janaf/janafThermoI.H
@@ -237,6 +237,31 @@ inline Foam::scalar Foam::janafThermo<EquationOfState>::S
 }
 
 
+template<class EquationOfState>
+inline Foam::scalar Foam::janafThermo<EquationOfState>::dGdT
+(
+    const scalar p,
+    const scalar T
+) const
+{
+    const coeffArray& a = coeffs(T);
+    return -((a[0] + a[5]/T)/T + a[1]/2 + T*(a[2]/3 + T*(a[3]/4 + T*a[4]/5)));
+}
+
+
+template<class EquationOfState>
+inline Foam::scalar Foam::janafThermo<EquationOfState>::dCpdT
+(
+    const scalar p,
+    const scalar T
+) const
+{
+    const coeffArray& a = coeffs(T);
+    return
+        (((4*a[4]*T + 3*a[3])*T + 2*a[2])*T + a[1]);
+}
+
+
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
 template<class EquationOfState>
diff --git a/src/thermophysicalModels/specie/thermo/thermo/thermo.H b/src/thermophysicalModels/specie/thermo/thermo/thermo.H
index a5cdba6ec..3a644d5d8 100644
--- a/src/thermophysicalModels/specie/thermo/thermo/thermo.H
+++ b/src/thermophysicalModels/specie/thermo/thermo/thermo.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -22,7 +22,7 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::thermo
+    Foam::species::thermo
 
 Description
     Basic thermodynamics type based on the use of fitting functions for
@@ -318,6 +318,15 @@ public:
             ) const;
 
 
+        // Derivative term used for Jacobian
+
+            //- Derivative of B (acooding to Niemeyer et al.) w.r.t. temperature
+            inline scalar dKcdTbyKc(const scalar p, const scalar T) const;
+
+            //- Derivative of cp w.r.t. temperature
+            inline scalar dcpdT(const scalar p, const scalar T) const;
+
+
         // I-O
 
             //- Write to Ostream
diff --git a/src/thermophysicalModels/specie/thermo/thermo/thermoI.H b/src/thermophysicalModels/specie/thermo/thermo/thermoI.H
index 791042fe7..0fced0cee 100644
--- a/src/thermophysicalModels/specie/thermo/thermo/thermoI.H
+++ b/src/thermophysicalModels/specie/thermo/thermo/thermoI.H
@@ -444,6 +444,35 @@ inline Foam::scalar Foam::species::thermo<Thermo, Type>::TEa
 }
 
 
+template<class Thermo, template<class> class Type>
+inline Foam::scalar
+Foam::species::thermo<Thermo, Type>::dKcdTbyKc
+(
+    const scalar p,
+    const scalar T
+) const
+{
+    const scalar nm = this->Y()/this->W();
+
+    if (equal(nm, small))
+    {
+        return -this->dGdT(Pstd, T)*this->Y()/RR;
+    }
+    else
+    {
+        return -(nm/T + this->dGdT(Pstd, T)*this->Y()/RR);
+    }
+}
+
+
+template<class Thermo, template<class> class Type>
+inline Foam::scalar
+Foam::species::thermo<Thermo, Type>::dcpdT(const scalar p, const scalar T) const
+{
+    return this->dCpdT(p, T)*this->W();;
+}
+
+
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
 template<class Thermo, template<class> class Type>
diff --git a/tutorials/combustion/chemFoam/nc7h16/constant/chemistryProperties b/tutorials/combustion/chemFoam/nc7h16/constant/chemistryProperties
index 2dc50a724..9c75d7dff 100644
--- a/tutorials/combustion/chemFoam/nc7h16/constant/chemistryProperties
+++ b/tutorials/combustion/chemFoam/nc7h16/constant/chemistryProperties
@@ -33,7 +33,7 @@ EulerImplicitCoeffs
 odeCoeffs
 {
     solver          seulex;
-    absTol          1e-14;
+    absTol          1e-12;
     relTol          1e-1;
 }
 
diff --git a/tutorials/combustion/fireFoam/LES/oppositeBurningPanels/constant/panelRegion/chemistryProperties b/tutorials/combustion/fireFoam/LES/oppositeBurningPanels/constant/panelRegion/chemistryProperties
index df407c07f..46997d9d2 100644
--- a/tutorials/combustion/fireFoam/LES/oppositeBurningPanels/constant/panelRegion/chemistryProperties
+++ b/tutorials/combustion/fireFoam/LES/oppositeBurningPanels/constant/panelRegion/chemistryProperties
@@ -28,7 +28,7 @@ initialChemicalTimeStep 1e-07;
 odeCoeffs
 {
     solver          seulex;
-    absTol          1e-12;
+    absTol          1e-8;
     relTol          0.01;
 }
 
diff --git a/tutorials/combustion/reactingFoam/RAS/DLR_A_LTS/constant/chemistryProperties b/tutorials/combustion/reactingFoam/RAS/DLR_A_LTS/constant/chemistryProperties
index c250d3841..aafbdf700 100644
--- a/tutorials/combustion/reactingFoam/RAS/DLR_A_LTS/constant/chemistryProperties
+++ b/tutorials/combustion/reactingFoam/RAS/DLR_A_LTS/constant/chemistryProperties
@@ -35,9 +35,9 @@ initialChemicalTimeStep 1e-07;
 
 odeCoeffs
 {
-    solver          Rosenbrock34;  // Rosenbrock34, seulex or rodas23
-    absTol          1e-12;
-    relTol          1e-7;
+    solver          seulex;
+    absTol          1e-08;
+    relTol          0.1;
 }
 
 reduction
@@ -86,7 +86,7 @@ tabulation
     printNumRetrieve   off;
 
     // Tolerance used for retrieve and grow
-    tolerance   1e-4;
+    tolerance   0.003;
 
     // ISAT is the only method currently available
     method    ISAT;
@@ -97,24 +97,24 @@ tabulation
         otherSpecies 1;
         Temperature  1000;
         Pressure     1e15;
-        deltaT       0.5;
+        deltaT       1;
     }
 
     // Maximum number of leafs stored in the binary tree
-    maxNLeafs  2000;
+    maxNLeafs  5000;
 
     // Maximum life time of the leafs (in time steps) used in unsteady
     // simulations to force renewal of the stored chemPoints and keep the tree
     // small
-    chPMaxLifeTime 100;
+    chPMaxLifeTime 1000;
 
     // Maximum number of growth allowed on a chemPoint to avoid distorted
     // chemPoints
-    maxGrowth  10;
+    maxGrowth  100;
 
     // Number of time steps between analysis of the tree to remove old
     // chemPoints or try to balance it
-    checkEntireTreeInterval  5;
+    checkEntireTreeInterval  500;
 
     // Parameters used to decide whether to balance or not if the tree's depth
     // is larger than maxDepthFactor*log2(nLeafs) then balance the tree
@@ -123,7 +123,6 @@ tabulation
     // Try to balance the tree only if the size of the tree is greater
     minBalanceThreshold 30;
 
-
     // Activate the use of a MRU (most recently used) list
     MRURetrieve false;
 
diff --git a/tutorials/combustion/reactingFoam/RAS/DLR_A_LTS/constant/chemistryProperties.test b/tutorials/combustion/reactingFoam/RAS/DLR_A_LTS/constant/chemistryProperties.test
deleted file mode 100644
index cb5c5b610..000000000
--- a/tutorials/combustion/reactingFoam/RAS/DLR_A_LTS/constant/chemistryProperties.test
+++ /dev/null
@@ -1,142 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  dev                                   |
-|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       dictionary;
-    location    "constant";
-    object      chemistryProperties;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-chemistryType
-{
-    solver            ode;
-    method            TDAC;
-}
-
-chemistry           on;
-
-importantSpecies
-{
-    CH4;
-    H2O;
-    O2;
-    CO2;
-}
-
-initialChemicalTimeStep 1e-07;
-
-odeCoeffs
-{
-    solver          Rosenbrock34;  // Rosenbrock34, seulex or rodas23
-    absTol          1e-12;
-    relTol          1e-7;
-}
-
-reduction
-{
-    // Activate reduction
-    active  on;
-
-    // Switch logging of the reduction statistics and performance
-    log         on;
-
-    // Tolerance depends on the reduction method, see details for each method
-    tolerance   1e-4;
-
-    // Available methods: DRG, DAC, DRGEP, PFA, EFA
-    method DAC;
-
-    // Search initiating set (SIS) of species, needed for most methods
-    initialSet
-    {
-        CO;
-        CH4;
-        HO2;
-    }
-
-    // For DAC, option to automatically change the SIS switch from HO2 to H2O
-    // and CO to CO2, + disable fuel
-    automaticSIS    off;
-
-    // When automaticSIS, the method needs to know the fuel
-    fuelSpecies
-    {
-        CH4 1;
-    }
-}
-
-tabulation
-{
-    // Activate tabulation
-    active      on;
-
-    // Switch logging of the tabulation statistics and performance
-    log         on;
-
-    printProportion    off;
-
-    printNumRetrieve   off;
-
-    // Tolerance used for retrieve and grow
-    tolerance   0.003;
-
-    // ISAT is the only method currently available
-    method    ISAT;
-
-    // Scale factors used in the definition of the ellipsoid of accuracy
-    scaleFactor
-    {
-        otherSpecies 1;
-        Temperature  10000;
-        Pressure     1e15;
-        deltaT       1;
-    }
-
-    // Maximum number of leafs stored in the binary tree
-    maxNLeafs  5000;
-
-    // Maximum life time of the leafs (in time steps) used in unsteady
-    // simulations to force renewal of the stored chemPoints and keep the tree
-    // small
-    chPMaxLifeTime 1000;
-
-    // Maximum number of growth allowed on a chemPoint to avoid distorted
-    // chemPoints
-    maxGrowth  100;
-
-    // Number of time steps between analysis of the tree to remove old
-    // chemPoints or try to balance it
-    checkEntireTreeInterval  500;
-
-    // Parameters used to decide whether to balance or not if the tree's depth
-    // is larger than maxDepthFactor*log2(nLeafs) then balance the tree
-    maxDepthFactor   2;
-
-    // Try to balance the tree only if the size of the tree is greater
-    minBalanceThreshold 30;
-
-    // Activate the use of a MRU (most recently used) list
-    MRURetrieve false;
-
-    // Maximum size of the MRU list
-    maxMRUSize 0;
-
-    // Allow to grow points
-    growPoints  true;
-
-    // When mechanism reduction is used, new dimensions might be added
-    // maxNumNewDim set the maximum number of new dimensions added during a
-    // growth
-    maxNumNewDim 10;
-}
-
-
-// ************************************************************************* //
diff --git a/tutorials/combustion/reactingFoam/RAS/DLR_A_LTS/constant/reactionsGRI b/tutorials/combustion/reactingFoam/RAS/DLR_A_LTS/constant/reactionsGRI
index 8bd6be358..9e89ab5d9 100644
--- a/tutorials/combustion/reactingFoam/RAS/DLR_A_LTS/constant/reactionsGRI
+++ b/tutorials/combustion/reactingFoam/RAS/DLR_A_LTS/constant/reactionsGRI
@@ -3608,3 +3608,6 @@ reactions
         Ta              0;
     }
 }
+
+Tlow            250;
+Thigh           5000;
diff --git a/tutorials/combustion/reactingFoam/RAS/SandiaD_LTS/Allrun b/tutorials/combustion/reactingFoam/RAS/SandiaD_LTS/Allrun
index 2769119da..d8137ce0a 100755
--- a/tutorials/combustion/reactingFoam/RAS/SandiaD_LTS/Allrun
+++ b/tutorials/combustion/reactingFoam/RAS/SandiaD_LTS/Allrun
@@ -14,16 +14,16 @@ runApplication chemkinToFoam \
 runApplication blockMesh
 runApplication setFields
 
-
 # Run the application without chemistry until 1500 to let the flow field develop
+foamDictionary -entry "startTime" -set "0" system/controlDict
 foamDictionary -entry "writeInterval" -set "1500" system/controlDict
 foamDictionary -entry "endTime" -set "1500" system/controlDict
 foamDictionary -entry "chemistry" -set "off" constant/chemistryProperties
 
 runApplication $application
 
-
 # Run with chemistry until flame reach its full size
+foamDictionary -entry "startTime" -set "1500" system/controlDict
 foamDictionary -entry "writeInterval" -set "100" system/controlDict
 foamDictionary -entry "endTime" -set "5000" system/controlDict
 foamDictionary -entry "chemistry" -set "on" constant/chemistryProperties
diff --git a/tutorials/combustion/reactingFoam/RAS/SandiaD_LTS/constant/chemistryProperties b/tutorials/combustion/reactingFoam/RAS/SandiaD_LTS/constant/chemistryProperties
index f8c1f4c3d..cf015f487 100644
--- a/tutorials/combustion/reactingFoam/RAS/SandiaD_LTS/constant/chemistryProperties
+++ b/tutorials/combustion/reactingFoam/RAS/SandiaD_LTS/constant/chemistryProperties
@@ -20,7 +20,7 @@ chemistryType
     method          TDAC;
 }
 
-chemistry       off;
+chemistry       on;
 
 importantSpecies
 {
@@ -35,8 +35,8 @@ initialChemicalTimeStep 1e-07;
 odeCoeffs
 {
     solver          seulex;
-    absTol          1e-12;
-    relTol          1e-07;
+    absTol          1e-08;
+    relTol          0.1;
 }
 
 reduction
diff --git a/tutorials/combustion/reactingFoam/RAS/SandiaD_LTS/constant/reactionsGRI b/tutorials/combustion/reactingFoam/RAS/SandiaD_LTS/constant/reactionsGRI
index 8bd6be358..9e89ab5d9 100644
--- a/tutorials/combustion/reactingFoam/RAS/SandiaD_LTS/constant/reactionsGRI
+++ b/tutorials/combustion/reactingFoam/RAS/SandiaD_LTS/constant/reactionsGRI
@@ -3608,3 +3608,6 @@ reactions
         Ta              0;
     }
 }
+
+Tlow            250;
+Thigh           5000;
diff --git a/tutorials/combustion/reactingFoam/RAS/SandiaD_LTS/system/controlDict b/tutorials/combustion/reactingFoam/RAS/SandiaD_LTS/system/controlDict
index bf9887f7b..774eb04d0 100644
--- a/tutorials/combustion/reactingFoam/RAS/SandiaD_LTS/system/controlDict
+++ b/tutorials/combustion/reactingFoam/RAS/SandiaD_LTS/system/controlDict
@@ -16,19 +16,19 @@ FoamFile
 
 application     reactingFoam;
 
-startFrom       latestTime;
+startFrom       startTime;
 
-startTime       0;
+startTime       1500;
 
 stopAt          endTime;
 
-endTime         1500;
+endTime         5000;
 
 deltaT          1;
 
 writeControl    runTime;
 
-writeInterval   1500;
+writeInterval   100;
 
 purgeWrite      0;
 
diff --git a/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/constant/chemistryProperties b/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/constant/chemistryProperties
index 3812e6131..f49edd526 100644
--- a/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/constant/chemistryProperties
+++ b/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/constant/chemistryProperties
@@ -33,7 +33,7 @@ EulerImplicitCoeffs
 odeCoeffs
 {
     solver          Rosenbrock43;
-    absTol          1e-12;
+    absTol          1e-8;
     relTol          0.01;
 }
 
diff --git a/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2DLTS/constant/chemistryProperties b/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2DLTS/constant/chemistryProperties
index 3812e6131..f49edd526 100644
--- a/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2DLTS/constant/chemistryProperties
+++ b/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2DLTS/constant/chemistryProperties
@@ -33,7 +33,7 @@ EulerImplicitCoeffs
 odeCoeffs
 {
     solver          Rosenbrock43;
-    absTol          1e-12;
+    absTol          1e-8;
     relTol          0.01;
 }
 
diff --git a/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2DLTS_GRI_TDAC/constant/chemistryProperties b/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2DLTS_GRI_TDAC/constant/chemistryProperties
index 2604f1fd8..9fd9f1912 100644
--- a/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2DLTS_GRI_TDAC/constant/chemistryProperties
+++ b/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2DLTS_GRI_TDAC/constant/chemistryProperties
@@ -24,11 +24,12 @@ chemistryType
 chemistry       on;
 
 initialChemicalTimeStep 1e-7;
+//maxChemicalTimeStep 1E-3;
 
 odeCoeffs
 {
     solver          seulex;
-    absTol          1e-12;
+    absTol          1e-8;
     relTol          1e-1;
 }
 
@@ -78,7 +79,7 @@ tabulation
     printNumRetrieve   off;
 
     // Tolerance used for retrieve and grow
-    tolerance   1e-3;
+    tolerance   3e-3;
 
     // ISAT is the only method currently available
     method    ISAT;
@@ -87,7 +88,7 @@ tabulation
     scaleFactor
     {
         otherSpecies 1;
-        Temperature  25000;
+        Temperature  10000;
         Pressure     1e15;
         deltaT       1;
     }
diff --git a/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2DLTS_GRI_TDAC/constant/reactionsGRI b/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2DLTS_GRI_TDAC/constant/reactionsGRI
index de8ff0a51..694ea19b6 100644
--- a/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2DLTS_GRI_TDAC/constant/reactionsGRI
+++ b/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2DLTS_GRI_TDAC/constant/reactionsGRI
@@ -5588,3 +5588,6 @@ reactions
         }
     }
 }
+
+Thigh   5000;
+Tlow    200;
diff --git a/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2DLTS_GRI_TDAC/system/controlDict b/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2DLTS_GRI_TDAC/system/controlDict
index 792dbf249..0aff953b0 100644
--- a/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2DLTS_GRI_TDAC/system/controlDict
+++ b/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2DLTS_GRI_TDAC/system/controlDict
@@ -23,7 +23,7 @@ startTime       0;
 
 stopAt          endTime;
 
-endTime         1000;
+endTime         1500;
 
 deltaT          1;
 
diff --git a/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D_GRI/constant/chemistryProperties b/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D_GRI/constant/chemistryProperties
index d63f983b7..678d5e60f 100644
--- a/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D_GRI/constant/chemistryProperties
+++ b/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D_GRI/constant/chemistryProperties
@@ -7,29 +7,28 @@
 \*---------------------------------------------------------------------------*/
 FoamFile
 {
-    version     2.0;
-    format      ascii;
-    class       dictionary;
-    location    "constant";
-    object      chemistryProperties;
+    version         2;
+    format          ascii;
+    class           dictionary;
+    location        "constant";
+    object          chemistryProperties;
 }
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 chemistryType
 {
-    solver            ode;
+    solver          ode;
 }
 
 chemistry       on;
 
-initialChemicalTimeStep 1e-7;
-maxChemicalTimeStep     1e-4;
+initialChemicalTimeStep 1e-07;
 
 odeCoeffs
 {
     solver          seulex;
-    absTol          1e-12;
-    relTol          1e-1;
+    absTol          1e-8;
+    relTol          0.1;
 }
 
+
 // ************************************************************************* //
diff --git a/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D_GRI/constant/reactionsGRI b/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D_GRI/constant/reactionsGRI
index 276111048..694ea19b6 100644
--- a/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D_GRI/constant/reactionsGRI
+++ b/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D_GRI/constant/reactionsGRI
@@ -2349,7 +2349,7 @@ reactions
         A               1.2e+11;
         beta            -1;
         Ta              0;
-        coeffs          
+        coeffs
 53
 (
 (CH4 2)
@@ -2415,7 +2415,7 @@ reactions
         A               5e+11;
         beta            -1;
         Ta              0;
-        coeffs          
+        coeffs
 53
 (
 (CH4 2)
@@ -2481,7 +2481,7 @@ reactions
         A               2.8e+12;
         beta            -0.86;
         Ta              0;
-        coeffs          
+        coeffs
 53
 (
 (CH4 1)
@@ -2547,7 +2547,7 @@ reactions
         A               1e+12;
         beta            -1;
         Ta              0;
-        coeffs          
+        coeffs
 53
 (
 (CH4 2)
@@ -2613,7 +2613,7 @@ reactions
         A               2.2e+16;
         beta            -2;
         Ta              0;
-        coeffs          
+        coeffs
 53
 (
 (CH4 2)
@@ -2679,7 +2679,7 @@ reactions
         A               1.87e+14;
         beta            -1;
         Ta              8554.25;
-        coeffs          
+        coeffs
 53
 (
 (CH4 2)
@@ -2745,7 +2745,7 @@ reactions
         A               1.06e+14;
         beta            -1.41;
         Ta              0;
-        coeffs          
+        coeffs
 53
 (
 (CH4 2)
@@ -2811,7 +2811,7 @@ reactions
         A               1.3e+11;
         beta            -0.11;
         Ta              2505.89;
-        coeffs          
+        coeffs
 53
 (
 (CH4 2)
@@ -2877,7 +2877,7 @@ reactions
         A               4.48e+13;
         beta            -1.32;
         Ta              372.362;
-        coeffs          
+        coeffs
 53
 (
 (CH4 2)
@@ -2943,7 +2943,7 @@ reactions
         A               3.1e+11;
         beta            0;
         Ta              27197.5;
-        coeffs          
+        coeffs
 53
 (
 (CH4 2)
@@ -3009,7 +3009,7 @@ reactions
         A               1.04e+26;
         beta            -3.3;
         Ta              63704;
-        coeffs          
+        coeffs
 53
 (
 (CH4 2)
@@ -3075,7 +3075,7 @@ reactions
         A               1.18e+13;
         beta            0;
         Ta              42630.4;
-        coeffs          
+        coeffs
 53
 (
 (CH4 2)
@@ -3155,7 +3155,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -3240,7 +3240,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -3325,7 +3325,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 3)
@@ -3410,7 +3410,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -3495,7 +3495,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -3580,7 +3580,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -3665,7 +3665,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -3750,7 +3750,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -3835,7 +3835,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -3920,7 +3920,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -4005,7 +4005,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -4090,7 +4090,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -4175,7 +4175,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -4260,7 +4260,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -4345,7 +4345,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -4430,7 +4430,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -4515,7 +4515,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -4600,7 +4600,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -4685,7 +4685,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -4770,7 +4770,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -4851,7 +4851,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -4932,7 +4932,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -5017,7 +5017,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -5102,7 +5102,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -5187,7 +5187,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -5272,7 +5272,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -5357,7 +5357,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -5442,7 +5442,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -5527,7 +5527,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -5588,3 +5588,6 @@ reactions
         }
     }
 }
+
+Thigh   5000;
+Tlow    200;
diff --git a/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D_GRI_TDAC/constant/chemistryProperties b/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D_GRI_TDAC/constant/chemistryProperties
index 3a230eb6f..9fd9f1912 100644
--- a/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D_GRI_TDAC/constant/chemistryProperties
+++ b/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D_GRI_TDAC/constant/chemistryProperties
@@ -24,12 +24,12 @@ chemistryType
 chemistry       on;
 
 initialChemicalTimeStep 1e-7;
-maxChemicalTimeStep     1e-4;
+//maxChemicalTimeStep 1E-3;
 
 odeCoeffs
 {
     solver          seulex;
-    absTol          1e-12;
+    absTol          1e-8;
     relTol          1e-1;
 }
 
@@ -79,7 +79,7 @@ tabulation
     printNumRetrieve   off;
 
     // Tolerance used for retrieve and grow
-    tolerance   1e-4;
+    tolerance   3e-3;
 
     // ISAT is the only method currently available
     method    ISAT;
@@ -88,7 +88,7 @@ tabulation
     scaleFactor
     {
         otherSpecies 1;
-        Temperature  2500;
+        Temperature  10000;
         Pressure     1e15;
         deltaT       1;
     }
diff --git a/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D_GRI_TDAC/constant/chemistryProperties.new b/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D_GRI_TDAC/constant/chemistryProperties.new
index 0b4586f51..330d01023 100644
--- a/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D_GRI_TDAC/constant/chemistryProperties.new
+++ b/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D_GRI_TDAC/constant/chemistryProperties.new
@@ -28,7 +28,7 @@ initialChemicalTimeStep 1e-7;
 odeCoeffs
 {
     solver          seulex;
-    absTol          1e-12;
+    absTol          1e-8;
     relTol          1e-1;
 }
 
diff --git a/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D_GRI_TDAC/constant/reactionsGRI b/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D_GRI_TDAC/constant/reactionsGRI
index 276111048..694ea19b6 100644
--- a/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D_GRI_TDAC/constant/reactionsGRI
+++ b/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D_GRI_TDAC/constant/reactionsGRI
@@ -2349,7 +2349,7 @@ reactions
         A               1.2e+11;
         beta            -1;
         Ta              0;
-        coeffs          
+        coeffs
 53
 (
 (CH4 2)
@@ -2415,7 +2415,7 @@ reactions
         A               5e+11;
         beta            -1;
         Ta              0;
-        coeffs          
+        coeffs
 53
 (
 (CH4 2)
@@ -2481,7 +2481,7 @@ reactions
         A               2.8e+12;
         beta            -0.86;
         Ta              0;
-        coeffs          
+        coeffs
 53
 (
 (CH4 1)
@@ -2547,7 +2547,7 @@ reactions
         A               1e+12;
         beta            -1;
         Ta              0;
-        coeffs          
+        coeffs
 53
 (
 (CH4 2)
@@ -2613,7 +2613,7 @@ reactions
         A               2.2e+16;
         beta            -2;
         Ta              0;
-        coeffs          
+        coeffs
 53
 (
 (CH4 2)
@@ -2679,7 +2679,7 @@ reactions
         A               1.87e+14;
         beta            -1;
         Ta              8554.25;
-        coeffs          
+        coeffs
 53
 (
 (CH4 2)
@@ -2745,7 +2745,7 @@ reactions
         A               1.06e+14;
         beta            -1.41;
         Ta              0;
-        coeffs          
+        coeffs
 53
 (
 (CH4 2)
@@ -2811,7 +2811,7 @@ reactions
         A               1.3e+11;
         beta            -0.11;
         Ta              2505.89;
-        coeffs          
+        coeffs
 53
 (
 (CH4 2)
@@ -2877,7 +2877,7 @@ reactions
         A               4.48e+13;
         beta            -1.32;
         Ta              372.362;
-        coeffs          
+        coeffs
 53
 (
 (CH4 2)
@@ -2943,7 +2943,7 @@ reactions
         A               3.1e+11;
         beta            0;
         Ta              27197.5;
-        coeffs          
+        coeffs
 53
 (
 (CH4 2)
@@ -3009,7 +3009,7 @@ reactions
         A               1.04e+26;
         beta            -3.3;
         Ta              63704;
-        coeffs          
+        coeffs
 53
 (
 (CH4 2)
@@ -3075,7 +3075,7 @@ reactions
         A               1.18e+13;
         beta            0;
         Ta              42630.4;
-        coeffs          
+        coeffs
 53
 (
 (CH4 2)
@@ -3155,7 +3155,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -3240,7 +3240,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -3325,7 +3325,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 3)
@@ -3410,7 +3410,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -3495,7 +3495,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -3580,7 +3580,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -3665,7 +3665,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -3750,7 +3750,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -3835,7 +3835,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -3920,7 +3920,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -4005,7 +4005,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -4090,7 +4090,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -4175,7 +4175,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -4260,7 +4260,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -4345,7 +4345,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -4430,7 +4430,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -4515,7 +4515,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -4600,7 +4600,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -4685,7 +4685,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -4770,7 +4770,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -4851,7 +4851,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -4932,7 +4932,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -5017,7 +5017,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -5102,7 +5102,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -5187,7 +5187,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -5272,7 +5272,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -5357,7 +5357,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -5442,7 +5442,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -5527,7 +5527,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -5588,3 +5588,6 @@ reactions
         }
     }
 }
+
+Thigh   5000;
+Tlow    200;
diff --git a/tutorials/lagrangian/reactingParcelFoam/counterFlowFlame2DLTS/constant/chemistryProperties b/tutorials/lagrangian/reactingParcelFoam/counterFlowFlame2DLTS/constant/chemistryProperties
index 3812e6131..f49edd526 100644
--- a/tutorials/lagrangian/reactingParcelFoam/counterFlowFlame2DLTS/constant/chemistryProperties
+++ b/tutorials/lagrangian/reactingParcelFoam/counterFlowFlame2DLTS/constant/chemistryProperties
@@ -33,7 +33,7 @@ EulerImplicitCoeffs
 odeCoeffs
 {
     solver          Rosenbrock43;
-    absTol          1e-12;
+    absTol          1e-8;
     relTol          0.01;
 }
 
diff --git a/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/bubbleColumnEvaporatingReacting/constant/chemistryProperties.gas b/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/bubbleColumnEvaporatingReacting/constant/chemistryProperties.gas
index d328270d6..d27415ed4 100644
--- a/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/bubbleColumnEvaporatingReacting/constant/chemistryProperties.gas
+++ b/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/bubbleColumnEvaporatingReacting/constant/chemistryProperties.gas
@@ -33,7 +33,7 @@ EulerImplicitCoeffs
 odeCoeffs
 {
     solver          Rosenbrock43;
-    absTol          1e-12;
+    absTol          1e-8;
     relTol          0.01;
 }