diff --git a/dumux/freeflow/staggered/fluxvariables.hh b/dumux/freeflow/staggered/fluxvariables.hh
index bb9c635422a61e30aa16c1adb207b570df9b82ff..a344047526f322bddf33947bf4e79e5f798aebf2 100644
--- a/dumux/freeflow/staggered/fluxvariables.hh
+++ b/dumux/freeflow/staggered/fluxvariables.hh
@@ -38,7 +38,7 @@ NEW_PROP_TAG(EnableInertiaTerms);
 }
 
 // forward declaration
-template<class TypeTag, bool enableComponentTransport, bool enableEnergyBalance>
+template<class TypeTag, bool enableComponentTransport>
 class FreeFlowFluxVariablesImpl;
 
 /*!
@@ -48,8 +48,7 @@ class FreeFlowFluxVariablesImpl;
  * \note  Not all specializations are currently implemented
  */
 template<class TypeTag>
-using FreeFlowFluxVariables = FreeFlowFluxVariablesImpl<TypeTag, GET_PROP_VALUE(TypeTag, EnableComponentTransport),
-                                                                 GET_PROP_VALUE(TypeTag, EnableEnergyBalanceStokes)>;
+using FreeFlowFluxVariables = FreeFlowFluxVariablesImpl<TypeTag, GET_PROP_VALUE(TypeTag, EnableComponentTransport)>;
 
 /*!
  * \ingroup Discretization
@@ -58,7 +57,7 @@ using FreeFlowFluxVariables = FreeFlowFluxVariablesImpl<TypeTag, GET_PROP_VALUE(
  */
 // specialization for immiscible, isothermal flow
 template<class TypeTag>
-class FreeFlowFluxVariablesImpl<TypeTag, false, false>
+class FreeFlowFluxVariablesImpl<TypeTag, false>
 : public FluxVariablesBase<TypeTag>
 {
     using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
diff --git a/dumux/freeflow/staggered/localresidual.hh b/dumux/freeflow/staggered/localresidual.hh
index 3b18a7bfaae6d754cd1cb35392e7f6c64bfc24df..54a5030441cde0c9044591ec7aa4e81a361aab58 100644
--- a/dumux/freeflow/staggered/localresidual.hh
+++ b/dumux/freeflow/staggered/localresidual.hh
@@ -90,6 +90,7 @@ class StaggeredNavierStokesResidualImpl<TypeTag, false> : public Dumux::Staggere
     using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
     using FluxVariables = typename GET_PROP_TYPE(TypeTag, FluxVariables);
     using EnergyLocalResidual = typename GET_PROP_TYPE(TypeTag, EnergyLocalResidual);
+    using EnergyFluxVariables = typename GET_PROP_TYPE(TypeTag, EnergyFluxVariables);
 
 
     using DofTypeIndices = typename GET_PROP(TypeTag, DofTypeIndices);
@@ -129,8 +130,12 @@ public:
                                                         const ElementFluxVariablesCache& elemFluxVarsCache)
     {
         FluxVariables fluxVars;
-        return fluxVars.computeFluxForCellCenter(this->problem(), element, fvGeometry, elemVolVars,
+        CellCenterPrimaryVariables flux = fluxVars.computeFluxForCellCenter(this->problem(), element, fvGeometry, elemVolVars,
                                                  globalFaceVars, scvf, elemFluxVarsCache[scvf]);
+
+        EnergyFluxVariables::energyFlux(flux, this->problem(), element, fvGeometry, elemVolVars, globalFaceVars, scvf, elemFluxVarsCache[scvf]);
+
+        return flux;
     }
 
     CellCenterPrimaryVariables computeSourceForCellCenter(const Element &element,
diff --git a/dumux/freeflow/staggeredni/fluxvariables.hh b/dumux/freeflow/staggeredni/fluxvariables.hh
index 661f02221ec3df51986a4218fe451c20f79dc4c3..346fe740a542fff78e7488929014dcf6905fcc91 100644
--- a/dumux/freeflow/staggeredni/fluxvariables.hh
+++ b/dumux/freeflow/staggeredni/fluxvariables.hh
@@ -24,8 +24,6 @@
 #define DUMUX_FREELOW_IMPLICIT_NI_FLUXVARIABLES_HH
 
 #include <dumux/implicit/properties.hh>
-#include <dumux/discretization/fluxvariablesbase.hh>
-#include "../staggered/fluxvariables.hh"
 
 namespace Dumux
 {
@@ -34,7 +32,6 @@ namespace Properties
 {
 // forward declaration
 NEW_PROP_TAG(EnableEnergyBalanceStokes);
-NEW_PROP_TAG(EnableInertiaTerms);
 }
 
 /*!
@@ -44,9 +41,45 @@ NEW_PROP_TAG(EnableInertiaTerms);
  * \note  Not all specializations are currently implemented
  */
 
-// specialization for immiscible, non-isothermal flow
+// forward declaration
+template<class TypeTag, bool enableEnergyBalance>
+class FreeFlowEnergyFluxVariablesImplementation;
+
 template<class TypeTag>
-class FreeFlowFluxVariablesImpl<TypeTag, false, true> : public FreeFlowFluxVariablesImpl<TypeTag, false, false>
+using FreeFlowEnergyFluxVariables = FreeFlowEnergyFluxVariablesImplementation<TypeTag, GET_PROP_VALUE(TypeTag, EnableEnergyBalanceStokes)>;
+
+// specialization for isothermal flow
+template<class TypeTag>
+class FreeFlowEnergyFluxVariablesImplementation<TypeTag, false>
+{
+    using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+    using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+    using Element = typename GridView::template Codim<0>::Entity;
+    using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+    using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
+    using GlobalFaceVars = typename GET_PROP_TYPE(TypeTag, GlobalFaceVars);
+    using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+    using FluxVariablesCache = typename GET_PROP_TYPE(TypeTag, FluxVariablesCache);
+    using CellCenterPrimaryVariables = typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables);
+
+public:
+
+    static void energyFlux(CellCenterPrimaryVariables& flux,
+                           const Problem& problem,
+                           const Element &element,
+                           const FVElementGeometry& fvGeometry,
+                           const ElementVolumeVariables& elemVolVars,
+                           const GlobalFaceVars& globalFaceVars,
+                           const SubControlVolumeFace &scvf,
+                           const FluxVariablesCache& fluxVarsCache)
+    { }
+
+};
+
+// specialization for non-isothermal flow
+template<class TypeTag>
+class FreeFlowEnergyFluxVariablesImplementation<TypeTag, true>
 {
     using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
     using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
@@ -63,40 +96,21 @@ class FreeFlowFluxVariablesImpl<TypeTag, false, true> : public FreeFlowFluxVaria
     using HeatConductionType = typename GET_PROP_TYPE(TypeTag, HeatConductionType);
     static constexpr bool useMoles = GET_PROP_VALUE(TypeTag, UseMoles);
 
-    using ParentType = FreeFlowFluxVariablesImpl<TypeTag, false, false>;
-
-    enum {
-         // grid and world dimension
-        dim = GridView::dimension,
-        dimWorld = GridView::dimensionworld,
-
-        energyBalanceIdx = Indices::energyBalanceIdx
-    };
+    enum { energyBalanceIdx = Indices::energyBalanceIdx };
 
 public:
-    /*!
-    * \brief Returns the fluxes for the cell center primary variables
-    * \param problem The problem
-    * \param element The element
-    * \param fvGeometry The finite-volume geometry
-    * \param elemVolVars All volume variables for the element
-    * \param globalFaceVars The face variables
-    * \param scvf The sub control volume face
-    * \param fluxVarsCache The flux variables cache
-    */
-    CellCenterPrimaryVariables computeFluxForCellCenter(const Problem& problem,
-                                                        const Element &element,
-                                                        const FVElementGeometry& fvGeometry,
-                                                        const ElementVolumeVariables& elemVolVars,
-                                                        const GlobalFaceVars& globalFaceVars,
-                                                        const SubControlVolumeFace &scvf,
-                                                        const FluxVariablesCache& fluxVarsCache)
+
+    static void energyFlux(CellCenterPrimaryVariables& flux,
+                           const Problem& problem,
+                           const Element &element,
+                           const FVElementGeometry& fvGeometry,
+                           const ElementVolumeVariables& elemVolVars,
+                           const GlobalFaceVars& globalFaceVars,
+                           const SubControlVolumeFace &scvf,
+                           const FluxVariablesCache& fluxVarsCache)
     {
-        CellCenterPrimaryVariables flux(0.0);
-        flux = ParentType::computeFluxForCellCenter(problem, element, fvGeometry, elemVolVars, globalFaceVars, scvf, fluxVarsCache);
-        flux += advectiveFluxForCellCenter_(problem, fvGeometry, elemVolVars, globalFaceVars, scvf);
-        flux += HeatConductionType::diffusiveFluxForCellCenter(problem, element, fvGeometry, elemVolVars, scvf);
-        return flux;
+        flux[energyBalanceIdx] += advectiveFluxForCellCenter_(problem, fvGeometry, elemVolVars, globalFaceVars, scvf);
+        flux[energyBalanceIdx] += HeatConductionType::diffusiveFluxForCellCenter(problem, element, fvGeometry, elemVolVars, scvf);
     }
 
 private:
@@ -108,13 +122,13 @@ private:
     * \param globalFaceVars The face variables
     * \param scvf The sub control volume face
     */
-    CellCenterPrimaryVariables advectiveFluxForCellCenter_(const Problem& problem,
-                                                          const FVElementGeometry& fvGeometry,
-                                                          const ElementVolumeVariables& elemVolVars,
-                                                          const GlobalFaceVars& globalFaceVars,
-                                                          const SubControlVolumeFace &scvf)
+    static Scalar advectiveFluxForCellCenter_(const Problem& problem,
+                                              const FVElementGeometry& fvGeometry,
+                                              const ElementVolumeVariables& elemVolVars,
+                                              const GlobalFaceVars& globalFaceVars,
+                                              const SubControlVolumeFace &scvf)
     {
-        CellCenterPrimaryVariables flux(0.0);
+        Scalar flux(0.0);
 
         const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx());
         const auto& insideVolVars = elemVolVars[insideScv];
@@ -143,7 +157,7 @@ private:
         const Scalar upstreamEnthalpy = upstreamVolVars.enthalpy();
         const Scalar downstreamEnthalpy = downstreamVolVars.enthalpy();
 
-        flux[energyBalanceIdx] = (upWindWeight * upstreamDensity * upstreamEnthalpy
+        flux = (upWindWeight * upstreamDensity * upstreamEnthalpy
                                  + (1.0 - upWindWeight) * downstreamDensity * downstreamEnthalpy)
                                  * velocity;
         flux *= scvf.area() * sign(scvf.outerNormalScalar());