diff --git a/applications/solvers/cfdemSolverRhoPimple/EEqn.H b/applications/solvers/cfdemSolverRhoPimple/EEqn.H
index b2a30416..08af6a05 100644
--- a/applications/solvers/cfdemSolverRhoPimple/EEqn.H
+++ b/applications/solvers/cfdemSolverRhoPimple/EEqn.H
@@ -23,6 +23,9 @@
 
     Cpv = he.name() == "e" ? thermo.Cv() : thermo.Cp();
 
+    // correct source for the thermodynamic reference temperature
+    dimensionedScalar Tref("Tref", dimTemperature, T[0]-he[0]/(Cpv[0]+SMALL));
+    Qsource += QCoeff*Tref;
 
     fvScalarMatrix EEqn
     (
diff --git a/applications/solvers/cfdemSolverRhoPimpleChem/EEqn.H b/applications/solvers/cfdemSolverRhoPimpleChem/EEqn.H
index edb4081f..5dfd1a5a 100644
--- a/applications/solvers/cfdemSolverRhoPimpleChem/EEqn.H
+++ b/applications/solvers/cfdemSolverRhoPimpleChem/EEqn.H
@@ -9,6 +9,10 @@ particleCloud.energyCoefficients(QCoeff);
 thCond=particleCloud.thermCondM().thermCond();
 Cpv = he.name() == "e" ? thermo.Cv() : thermo.Cp();
 
+// correct source for the thermodynamic reference temperature
+dimensionedScalar Tref("Tref", dimTemperature, T[0]-he[0]/(Cpv[0]+SMALL));
+Qsource += QCoeff*Tref;
+
 fvScalarMatrix EEqn
 (
     fvm::ddt(rhoeps, he) + fvm::div(phi, he)
diff --git a/applications/solvers/cfdemSolverRhoSimple/EEqn.H b/applications/solvers/cfdemSolverRhoSimple/EEqn.H
index 3dbb5c3f..53630bfb 100644
--- a/applications/solvers/cfdemSolverRhoSimple/EEqn.H
+++ b/applications/solvers/cfdemSolverRhoSimple/EEqn.H
@@ -25,6 +25,9 @@
 
     Cpv = he.name() == "e" ? thermo.Cv() : thermo.Cp();
 
+    // correct source for the thermodynamic reference temperature
+    dimensionedScalar Tref("Tref", dimTemperature, T[0]-he[0]/(Cpv[0]+SMALL));
+    Qsource += QCoeff*Tref;
 
     fvScalarMatrix EEqn
     (