diff --git a/applications/solvers/cfdemSolverRhoPimple/EEqn.H b/applications/solvers/cfdemSolverRhoPimple/EEqn.H
index 022e933d..afce5d6d 100644
--- a/applications/solvers/cfdemSolverRhoPimple/EEqn.H
+++ b/applications/solvers/cfdemSolverRhoPimple/EEqn.H
@@ -20,21 +20,23 @@
 
     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;
+// For implict T terms in the energy/enthalpy transport equation, use
+// (he_n+1 - he_n) / (T_n+1 - T_n) = Cpv to eliminate T_n+1 with he_n+1.
+// This formula is valid for ideal gases with e=e(T) and h=h(T). For
+// incompressible fluids, e=e(T) holds, too, but enthalpy would need correction
+// terms accounting for pressure variations.
 
     fvScalarMatrix EEqn
     (
         fvm::ddt(rhoeps, he) + fvm::div(phi, he)
       + addSource
-     // net heat transfer from particles to fluid
       - Qsource
+      - QCoeff*T
       - fvm::Sp(QCoeff/Cpv, he)
-     // thermal conduction of the fluid with effective conductivity
+      + QCoeff/Cpv*he
+      - fvc::laplacian(voidfraction*thCond,T)
       - fvm::laplacian(voidfraction*thCond/Cpv,he)
-      // + particle-fluid energy transfer due to work
-      // + fluid energy dissipation due to shearing
+      + fvc::laplacian(voidfraction*thCond/Cpv,he)
      ==
         fvOptions(rho, he)
     );
diff --git a/applications/solvers/rcfdemSolverCoupledHeattransfer/TEqImp.H b/applications/solvers/rcfdemSolverCoupledHeattransfer/TEqImp.H
index 8c13a1e4..c79de93c 100644
--- a/applications/solvers/rcfdemSolverCoupledHeattransfer/TEqImp.H
+++ b/applications/solvers/rcfdemSolverCoupledHeattransfer/TEqImp.H
@@ -10,6 +10,7 @@
     // main contribution due to gas expansion, not due to transport of kinetic energy
     // fvc::ddt(rhoeps, K) + fvc::div(phiRec, K)
 
+    // assuming constant Cv such that e = Cv * T
     fvScalarMatrix TEqn =
     (
         fvm::ddt(rhoeps, T)
diff --git a/applications/solvers/rcfdemSolverRhoSteadyPimple/EEqn.H b/applications/solvers/rcfdemSolverRhoSteadyPimple/EEqn.H
index d9b32f48..d4011345 100644
--- a/applications/solvers/rcfdemSolverRhoSteadyPimple/EEqn.H
+++ b/applications/solvers/rcfdemSolverRhoSteadyPimple/EEqn.H
@@ -22,19 +22,33 @@
 
     Cpv = he.name() == "e" ? thermo.Cv() : thermo.Cp();
 
+// For implict T terms in the energy/enthalpy transport equation, use
+// (he_n+1 - he_n) / (T_n+1 - T_n) = Cpv to eliminate T_n+1 with he_n+1.
+// This formula is valid for ideal gases with e=e(T) and h=h(T). For
+// incompressible fluids, e=e(T) holds, too, but enthalpy would need correction
+// terms accounting for pressure variations.
 
     fvScalarMatrix EEqn
     (
         fvm::div(phi, he)
       + addSource
       - Qsource
+      - QCoeff*T
       - fvm::Sp(QCoeff/Cpv, he)
- //     - fvm::laplacian(voidfractionRec*kf/Cpv,he)
+      + QCoeff/Cpv*he
+      - fvc::laplacian(voidfractionRec*thCond,T)
       - fvm::laplacian(voidfractionRec*thCond/Cpv,he)
+      + fvc::laplacian(voidfractionRec*thCond/Cpv,he)
      ==
         fvOptions(rho, he)
     );
 
+    if (transientEEqn)
+    {
+        EEqn += fvm::ddt(rho,voidfractionRec,he);
+    }
+
+
     EEqn.relax();
 
     fvOptions.constrain(EEqn);
diff --git a/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/EEqn.H b/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/EEqn.H
index 692f4873..576fe6a4 100644
--- a/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/EEqn.H
+++ b/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/EEqn.H
@@ -22,20 +22,32 @@
 
     Cpv = he.name() == "e" ? thermo.Cv() : thermo.Cp();
 
+// For implict T terms in the energy/enthalpy transport equation, use
+// (he_n+1 - he_n) / (T_n+1 - T_n) = Cpv to eliminate T_n+1 with he_n+1.
+// This formula is valid for ideal gases with e=e(T) and h=h(T). For
+// incompressible fluids, e=e(T) holds, too, but enthalpy would need correction
+// terms accounting for pressure variations.
 
     fvScalarMatrix EEqn
     (
-        fvm::div(phi, he)
+      fvm::div(phi, he)
       + addSource
-    // net heat transfer from particles to fluid
       - Qsource
+      - QCoeff*T
       - fvm::Sp(QCoeff/Cpv, he)
- //     - fvm::laplacian(voidfractionRec*kf/Cpv,he)
+      + QCoeff/Cpv*he
+      - fvc::laplacian(voidfractionRec*thCond,T)
       - fvm::laplacian(voidfractionRec*thCond/Cpv,he)
+      + fvc::laplacian(voidfractionRec*thCond/Cpv,he)
      ==
         fvOptions(rho, he)
     );
 
+    if (transientEEqn)
+    {
+        EEqn += fvm::ddt(rho,voidfractionRec,he);
+    }
+
     EEqn.relax();
 
     fvOptions.constrain(EEqn);