diff --git a/applications/solvers/lagrangian/sprayFoam/Make/options b/applications/solvers/lagrangian/sprayFoam/Make/options
index daee4c3339..bcbb2fb5cb 100644
--- a/applications/solvers/lagrangian/sprayFoam/Make/options
+++ b/applications/solvers/lagrangian/sprayFoam/Make/options
@@ -1,4 +1,5 @@
 EXE_INC = \
+    -I../reactingParcelFoam \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I${LIB_SRC}/meshTools/lnInclude \
     -I${LIB_SRC}/sampling/lnInclude \
diff --git a/applications/solvers/lagrangian/sprayFoam/YEqn.H b/applications/solvers/lagrangian/sprayFoam/YEqn.H
deleted file mode 100644
index 6cc645bde5..0000000000
--- a/applications/solvers/lagrangian/sprayFoam/YEqn.H
+++ /dev/null
@@ -1,53 +0,0 @@
-tmp<fv::convectionScheme<scalar> > mvConvection
-(
-    fv::convectionScheme<scalar>::New
-    (
-        mesh,
-        fields,
-        phi,
-        mesh.divScheme("div(phi,Yi_hs)")
-    )
-);
-
-
-{
-    combustion->correct();
-    dQ = combustion->dQ();
-    label inertIndex = -1;
-    volScalarField Yt(0.0*Y[0]);
-
-    forAll(Y, i)
-    {
-        if (Y[i].name() != inertSpecie)
-        {
-            volScalarField& Yi = Y[i];
-
-            fvScalarMatrix YiEqn
-            (
-                fvm::ddt(rho, Yi)
-              + mvConvection->fvmDiv(phi, Yi)
-              - fvm::laplacian(turbulence->muEff(), Yi)
-              ==
-                parcels.SYi(i, Yi)
-              + combustion->R(Yi)
-              + sources(rho, Yi)
-            );
-
-            YiEqn.relax();
-
-            sources.constrain(YiEqn);
-
-            YiEqn.solve(mesh.solver("Yi"));
-
-            Yi.max(0.0);
-            Yt += Yi;
-        }
-        else
-        {
-            inertIndex = i;
-        }
-    }
-
-    Y[inertIndex] = scalar(1) - Yt;
-    Y[inertIndex].max(0.0);
-}
diff --git a/applications/solvers/lagrangian/sprayFoam/createFields.H b/applications/solvers/lagrangian/sprayFoam/createFields.H
index ead3dac75f..612ea4f3c9 100644
--- a/applications/solvers/lagrangian/sprayFoam/createFields.H
+++ b/applications/solvers/lagrangian/sprayFoam/createFields.H
@@ -6,6 +6,7 @@
     );
 
     psiReactionThermo& thermo = combustion->thermo();
+    thermo.validate(args.executable(), "h", "e");
 
     SLGThermo slgThermo(mesh, thermo);
 
@@ -23,7 +24,6 @@
     }
 
     volScalarField& p = thermo.p();
-    volScalarField& hs = thermo.he();
     const volScalarField& T = thermo.T();
     const volScalarField& psi = thermo.psi();
 
@@ -92,7 +92,7 @@
     {
         fields.add(Y[i]);
     }
-    fields.add(hs);
+    fields.add(thermo.he());
 
     volScalarField dQ
     (
diff --git a/applications/solvers/lagrangian/sprayFoam/hsEqn.H b/applications/solvers/lagrangian/sprayFoam/hsEqn.H
deleted file mode 100644
index c7c1ec0957..0000000000
--- a/applications/solvers/lagrangian/sprayFoam/hsEqn.H
+++ /dev/null
@@ -1,28 +0,0 @@
-{
-    fvScalarMatrix hsEqn
-    (
-        fvm::ddt(rho, hs)
-      + mvConvection->fvmDiv(phi, hs)
-      - fvm::laplacian(turbulence->alphaEff(), hs)
-     ==
-        dpdt
-      - (fvc::ddt(rho, K) + fvc::div(phi, K))
-      + parcels.Sh(hs)
-      + radiation->Sh(thermo)
-      + combustion->Sh()
-      + sources(rho, hs)
-    );
-
-    hsEqn.relax();
-
-    sources.constrain(hsEqn);
-
-    hsEqn.solve();
-
-    thermo.correct();
-
-    radiation->correct();
-
-    Info<< "T gas min/max   = " << min(T).value() << ", "
-        << max(T).value() << endl;
-}
diff --git a/applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/Make/options b/applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/Make/options
index 5c55d693a3..d2d83e3920 100644
--- a/applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/Make/options
+++ b/applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/Make/options
@@ -1,5 +1,6 @@
 EXE_INC = \
     -I.. \
+    -I../../reactingParcelFoam \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I${LIB_SRC}/meshTools/lnInclude \
     -I${LIB_SRC}/sampling/lnInclude \
diff --git a/applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/createClouds.H b/applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/createClouds.H
deleted file mode 100644
index ee0985ff70..0000000000
--- a/applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/createClouds.H
+++ /dev/null
@@ -1,9 +0,0 @@
-Info<< "\nConstructing reacting cloud" << endl;
-basicSprayCloud parcels
-(
-    "sprayCloud",
-    rho,
-    U,
-    g,
-    slgThermo
-);
diff --git a/applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/sprayEngineFoam.C b/applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/sprayEngineFoam.C
index baf1e40dc7..06ed6e1c2b 100644
--- a/applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/sprayEngineFoam.C
+++ b/applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/sprayEngineFoam.C
@@ -60,6 +60,7 @@ int main(int argc, char *argv[])
     #include "startSummary.H"
 
     pimpleControl pimple(mesh);
+    bool solveSpecies = true;
 
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -86,7 +87,7 @@ int main(int argc, char *argv[])
         {
             #include "UEqn.H"
             #include "YEqn.H"
-            #include "hsEqn.H"
+            #include "EEqn.H"
 
             // --- Pressure corrector loop
             while (pimple.correct())
diff --git a/applications/solvers/lagrangian/sprayFoam/sprayFoam.C b/applications/solvers/lagrangian/sprayFoam/sprayFoam.C
index fce272d042..78713032f3 100644
--- a/applications/solvers/lagrangian/sprayFoam/sprayFoam.C
+++ b/applications/solvers/lagrangian/sprayFoam/sprayFoam.C
@@ -57,6 +57,7 @@ int main(int argc, char *argv[])
     #include "setInitialDeltaT.H"
 
     pimpleControl pimple(mesh);
+    bool solveSpecies = true;
 
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -81,7 +82,7 @@ int main(int argc, char *argv[])
         {
             #include "UEqn.H"
             #include "YEqn.H"
-            #include "hsEqn.H"
+            #include "EEqn.H"
 
             // --- Pressure corrector loop
             while (pimple.correct())
diff --git a/tutorials/lagrangian/sprayFoam/aachenBomb/system/fvSchemes b/tutorials/lagrangian/sprayFoam/aachenBomb/system/fvSchemes
index c9c18972f1..03bad88200 100644
--- a/tutorials/lagrangian/sprayFoam/aachenBomb/system/fvSchemes
+++ b/tutorials/lagrangian/sprayFoam/aachenBomb/system/fvSchemes
@@ -36,7 +36,7 @@ divSchemes
     div(phi,epsilon) Gauss upwind;
     div(U)          Gauss linear;
     div((muEff*dev2(T(grad(U))))) Gauss linear;
-    div(phi,Yi_hs)  Gauss upwind;
+    div(phi,Yi_h)   Gauss upwind;
 }
 
 laplacianSchemes