diff --git a/dumux/discretization/staggered/freeflow/fourierslaw.hh b/dumux/discretization/staggered/freeflow/fourierslaw.hh
index f27080c9836465d1336da04845f1d4d1796732f2..d83273a35f76b4c8f463b1ffd282a7807d0e0acf 100644
--- a/dumux/discretization/staggered/freeflow/fourierslaw.hh
+++ b/dumux/discretization/staggered/freeflow/fourierslaw.hh
@@ -44,7 +44,6 @@ template <class TypeTag>
 class FouriersLawImplementation<TypeTag, DiscretizationMethod::staggered >
 {
     using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
-    using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
     using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
     using SubControlVolume = typename FVElementGeometry::SubControlVolume;
     using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
@@ -63,22 +62,14 @@ public:
     //! We don't cache anything for this law
     using Cache = FluxVariablesCaching::EmptyDiffusionCache;
 
-    //! calculate the molecular diffusive fluxes
-    static Scalar diffusiveFluxForCellCenter(const Problem& problem,
-                                             const Element& element,
-                                             const FVElementGeometry& fvGeometry,
-                                             const ElementVolumeVariables& elemVolVars,
-                                             const SubControlVolumeFace &scvf)
+    //! calculate the diffusive energy fluxes
+    static Scalar flux(const Element& element,
+                       const FVElementGeometry& fvGeometry,
+                       const ElementVolumeVariables& elemVolVars,
+                       const SubControlVolumeFace &scvf)
     {
         Scalar flux(0.0);
 
-        if(scvf.boundary())
-        {
-            const auto bcTypes = problem.boundaryTypesAtPos(scvf.center());
-            if(bcTypes.isOutflow(energyBalanceIdx) || bcTypes.isNeumann(energyBalanceIdx))
-                return flux;
-        }
-
         const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx());
         const auto& outsideScv = fvGeometry.scv(scvf.outsideScvIdx());
         const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
diff --git a/dumux/freeflow/compositional/staggered/localresidual.hh b/dumux/freeflow/compositional/staggered/localresidual.hh
index 099589dcfcc8bce850d80eb8ed6a3d117d65910b..2e9525c60541ef28cc80ce1f3e25ca7e9d63883a 100644
--- a/dumux/freeflow/compositional/staggered/localresidual.hh
+++ b/dumux/freeflow/compositional/staggered/localresidual.hh
@@ -59,6 +59,8 @@ class FreeflowNCResidualImpl<TypeTag, BaseLocalResidual, DiscretizationMethod::s
     static constexpr bool useMoles = GET_PROP_VALUE(TypeTag, UseMoles);
     static constexpr auto cellCenterOffset = ParentType::cellCenterOffset;
 
+    using EnergyLocalResidual = typename ParentType::EnergyLocalResidual;
+
 public:
     using ParentType::ParentType;
 
@@ -88,8 +90,7 @@ public:
         if(Indices::replaceCompEqIdx < numComponents)
             storage[Indices::replaceCompEqIdx] = density;
 
-        this->computeStorageForCellCenterNonIsothermal_(std::integral_constant<bool, ModelTraits::enableEnergyBalance() >(),
-                                                        problem, scv, volVars, storage);
+        EnergyLocalResidual::fluidPhaseStorage(storage, volVars);
 
         return storage;
     }
diff --git a/dumux/freeflow/compositional/volumevariables.hh b/dumux/freeflow/compositional/volumevariables.hh
index 92576fc064a3e61533865a2c481949ee360139e2..be8632b0d0947c0263a7ccfda52290ad0e66a857 100644
--- a/dumux/freeflow/compositional/volumevariables.hh
+++ b/dumux/freeflow/compositional/volumevariables.hh
@@ -42,9 +42,8 @@ class FreeflowNCVolumeVariables : public FreeFlowVolumeVariables< Traits, Freefl
     using ParentType = FreeFlowVolumeVariables<Traits, ThisType>;
 
     using Scalar = typename Traits::PrimaryVariables::value_type;
-    using Indices = typename Traits::ModelTraits::Indices;
 
-    static constexpr int fluidSystemPhaseIdx = Indices::fluidSystemPhaseIdx;
+    static constexpr int fluidSystemPhaseIdx = Traits::ModelTraits::Indices::fluidSystemPhaseIdx;
     static constexpr int numComponents = Traits::ModelTraits::numComponents();
 
     static constexpr bool useMoles = Traits::ModelTraits::useMoles();
@@ -54,6 +53,8 @@ public:
     using FluidSystem = typename Traits::FluidSystem;
     //! export the fluid state type
     using FluidState = typename Traits::FluidState;
+    //! export the indices type
+    using Indices = typename Traits::ModelTraits::Indices;
 
     /*!
      * \brief Update all quantities for a given control volume
diff --git a/dumux/freeflow/navierstokes/model.hh b/dumux/freeflow/navierstokes/model.hh
index ad0c3983ae4f78dc3f5f97ac28388ce7c10b0105..3d4ef755033d4ed1764c1829c7384491e7fb834a 100644
--- a/dumux/freeflow/navierstokes/model.hh
+++ b/dumux/freeflow/navierstokes/model.hh
@@ -232,9 +232,6 @@ public:
      using type = FreeflowNonIsothermalVtkOutputFields<IsothermalFields, ModelTraits>;
 };
 
-//! Use Fourier's Law as default heat conduction type
-SET_TYPE_PROP(NavierStokesNI, HeatConductionType, FouriersLaw<TypeTag>);
-
  // \}
 }
 
diff --git a/dumux/freeflow/navierstokes/staggered/fluxvariables.hh b/dumux/freeflow/navierstokes/staggered/fluxvariables.hh
index f3191ef751c020b39e9ea5a89c9754a01f0ec06e..afd4f100ee65dea8a41aee464a2e5700cc8e35f6 100644
--- a/dumux/freeflow/navierstokes/staggered/fluxvariables.hh
+++ b/dumux/freeflow/navierstokes/staggered/fluxvariables.hh
@@ -84,6 +84,8 @@ class NavierStokesFluxVariablesImpl<TypeTag, DiscretizationMethod::staggered>
 
 public:
 
+    using HeatConductionType = typename GET_PROP_TYPE(TypeTag, HeatConductionType);
+
     /*!
     * \brief Returns the advective flux over a sub control volume face.
     * \param elemVolVars All volume variables for the element
diff --git a/dumux/freeflow/navierstokes/staggered/localresidual.hh b/dumux/freeflow/navierstokes/staggered/localresidual.hh
index d96aa7a3bf4b35e539f28ab4fae4b85434e23b3e..3f66da3e1ac5a3ba76688f6345c6483172065513 100644
--- a/dumux/freeflow/navierstokes/staggered/localresidual.hh
+++ b/dumux/freeflow/navierstokes/staggered/localresidual.hh
@@ -28,6 +28,7 @@
 #include <dumux/discretization/methods.hh>
 #include <dumux/assembly/staggeredlocalresidual.hh>
 #include <dune/common/hybridutilities.hh>
+#include <dumux/freeflow/nonisothermal/localresidual.hh>
 
 namespace Dumux
 {
@@ -79,6 +80,7 @@ class NavierStokesResidualImpl<TypeTag, DiscretizationMethod::staggered>
     using ModelTraits = typename GET_PROP_TYPE(TypeTag, ModelTraits);
 
 public:
+    using EnergyLocalResidual = FreeFlowEnergyLocalResidual<FVGridGeometry, FluxVariables, ModelTraits::enableEnergyBalance()>;
 
     // account for the offset of the cell center privars within the PrimaryVariables container
     static constexpr auto cellCenterOffset = ModelTraits::numEq() - CellCenterPrimaryVariables::dimension;
@@ -100,8 +102,7 @@ public:
         CellCenterPrimaryVariables flux = fluxVars.computeFluxForCellCenter(problem, element, fvGeometry, elemVolVars,
                                                  elemFaceVars, scvf, elemFluxVarsCache[scvf]);
 
-        computeFluxForCellCenterNonIsothermal_(std::integral_constant<bool, ModelTraits::enableEnergyBalance()>(),
-                                               problem, element, fvGeometry, elemVolVars, elemFaceVars, scvf, elemFluxVarsCache, flux);
+        EnergyLocalResidual::heatFlux(flux, problem, element, fvGeometry, elemVolVars, elemFaceVars, scvf);
 
         return flux;
     }
@@ -135,8 +136,7 @@ public:
         CellCenterPrimaryVariables storage;
         storage[Indices::conti0EqIdx - ModelTraits::dim()] = volVars.density();
 
-        computeStorageForCellCenterNonIsothermal_(std::integral_constant<bool, ModelTraits::enableEnergyBalance() >(),
-                                                  problem, scv, volVars, storage);
+        EnergyLocalResidual::fluidPhaseStorage(storage, volVars);
 
         return storage;
     }
@@ -204,24 +204,6 @@ public:
 
 protected:
 
-    //  /*!
-    //  * \brief Evaluate boundary conditions
-    //  */
-    // template<class ElementBoundaryTypes>
-    // void evalBoundary_(const Element& element,
-    //                    const FVElementGeometry& fvGeometry,
-    //                    const ElementVolumeVariables& elemVolVars,
-    //                    const ElementFaceVariables& elemFaceVars,
-    //                    const ElementBoundaryTypes& elemBcTypes,
-    //                    const ElementFluxVariablesCache& elemFluxVarsCache)
-    // {
-    //     evalBoundaryForCellCenter_(element, fvGeometry, elemVolVars, elemFaceVars, elemBcTypes, elemFluxVarsCache);
-    //     for (auto&& scvf : scvfs(fvGeometry))
-    //     {
-    //         evalBoundaryForFace_(element, fvGeometry, scvf, elemVolVars, elemFaceVars, elemBcTypes, elemFluxVarsCache);
-    //     }
-    // }
-
      /*!
      * \brief Evaluate boundary conditions for a cell center dof
      */
@@ -245,6 +227,7 @@ protected:
                 // handle the actual boundary conditions:
                 const auto bcTypes = problem.boundaryTypes(element, scvf);
 
+                EnergyLocalResidual::heatFlux(boundaryFlux, problem, element, fvGeometry, elemVolVars, elemFaceVars, scvf);
                 if(bcTypes.hasNeumann())
                 {
                     static constexpr auto numEqCellCenter = CellCenterResidual::dimension;
@@ -327,52 +310,6 @@ protected:
         }
     }
 
-    //! Evaluate energy fluxes entering or leaving the cell center control volume for non isothermal models
-    void computeFluxForCellCenterNonIsothermal_(std::true_type,
-                                                const Problem& problem,
-                                                const Element &element,
-                                                const FVElementGeometry& fvGeometry,
-                                                const ElementVolumeVariables& elemVolVars,
-                                                const ElementFaceVariables& elemFaceVars,
-                                                const SubControlVolumeFace &scvf,
-                                                const ElementFluxVariablesCache& elemFluxVarsCache,
-                                                CellCenterPrimaryVariables& flux) const
-    {
-        // if we are on an inflow/outflow boundary, use the volVars of the element itself
-        // TODO: catch neumann and outflow in localResidual's evalBoundary_()
-        bool isOutflow = false;
-        if(scvf.boundary())
-        {
-            const auto bcTypes = problem.boundaryTypesAtPos(scvf.center());
-                if(bcTypes.isOutflow(Indices::energyBalanceIdx))
-                    isOutflow = true;
-        }
-
-        auto upwindTerm = [](const auto& volVars) { return volVars.density() * volVars.enthalpy(); };
-        using HeatConductionType = typename GET_PROP_TYPE(TypeTag, HeatConductionType);
-
-        flux[Indices::energyBalanceIdx - cellCenterOffset] = FluxVariables::advectiveFluxForCellCenter(elemVolVars, elemFaceVars, scvf, upwindTerm, isOutflow);
-        flux[Indices::energyBalanceIdx - cellCenterOffset] += HeatConductionType::diffusiveFluxForCellCenter(problem, element, fvGeometry, elemVolVars, scvf);
-    }
-
-    //! Evaluate energy fluxes entering or leaving the cell center control volume for non isothermal models
-    template <typename... Args>
-    void computeFluxForCellCenterNonIsothermal_(std::false_type, Args&&... args) const {}
-
-    //! Evaluate energy storage for non isothermal models
-    void computeStorageForCellCenterNonIsothermal_(std::true_type,
-                                                   const Problem& problem,
-                                                   const SubControlVolume& scv,
-                                                   const VolumeVariables& volVars,
-                                                   CellCenterPrimaryVariables& storage) const
-    {
-        storage[Indices::energyBalanceIdx - cellCenterOffset] = volVars.density() * volVars.internalEnergy();
-    }
-
-    //! Evaluate energy storage for isothermal models
-    template <typename... Args>
-    void computeStorageForCellCenterNonIsothermal_(std::false_type, Args&&... args) const {}
-
 private:
 
     //! Returns the implementation of the problem (i.e. static polymorphism)
diff --git a/dumux/freeflow/navierstokes/volumevariables.hh b/dumux/freeflow/navierstokes/volumevariables.hh
index 87e1c41612b66160499025356e2df94a176e091b..8bc831ae3e194d0c09db20c3c518f02a71ac0e21 100644
--- a/dumux/freeflow/navierstokes/volumevariables.hh
+++ b/dumux/freeflow/navierstokes/volumevariables.hh
@@ -40,15 +40,16 @@ class NavierStokesVolumeVariables : public FreeFlowVolumeVariables< Traits, Navi
     using ParentType = FreeFlowVolumeVariables<Traits, ThisType>;
 
     using Scalar = typename Traits::PrimaryVariables::value_type;
-    using Indices = typename Traits::ModelTraits::Indices;
 
-    static constexpr int fluidSystemPhaseIdx = Indices::fluidSystemPhaseIdx;
+    static constexpr int fluidSystemPhaseIdx = Traits::ModelTraits::Indices::fluidSystemPhaseIdx;
 
 public:
     //! export the underlying fluid system
     using FluidSystem = typename Traits::FluidSystem;
     //! export the fluid state type
     using FluidState = typename Traits::FluidState;
+    //! export the indices type
+    using Indices = typename Traits::ModelTraits::Indices;
 
     /*!
      * \brief Update all quantities for a given control volume
diff --git a/dumux/freeflow/nonisothermal/localresidual.hh b/dumux/freeflow/nonisothermal/localresidual.hh
new file mode 100644
index 0000000000000000000000000000000000000000..ca84910ccea765e3c8e4efb2032c5497cfa50ba3
--- /dev/null
+++ b/dumux/freeflow/nonisothermal/localresidual.hh
@@ -0,0 +1,135 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ *   See the file COPYING for full copying permissions.                      *
+ *                                                                           *
+ *   This program is free software: you can redistribute it and/or modify    *
+ *   it under the terms of the GNU General Public License as published by    *
+ *   the Free Software Foundation, either version 2 of the License, or       *
+ *   (at your option) any later version.                                     *
+ *                                                                           *
+ *   This program is distributed in the hope that it will be useful,         *
+ *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
+ *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the            *
+ *   GNU General Public License for more details.                            *
+ *                                                                           *
+ *   You should have received a copy of the GNU General Public License       *
+ *   along with this program.  If not, see <http://www.gnu.org/licenses/>.   *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup FreeflowNIModel
+ * \copydoc Dumux::FreeFlowEnergyLocalResidual
+ */
+#ifndef DUMUX_FREE_FLOW_ENERGY_LOCAL_RESIDUAL_HH
+#define DUMUX_FREE_FLOW_ENERGY_LOCAL_RESIDUAL_HH
+
+#include <dumux/discretization/methods.hh>
+
+namespace Dumux {
+
+// forward declaration
+template<class FVGridGeometry, class FluxVariables, DiscretizationMethod discMethod, bool enableEneryBalance>
+class FreeFlowEnergyLocalResidualImplementation;
+
+/*!
+ * \ingroup FreeflowNIModel
+ * \brief Element-wise calculation of the local residual for non-isothermal
+ *        free-flow models
+ */
+template<class FVGridGeometry, class FluxVariables, bool enableEneryBalance>
+using FreeFlowEnergyLocalResidual =
+      FreeFlowEnergyLocalResidualImplementation<FVGridGeometry,
+                                                FluxVariables,
+                                                FVGridGeometry::discMethod,
+                                                enableEneryBalance>;
+
+/*!
+ * \ingroup FreeflowNIModel
+ * \brief Specialization for isothermal models, does nothing
+ */
+template<class FVGridGeometry, class FluxVariables, DiscretizationMethod discMethod>
+class FreeFlowEnergyLocalResidualImplementation<FVGridGeometry, FluxVariables, discMethod, false>
+{
+public:
+
+    //! do nothing for the isothermal case
+    template <typename... Args>
+    static void fluidPhaseStorage(Args&&... args)
+    {}
+
+    //! do nothing for the isothermal case
+    template <typename... Args>
+    static void heatConvectionFlux(Args&&... args)
+    {}
+
+    //! do nothing for the isothermal case
+    template <typename... Args>
+    static void heatFlux(Args&&... args)
+    {}
+};
+
+/*!
+ * \ingroup FreeflowNIModel
+ * \brief Specialization for staggered non-isothermal models
+ */
+template<class FVGridGeometry, class FluxVariables>
+class FreeFlowEnergyLocalResidualImplementation<FVGridGeometry,
+                                                FluxVariables,
+                                                DiscretizationMethod::staggered,
+                                                true>
+{
+    using Element = typename FVGridGeometry::GridView::template Codim<0>::Entity;
+    using FVElementGeometry = typename FVGridGeometry::LocalView;
+    using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
+    using HeatConductionType = typename FluxVariables::HeatConductionType;
+
+public:
+
+    //! The energy storage in the fluid phase
+    template<class NumEqVector, class VolumeVariables>
+    static void fluidPhaseStorage(NumEqVector& storage,
+                                  const VolumeVariables& volVars)
+    {
+        static constexpr auto localEnergyBalanceIdx = NumEqVector::dimension - 1;
+        storage[localEnergyBalanceIdx] += volVars.density() * volVars.internalEnergy();
+    }
+
+    //! The convective and conductive heat fluxes in the fluid phase
+    template<class NumEqVector, class Problem, class ElementVolumeVariables, class ElementFaceVariables>
+    static void heatFlux(NumEqVector& flux,
+                         const Problem& problem,
+                         const Element &element,
+                         const FVElementGeometry& fvGeometry,
+                         const ElementVolumeVariables& elemVolVars,
+                         const ElementFaceVariables& elemFaceVars,
+                         const SubControlVolumeFace& scvf)
+    {
+        static constexpr auto localEnergyBalanceIdx = NumEqVector::dimension - 1;
+        using Indices = typename ElementVolumeVariables::VolumeVariables::Indices;
+
+        bool isOutflow = false;
+        if(scvf.boundary())
+        {
+            const auto bcTypes = problem.boundaryTypesAtPos(scvf.center());
+            if(bcTypes.isOutflow(Indices::energyBalanceIdx))
+                isOutflow = true;
+        }
+
+        auto upwindTerm = [](const auto& volVars) { return volVars.density() * volVars.enthalpy(); };
+        flux[localEnergyBalanceIdx] += FluxVariables::advectiveFluxForCellCenter(elemVolVars,
+                                                                            elemFaceVars,
+                                                                            scvf,
+                                                                            upwindTerm,
+                                                                            isOutflow);
+        if(!isOutflow)
+            flux[localEnergyBalanceIdx] += HeatConductionType::flux(element,
+                                                               fvGeometry,
+                                                               elemVolVars,
+                                                               scvf);
+    }
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/dumux/freeflow/properties.hh b/dumux/freeflow/properties.hh
index 0138d3ba91d9706b022b95d7c263419703cfd3f7..772e450e0481e16344c9da8586ac8f72d3893ac9 100644
--- a/dumux/freeflow/properties.hh
+++ b/dumux/freeflow/properties.hh
@@ -28,6 +28,7 @@
 
 #include <dumux/common/properties.hh>
 #include <dumux/common/properties/model.hh>
+#include <dumux/discretization/fourierslaw.hh>
 
 namespace Dumux
 {
@@ -36,6 +37,9 @@ namespace Properties
 //! Type tag for free-flow models
 NEW_TYPE_TAG(FreeFlow, INHERITS_FROM(ModelProperties));
 
+//! Use Fourier's Law as default heat conduction type
+SET_TYPE_PROP(FreeFlow, HeatConductionType, FouriersLaw<TypeTag>);
+
 } // namespace Properties
 } // namespace Dumux
 
diff --git a/dumux/freeflow/rans/twoeq/lowrekepsilon/volumevariables.hh b/dumux/freeflow/rans/twoeq/lowrekepsilon/volumevariables.hh
index af2d6a34fb5234372b101dde5a1a7bea1fd0f9d3..7cc2e6a2a6758959f03da7b34881037ba3e5c7ed 100644
--- a/dumux/freeflow/rans/twoeq/lowrekepsilon/volumevariables.hh
+++ b/dumux/freeflow/rans/twoeq/lowrekepsilon/volumevariables.hh
@@ -47,14 +47,15 @@ class LowReKEpsilonVolumeVariables
     using NavierStokesParentType = NSVolumeVariables;
 
     using Scalar = typename Traits::PrimaryVariables::value_type;
-    using Indices = typename Traits::ModelTraits::Indices;
 
     static constexpr bool enableEnergyBalance = Traits::ModelTraits::enableEnergyBalance();
-    static constexpr int fluidSystemPhaseIdx = Indices::fluidSystemPhaseIdx;
+    static constexpr int fluidSystemPhaseIdx = Traits::ModelTraits::Indices::fluidSystemPhaseIdx;
 
 public:
     //! export the underlying fluid system
     using FluidSystem = typename Traits::FluidSystem;
+    //! export the indices type
+    using Indices = typename Traits::ModelTraits::Indices;
 
     /*!
      * \brief Update all quantities for a given control volume
diff --git a/dumux/freeflow/rans/zeroeq/volumevariables.hh b/dumux/freeflow/rans/zeroeq/volumevariables.hh
index 4c85ac557cde9667419ae72e9c6042d7554f1331..aabe2bcbce0c6bb245b3fc8d916c7a5fb203b39b 100644
--- a/dumux/freeflow/rans/zeroeq/volumevariables.hh
+++ b/dumux/freeflow/rans/zeroeq/volumevariables.hh
@@ -46,16 +46,17 @@ class ZeroEqVolumeVariables
     using NavierStokesParentType = NSVolumeVariables;
 
     using Scalar = typename Traits::PrimaryVariables::value_type;
-    using Indices = typename Traits::ModelTraits::Indices;
 
     static constexpr bool enableEnergyBalance = Traits::ModelTraits::enableEnergyBalance();
-    static constexpr int fluidSystemPhaseIdx = Indices::fluidSystemPhaseIdx;
+    static constexpr int fluidSystemPhaseIdx = Traits::ModelTraits::Indices::fluidSystemPhaseIdx;
 
 public:
     //! export the underlying fluid system
     using FluidSystem = typename Traits::FluidSystem;
     //! export the fluid state type
     using FluidState = typename Traits::FluidState;
+    //! export the indices type
+    using Indices = typename Traits::ModelTraits::Indices;
 
     /*!
      * \brief Update all quantities for a given control volume