diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.C b/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.C
index 56407ef4..422cd839 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.C
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.C
@@ -51,7 +51,8 @@ Foam::chemistryModel<CompType, ThermoType>::chemistryModel
     nSpecie_(Y_.size()),
     nReaction_(reactions_.size()),
 
-    RR_(nSpecie_)
+    RR_(nSpecie_),
+    eRR_(nSpecie_)
 {
     // create the fields for the chemistry sources
     forAll(RR_, fieldI)
@@ -73,10 +74,45 @@ Foam::chemistryModel<CompType, ThermoType>::chemistryModel
                 dimensionedScalar("zero", dimMass/dimVolume/dimTime, 0.0)
             )
         );
+
+        eRR_.set
+        (
+            fieldI,
+            new DimensionedField<scalar, volMesh>
+            (
+                IOobject
+                (
+                    "eRR." + Y_[fieldI].name(),
+                    mesh.time().timeName(),
+                    mesh,
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE
+                ),
+                mesh,
+                dimensionedScalar("zero", dimMass/dimVolume/dimTime, 0.0)
+            )
+        );
     }
 
     Info<< "chemistryModel: Number of species = " << nSpecie_
         << " and reactions = " << nReaction_ << endl;
+
+    label neReactions_ = 0;
+
+    for (label i=0; i<nReaction_; i++)
+    {
+        // if (reactions_[i].species().contains("E-"))
+        if (reactions_[i].hasElectron())
+        {
+            neReactions_++;
+        }
+        else
+        {
+        }
+    }
+
+    Info<< "chemistryModel: Number electron reactions = " << neReactions_ << endl;
+
 }
 
 
@@ -132,6 +168,53 @@ Foam::chemistryModel<CompType, ThermoType>::omega
 }
 
 
+template<class CompType, class ThermoType>
+Foam::tmp<Foam::scalarField>
+Foam::chemistryModel<CompType, ThermoType>::eomega
+(
+    const scalarField& c,
+    const scalar T,
+    const scalar p
+) const
+{
+    scalar pf, cf, pr, cr;
+    label lRef, rRef;
+
+    tmp<scalarField> tom(new scalarField(nEqns(), 0.0));
+    scalarField& om = tom();
+
+    forAll(reactions_, i)
+    {
+        const Reaction<ThermoType>& R = reactions_[i];
+
+        // if (R.species().contains("E-"))
+        if (R.hasElectron())
+        {
+            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;
+                om[si] -= sl*omegai;
+            }
+
+            forAll(R.rhs(), s)
+            {
+                const label si = R.rhs()[s].index;
+                const scalar sr = R.rhs()[s].stoichCoeff;
+                om[si] += sr*omegai;
+            }
+        }
+    }
+
+    return tom;
+}
+
+
 template<class CompType, class ThermoType>
 Foam::scalar Foam::chemistryModel<CompType, ThermoType>::omegaI
 (
@@ -576,7 +659,8 @@ Foam::chemistryModel<CompType, ThermoType>::Sh() const
             forAll(Sh, cellI)
             {
                 const scalar hi = specieThermo_[i].Hc();
-                Sh[cellI] -= hi*RR_[i][cellI];
+                Sh[cellI] -= hi*(RR_[i][cellI]-eRR_[i][cellI]);
+                // Sh[cellI] -= hi*(RR_[i][cellI]);
             }
         }
     }
@@ -769,7 +853,8 @@ void Foam::chemistryModel<CompType, ThermoType>::calculate()
 
     for (label i=0; i<nReaction_; i++)
     {
-        if (reactions_[i].species().contains("E-"))
+        // if (reactions_[i].species().contains("E-"))
+        if (reactions_[i].hasElectron())
         {
             eRxnLabels.append(i);
             neReactions_++;
@@ -810,16 +895,19 @@ void Foam::chemistryModel<CompType, ThermoType>::calculate()
         }
 
         const scalarField dcdt(omega(c, Ti, pi));
+        const scalarField dcdte(eomega(c, Ti, pi));
 
         for (label i=0; i<nSpecie_; i++)
         {
             if (Y_[i].name() == "E-")
             {
                 RR_[i][celli] = dcdt[i]*thermoType::NA;
+                eRR_[i][celli] = dcdt[i]*thermoType::NA;
             }
             else
             {
                 RR_[i][celli] = dcdt[i]*specieThermo_[i].W();
+                eRR_[i][celli] = dcdte[i]*specieThermo_[i].W();
             }
         }
     }
@@ -910,7 +998,8 @@ Foam::scalar Foam::chemistryModel<CompType, ThermoType>::solve
 
     for (label i=0; i<nReaction_; i++)
     {
-        if (reactions_[i].species().contains("E-"))
+        // if (reactions_[i].species().contains("E-"))
+        if (reactions_[i].hasElectron())
         {
             eRxnLabels.append(i);
             neReactions_++;
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.H b/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.H
index 368090a2..c3ede47a 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.H
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.H
@@ -100,6 +100,9 @@ protected:
         //- List of reaction rate per specie [kg/m3/s]
         PtrList<DimensionedField<scalar, volMesh> > RR_;
 
+        //- List of source rate per specie by electron reactions [kg/m3/s]
+        PtrList<DimensionedField<scalar, volMesh> > eRR_;
+
 
     // Protected Member Functions
 
@@ -107,6 +110,10 @@ protected:
         //  (e.g. for multi-chemistry model)
         inline PtrList<DimensionedField<scalar, volMesh> >& RR();
 
+        //- Write access to chemical source terms
+        //  (e.g. for multi-chemistry model)
+        inline PtrList<DimensionedField<scalar, volMesh> >& eRR();
+
 
 public:
 
@@ -146,6 +153,14 @@ public:
             const scalar p
         ) const;
 
+        //- rate of change in concentration, for each species, by electron reactions
+        tmp<scalarField> eomega
+        (
+            const scalarField& c,
+            const scalar T,
+            const scalar p
+        ) const;
+
         //- Return the reaction rate for reaction r and the reference
         //  species and charateristic times
         virtual scalar omega
@@ -198,6 +213,18 @@ public:
                 const label i
             );
 
+            //- Return const access to the chemical source terms for specie, i
+            inline const DimensionedField<scalar, volMesh>& eRR
+            (
+                const label i
+            ) const;
+
+            //- Return non const access to chemical source terms [kg/m3/s]
+            virtual DimensionedField<scalar, volMesh>& eRR
+            (
+                const label i
+            );
+
             //- Return reaction rate of the specieI in reactionI
             virtual tmp<DimensionedField<scalar, volMesh> > calculateRR
             (
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModelI.H b/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModelI.H
index c009a11e..ec8791a2 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModelI.H
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModelI.H
@@ -89,4 +89,33 @@ Foam::chemistryModel<CompType, ThermoType>::RR
 }
 
 
+template<class CompType, class ThermoType>
+inline Foam::PtrList<Foam::DimensionedField<Foam::scalar, Foam::volMesh> >&
+Foam::chemistryModel<CompType, ThermoType>::eRR()
+{
+    return eRR_;
+}
+
+
+template<class CompType, class ThermoType>
+inline const Foam::DimensionedField<Foam::scalar, Foam::volMesh>&
+Foam::chemistryModel<CompType, ThermoType>::eRR
+(
+    const label i
+) const
+{
+    return eRR_[i];
+}
+
+template<class CompType, class ThermoType>
+Foam::DimensionedField<Foam::scalar, Foam::volMesh>&
+Foam::chemistryModel<CompType, ThermoType>::eRR
+(
+    const label i
+)
+{
+    return eRR_[i];
+}
+
+
 // ************************************************************************* //
diff --git a/src/thermophysicalModels/specie/specie/specie.C b/src/thermophysicalModels/specie/specie/specie.C
index 52799c4b..95a7ea96 100644
--- a/src/thermophysicalModels/specie/specie/specie.C
+++ b/src/thermophysicalModels/specie/specie/specie.C
@@ -44,6 +44,9 @@ const Foam::scalar Foam::specie::e = constant::electromagnetic::e.value();
 //- Avogadro number (default in [1/mol])
 const Foam::scalar Foam::specie::NA = constant::physicoChemical::NA.value()*1000;
 
+//- Boltzmann constant (default in [J/K])
+const Foam::scalar Foam::specie::k = constant::physicoChemical::k.value();
+
 namespace Foam
 {
 defineTypeNameAndDebug(specie, 0);
diff --git a/src/thermophysicalModels/specie/specie/specie.H b/src/thermophysicalModels/specie/specie/specie.H
index 3b1a4838..8e555f3d 100644
--- a/src/thermophysicalModels/specie/specie/specie.H
+++ b/src/thermophysicalModels/specie/specie/specie.H
@@ -94,6 +94,9 @@ public:
             //- Avogadro number [1/kmol]
             static const scalar NA;
 
+            //- Boltzmann constant [J/K]
+            static const scalar k;
+
 
     // Constructors