diff --git a/dumux/porousmediumflow/2p/implicit/volumevariables.hh b/dumux/porousmediumflow/2p/implicit/volumevariables.hh
index c7fc61d2eb0e86347ed63de12a71c54f995e7423..0e90d7625366df7d1e539c3bbe78ec484824d582 100644
--- a/dumux/porousmediumflow/2p/implicit/volumevariables.hh
+++ b/dumux/porousmediumflow/2p/implicit/volumevariables.hh
@@ -103,8 +103,6 @@ public:
         porosity_ = problem.spatialParams().porosity(element, scv, elemSol);
         permeability_ = problem.spatialParams().permeability(element, scv, elemSol);
 
-        // energy related quantities not belonging to the fluid state
-        asImp_().updateEnergy_(elemSol, problem, element, scv);
     }
 
     /*!
@@ -116,7 +114,7 @@ public:
                                    const SubControlVolume& scv,
                                    FluidState& fluidState)
     {
-        Scalar t = Implementation::temperature_(elemSol, problem, element, scv);
+        Scalar t =  ParentType::temperature(elemSol, problem, element, scv);
         fluidState.setTemperature(t);
 
         const auto& materialParams = problem.spatialParams().materialLawParams(element, scv, elemSol);
@@ -156,7 +154,7 @@ public:
             fluidState.setDensity(phaseIdx, rho);
 
             // compute and set the enthalpy
-            Scalar h = Implementation::enthalpy_(fluidState, paramCache, phaseIdx);
+            Scalar h = ParentType::enthalpy(fluidState, paramCache, phaseIdx);
             fluidState.setEnthalpy(phaseIdx, h);
         }
     }
@@ -243,30 +241,6 @@ public:
     { return permeability_; }
 
 protected:
-    static Scalar temperature_(const ElementSolutionVector &elemSol,
-                               const Problem& problem,
-                               const Element &element,
-                               const SubControlVolume &scv)
-    {
-        return problem.temperatureAtPos(scv.dofPosition());
-    }
-
-    template<class ParameterCache>
-    static Scalar enthalpy_(const FluidState& fluidState,
-                            const ParameterCache& paramCache,
-                            const int phaseIdx)
-    {
-        return 0;
-    }
-
-    /*!
-     * \brief Called by update() to compute the energy related quantities.
-     */
-    void updateEnergy_(const ElementSolutionVector &elemSol,
-                       const Problem &problem,
-                       const Element &element,
-                       const SubControlVolume& scv)
-    {}
 
     FluidState fluidState_;
     Scalar porosity_;
diff --git a/test/porousmediumflow/2p/implicit/injectionproblem2pni.hh b/test/porousmediumflow/2p/implicit/injectionproblem2pni.hh
index eb5425cd02585522e78c3bb86029d873ab38af9d..db808df30fcfe262a4ad12aee2d460eca32f6f1d 100644
--- a/test/porousmediumflow/2p/implicit/injectionproblem2pni.hh
+++ b/test/porousmediumflow/2p/implicit/injectionproblem2pni.hh
@@ -29,9 +29,9 @@
 
 #include <dumux/porousmediumflow/2p/implicit/model.hh>
 #include <dumux/porousmediumflow/implicit/problem.hh>
-
+#include <dumux/implicit/cellcentered/tpfa/properties.hh>
 #include <dumux/material/fluidsystems/h2on2.hh>
-
+#include <dumux/material/components/n2.hh>
 // use the same spatial parameters as the injection problem of the
 // 2p2c test program
 #include "test/porousmediumflow/2p2c/implicit/injectionspatialparams.hh"
@@ -48,11 +48,11 @@ namespace Properties
 #if !ISOTHERMAL
 NEW_TYPE_TAG(InjectionProblem2PNI, INHERITS_FROM(TwoPNI, InjectionSpatialParams));
 NEW_TYPE_TAG(InjectionBoxProblem2PNI, INHERITS_FROM(BoxModel, InjectionProblem2PNI));
-NEW_TYPE_TAG(InjectionCCProblem2PNI, INHERITS_FROM(CCModel, InjectionProblem2PNI));
+NEW_TYPE_TAG(InjectionCCProblem2PNI, INHERITS_FROM(CCTpfaModel, InjectionProblem2PNI));
 #else
 NEW_TYPE_TAG(InjectionProblem2PNI, INHERITS_FROM(TwoP, InjectionSpatialParams));
 NEW_TYPE_TAG(InjectionBoxProblem2PNI, INHERITS_FROM(BoxModel, InjectionProblem2PNI));
-NEW_TYPE_TAG(InjectionCCProblem2PNI, INHERITS_FROM(CCModel, InjectionProblem2PNI));
+NEW_TYPE_TAG(InjectionCCProblem2PNI, INHERITS_FROM(CCTpfaModel, InjectionProblem2PNI));
 #endif
 
 // Set the grid type
@@ -118,11 +118,8 @@ class InjectionProblem2PNI : public ImplicitPorousMediaProblem<TypeTag>
     typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
     typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
 
-#if ISOTHERMAL
-    typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
-#else
     typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
-#endif
+
     enum {
         pressureIdx = Indices::pressureIdx,
         saturationIdx = Indices::saturationIdx,
@@ -139,7 +136,15 @@ class InjectionProblem2PNI : public ImplicitPorousMediaProblem<TypeTag>
     };
 
 
-    typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
+    using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+    using NeumannFluxes = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+    using Sources = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+    using ThermalConductivityModel = typename GET_PROP_TYPE(TypeTag, ThermalConductivityModel);
+    using GasN2 = N2<Scalar>;
+
+    using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
+    using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+
     typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
     typedef typename GET_PROP_TYPE(TypeTag, TimeManager) TimeManager;
 
@@ -207,10 +212,11 @@ public:
      *               \f$ [ \textnormal{unit of primary variable} / (m^\textrm{dim} \cdot s )] \f$
      * \param globalPos The global position
      */
-    void sourceAtPos(PrimaryVariables &values,
-                const GlobalPosition &globalPos) const
+    Sources sourceAtPos(const GlobalPosition &globalPos) const
     {
+        Sources values(0.0);
         values = 0;
+        return values;
     }
 
     // \}
@@ -227,19 +233,18 @@ public:
      * \param values Stores the value of the boundary type
      * \param globalPos The global position
      */
-    void boundaryTypesAtPos(BoundaryTypes &values,
-                            const GlobalPosition &globalPos) const
+    BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
     {
+        BoundaryTypes values;
         if (globalPos[0] < eps_)
             values.setAllDirichlet();
         else
             values.setAllNeumann();
-
+        return values;
 #if !ISOTHERMAL
-        // set a dirichlet value for the temperature, use the energy
-        // equation to set the value
         values.setDirichlet(temperatureIdx);
 #endif
+
     }
 
     /*!
@@ -250,14 +255,16 @@ public:
      *               \f$ [ \textnormal{unit of primary variable} ] \f$
      * \param globalPos The global position
      */
-    void dirichletAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const
+     PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
     {
+        PrimaryVariables values;
         Scalar densityW = 1000.0;
         values[pressureIdx] = 1e5 + (maxDepth_ - globalPos[1])*densityW*9.81;
         values[saturationIdx] = 0.0;
 #if !ISOTHERMAL
         values[temperatureIdx] = 283.0 + (maxDepth_ - globalPos[1])*0.03;
 #endif
+         return values;
     }
 
     /*!
@@ -275,25 +282,20 @@ public:
      * The \a values store the mass flux of each phase normal to the boundary.
      * Negative values indicate an inflow.
      */
-    void neumann(PrimaryVariables &values,
-                 const Element &element,
-                 const FVElementGeometry &fvGeometry,
-                 const Intersection &intersection,
-                 int scvIdx,
-                 int boundaryFaceIdx) const
+    NeumannFluxes neumann(const Element &element,
+                             const FVElementGeometry& fvGeometry,
+                             const ElementVolumeVariables& elemVolVars,
+                             const SubControlVolumeFace& scvf) const
     {
-        values = 0;
-
-        GlobalPosition globalPos;
-        if (isBox)
-            globalPos = element.geometry().corner(scvIdx);
-        else
-            globalPos = intersection.geometry().center();
+        NeumannFluxes values(0.0);
+        const auto globalPos = scvf.ipGlobal();
 
         if (globalPos[1] < 15 + eps_ && globalPos[1] > 7 - eps_) {
             // inject air. negative values mean injection
             values[contiNEqIdx] = -1e-3; // kg/(s*m^2)
+
         }
+         return values;
     }
 
     // \}
@@ -311,8 +313,9 @@ public:
      *               \f$ [ \textnormal{unit of primary variables} ] \f$
      * \param globalPos The global position
      */
-    void initialAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const
+    PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
     {
+        PrimaryVariables values;
         Scalar densityW = 1000.0;
         values[pressureIdx] = 1e5 + (maxDepth_ - globalPos[1])*densityW*9.81;
         values[saturationIdx] = 0.0;
@@ -322,6 +325,7 @@ public:
         if (globalPos[0] > 20 - eps_ && globalPos[0] < 30 + eps_ && globalPos[1] > 5 - eps_ && globalPos[1] < 35 + eps_)
             values[temperatureIdx] = 380;
 #endif // !ISOTHERMAL
+         return values;
     }
     // \}
 
diff --git a/test/porousmediumflow/2p/implicit/lensproblem.hh b/test/porousmediumflow/2p/implicit/lensproblem.hh
index 1e7d797f9f474214365e05a81b9a4bc48f842bc9..babb81109d5c4e597319a6f876545c6a78768361 100644
--- a/test/porousmediumflow/2p/implicit/lensproblem.hh
+++ b/test/porousmediumflow/2p/implicit/lensproblem.hh
@@ -193,6 +193,8 @@ class LensProblem : public GET_PROP_TYPE(TypeTag, BaseProblem)
     enum { adaptiveGrid = GET_PROP_VALUE(TypeTag, AdaptiveGrid) };
 
     typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
+    using NeumannFluxes = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+    using Sources = typename GET_PROP_TYPE(TypeTag, NumEqVector);
     typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
     typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
     typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
@@ -275,7 +277,7 @@ public:
      *               \f$ [ \textnormal{unit of primary variable} / (m^\textrm{dim} \cdot s )] \f$
      * \param globalPos The global position
      */
-    PrimaryVariables sourceAtPos(const GlobalPosition &globalPos) const
+    Sources sourceAtPos(const GlobalPosition &globalPos) const
     {
         return PrimaryVariables(0.0);
     }
@@ -348,9 +350,9 @@ public:
      * For this method, the \a values parameter stores the mass flux
      * in normal direction of each phase. Negative values mean influx.
      */
-    PrimaryVariables neumannAtPos(const GlobalPosition &globalPos) const
+    NeumannFluxes neumannAtPos(const GlobalPosition &globalPos) const
     {
-        PrimaryVariables values(0.0);
+        NeumannFluxes values(0.0);
         if (onInlet_(globalPos)) {
             values[contiNEqIdx] = -0.04; // kg / (m * s)
         }