diff --git a/applications/solvers/lagrangian/trackedReactingParcelFoam/chemistry.H b/applications/solvers/lagrangian/trackedReactingParcelFoam/chemistry.H
new file mode 100644
index 0000000000000000000000000000000000000000..07b1e9953b0db867186f6c668d27a9415a26c265
--- /dev/null
+++ b/applications/solvers/lagrangian/trackedReactingParcelFoam/chemistry.H
@@ -0,0 +1,25 @@
+{
+    Info << "Solving chemistry" << endl;
+
+    chemistry.solve
+    (
+        runTime.value() - runTime.deltaT().value(),
+        runTime.deltaT().value()
+    );
+
+    // turbulent time scale
+    if (turbulentReaction)
+    {
+        DimensionedField<scalar, volMesh> tk =
+            Cmix*sqrt(turbulence->muEff()/rho/turbulence->epsilon());
+        DimensionedField<scalar, volMesh> tc =
+            chemistry.tc()().dimensionedInternalField();
+
+        // Chalmers PaSR model
+        kappa = (runTime.deltaT() + tc)/(runTime.deltaT() + tc + tk);
+    }
+    else
+    {
+        kappa = 1.0;
+    }
+}
diff --git a/applications/solvers/lagrangian/trackedReactingParcelFoam/hEqn.H b/applications/solvers/lagrangian/trackedReactingParcelFoam/hEqn.H
new file mode 100644
index 0000000000000000000000000000000000000000..d686e452f46f053c7b52f7b19872f964d9c054ef
--- /dev/null
+++ b/applications/solvers/lagrangian/trackedReactingParcelFoam/hEqn.H
@@ -0,0 +1,20 @@
+{
+    fvScalarMatrix hEqn
+    (
+        fvm::ddt(rho, h)
+      + mvConvection->fvmDiv(phi, h)
+      - fvm::laplacian(turbulence->alphaEff(), h)
+     ==
+        DpDt
+     +  parcels.Sh()
+     +  radiation->Sh(thermo)
+    );
+
+    hEqn.relax();
+
+    hEqn.solve();
+
+    thermo.correct();
+
+    radiation->correct();
+}
diff --git a/applications/solvers/lagrangian/trackedReactingParcelFoam/readChemistryProperties.H b/applications/solvers/lagrangian/trackedReactingParcelFoam/readChemistryProperties.H
new file mode 100644
index 0000000000000000000000000000000000000000..1a60e6fb34645a004fd39321f7a54d3bd5b45381
--- /dev/null
+++ b/applications/solvers/lagrangian/trackedReactingParcelFoam/readChemistryProperties.H
@@ -0,0 +1,22 @@
+Info<< "Reading chemistry properties\n" << endl;
+
+IOdictionary chemistryProperties
+(
+    IOobject
+    (
+        "chemistryProperties",
+        runTime.constant(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::NO_WRITE
+    )
+);
+
+Switch turbulentReaction(chemistryProperties.lookup("turbulentReaction"));
+
+dimensionedScalar Cmix("Cmix", dimless, 1.0);
+
+if (turbulentReaction)
+{
+    chemistryProperties.lookup("Cmix") >> Cmix;
+}