diff --git a/dumux/porousmediumflow/3p3c/implicit/model.hh b/dumux/porousmediumflow/3p3c/implicit/model.hh
index 15f401091437774ae3ff4b846e21b28d94773650..63bb7080d2d0b7abf4277fef1972a6ba5a82c281 100644
--- a/dumux/porousmediumflow/3p3c/implicit/model.hh
+++ b/dumux/porousmediumflow/3p3c/implicit/model.hh
@@ -164,11 +164,15 @@ public:
         vtkOutputModule.addSecondaryVariable("rhon", [](const VolumeVariables& v){ return v.density(nPhaseIdx); });
         vtkOutputModule.addSecondaryVariable("rhog", [](const VolumeVariables& v){ return v.density(gPhaseIdx); });
 
+
         for (int i = 0; i < numPhases; ++i)
             for (int j = 0; j < numComponents; ++j)
                 vtkOutputModule.addSecondaryVariable("x^" + FluidSystem::componentName(j) + "_" + FluidSystem::phaseName(i),
                                                      [i,j](const VolumeVariables& v){ return v.moleFraction(i,j); });
 
+        vtkOutputModule.addSecondaryVariable("porosity", [](const VolumeVariables& v){ return v.porosity(); });
+        vtkOutputModule.addSecondaryVariable("permeability",
+                                             [](const VolumeVariables& v){ return v.permeability(); });
         vtkOutputModule.addSecondaryVariable("temperature", [](const VolumeVariables& v){ return v.temperature(); });
     }
 
diff --git a/dumux/porousmediumflow/3p3c/implicit/properties.hh b/dumux/porousmediumflow/3p3c/implicit/properties.hh
index 3c74683eb3821c7cf5936f8abb8d35615ae54201..84dc9e7986d675f02a753d1ef0a27abd322b573b 100644
--- a/dumux/porousmediumflow/3p3c/implicit/properties.hh
+++ b/dumux/porousmediumflow/3p3c/implicit/properties.hh
@@ -73,6 +73,8 @@ NEW_PROP_TAG(ProblemEnableGravity); //!< Returns whether gravity is considered i
 NEW_PROP_TAG(UseConstraintSolver); //!< Determines whether a constraint solver should be used explicitly
 NEW_PROP_TAG(SpatialParamsForchCoeff); //!< Property for the forchheimer coefficient
 NEW_PROP_TAG(TauTortuosity); //!< Tortuosity value (tau) used in macroscopic diffusion
+NEW_PROP_TAG(UseMoles);//!< Defines whether mole (true) or mass (false) fractions are used
+NEW_PROP_TAG(ReplaceCompEqIdx); //!< The index of the total mass balance equation, if one component balance is replaced (ReplaceCompEqIdx < NumComponents)
 }
 }
 
diff --git a/dumux/porousmediumflow/3p3c/implicit/propertydefaults.hh b/dumux/porousmediumflow/3p3c/implicit/propertydefaults.hh
index f5542001a02ff3885f14fb40b659efa4a61f6738..a1ff05ff0ea6efff7926c68a741c350253748bb4 100644
--- a/dumux/porousmediumflow/3p3c/implicit/propertydefaults.hh
+++ b/dumux/porousmediumflow/3p3c/implicit/propertydefaults.hh
@@ -34,10 +34,10 @@
 #include "volumevariables.hh"
 #include "properties.hh"
 #include "newtoncontroller.hh"
-#include "localresidual.hh"
 #include "primaryvariableswitch.hh"
 
 #include <dumux/porousmediumflow/compositional/switchableprimaryvariables.hh>
+#include <dumux/porousmediumflow/compositional/localresidual.hh>
 #include <dumux/porousmediumflow/nonisothermal/implicit/propertydefaults.hh>
 #include <dumux/material/fluidmatrixinteractions/diffusivitymillingtonquirk.hh>
 #include <dumux/porousmediumflow/implicit/darcyfluxvariables.hh>
@@ -87,6 +87,8 @@ SET_PROP(ThreePThreeC, NumPhases)
                   "Only fluid systems with 3 phases are supported by the 3p3c model!");
 };
 
+//! Set as default that no component mass balance is replaced by the total mass balance
+SET_INT_PROP(ThreePThreeC, ReplaceCompEqIdx, 100);
 /*!
  * \brief The fluid state which is used by the volume variables to
  *        store the thermodynamic state. This should be chosen
@@ -110,7 +112,7 @@ SET_INT_PROP(ThreePThreeC, NumEq, 3); //!< set the number of equations to 2
 SET_TYPE_PROP(ThreePThreeC, MaterialLawParams, typename GET_PROP_TYPE(TypeTag, MaterialLaw)::Params);
 
 //! The local residual function of the conservation equations
-SET_TYPE_PROP(ThreePThreeC, LocalResidual, ThreePThreeCLocalResidual<TypeTag>);
+SET_TYPE_PROP(ThreePThreeC, LocalResidual, CompositionalLocalResidual<TypeTag>);
 
 //! Enable advection
 SET_BOOL_PROP(ThreePThreeC, EnableAdvection, true);
@@ -174,6 +176,9 @@ SET_SCALAR_PROP(BoxModel, SpatialParamsForchCoeff, 0.55);
  */
 SET_SCALAR_PROP(ThreePThreeC, TauTortuosity, 0.5);
 
+//! Use mole fractions in the balance equations by default
+SET_BOOL_PROP(ThreePThreeC, UseMoles, true);
+
 //! Somerton is used as default model to compute the effective thermal heat conductivity
 SET_PROP(ThreePThreeCNI, ThermalConductivityModel)
 {
@@ -198,7 +203,7 @@ SET_TYPE_PROP(ThreePThreeCNI, IsothermalModel, ThreePThreeCModel<TypeTag>);
 SET_TYPE_PROP(ThreePThreeCNI, IsothermalVolumeVariables, ThreePThreeCVolumeVariables<TypeTag>);
 
 //set isothermal LocalResidual
-SET_TYPE_PROP(ThreePThreeCNI, IsothermalLocalResidual, ThreePThreeCLocalResidual<TypeTag>);
+SET_TYPE_PROP(ThreePThreeCNI, IsothermalLocalResidual, CompositionalLocalResidual<TypeTag>);
 
 //set isothermal Indices
 SET_PROP(ThreePThreeCNI, IsothermalIndices)
diff --git a/dumux/porousmediumflow/3p3c/implicit/volumevariables.hh b/dumux/porousmediumflow/3p3c/implicit/volumevariables.hh
index 1168f2808e6b350c96692dc207800f9a3ec2c4d9..54537f004c6cd9576b6358c12d76973af0f16650 100644
--- a/dumux/porousmediumflow/3p3c/implicit/volumevariables.hh
+++ b/dumux/porousmediumflow/3p3c/implicit/volumevariables.hh
@@ -56,6 +56,7 @@ class ThreePThreeCVolumeVariables : public ImplicitVolumeVariables<TypeTag>
     using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
     using MaterialLaw = typename GET_PROP_TYPE(TypeTag, MaterialLaw);
     using MaterialLawParams = typename GET_PROP_TYPE(TypeTag, MaterialLawParams);
+    using ElementSolutionVector = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
 
     // constraint solvers
     using SpatialParams = typename GET_PROP_TYPE(TypeTag, SpatialParams);
@@ -122,6 +123,7 @@ public:
         const MaterialLawParams &materialParams =
             problem.spatialParams().materialLawParams(element, scv, elemSol);
 
+
         Scalar temp = ParentType::temperature(elemSol, problem, element, scv);
         fluidState_.setTemperature(temp);
 
@@ -703,6 +705,8 @@ private:
             diffCoefficient_[phaseIdx][compIdx] = std::move(d);
         else if (compIdx > phaseIdx)
             diffCoefficient_[phaseIdx][compIdx-1] = std::move(d);
+        else if (phaseIdx == nPhaseIdx)
+            diffCoefficient_[phaseIdx][compIdx-1] = 0;
         else
             DUNE_THROW(Dune::InvalidStateException, "Diffusion coeffiecient for phaseIdx = compIdx doesn't exist");
     }
diff --git a/test/porousmediumflow/3p3c/implicit/columnxylolproblem.hh b/test/porousmediumflow/3p3c/implicit/columnxylolproblem.hh
index 0665d375f11e8b36d828f789c431faf51f2faa2e..c9df17677848c5e44a1fbef78a54ce2d082076dd 100644
--- a/test/porousmediumflow/3p3c/implicit/columnxylolproblem.hh
+++ b/test/porousmediumflow/3p3c/implicit/columnxylolproblem.hh
@@ -26,7 +26,7 @@
 #define DUMUX_COLUMNXYLOLPROBLEM_HH
 
 #include <dumux/material/fluidsystems/h2oairxylene.hh>
-
+#include <dumux/implicit/cellcentered/tpfa/properties.hh>
 #include <dumux/porousmediumflow/3p3c/implicit/model.hh>
 #include <dumux/porousmediumflow/implicit/problem.hh>
 
@@ -43,7 +43,7 @@ namespace Properties
 {
 NEW_TYPE_TAG(ColumnProblem, INHERITS_FROM(ThreePThreeCNI, ColumnSpatialParams));
 NEW_TYPE_TAG(ColumnBoxProblem, INHERITS_FROM(BoxModel, ColumnProblem));
-NEW_TYPE_TAG(ColumnCCProblem, INHERITS_FROM(CCModel, ColumnProblem));
+NEW_TYPE_TAG(ColumnCCProblem, INHERITS_FROM(CCTpfaModel, ColumnProblem));
 
 // Set the grid type
 SET_TYPE_PROP(ColumnProblem, Grid, Dune::YaspGrid<2>);
@@ -113,6 +113,7 @@ class ColumnProblem : public ImplicitPorousMediaProblem<TypeTag>
 
 
     typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
+    typedef typename GET_PROP_TYPE(TypeTag, NumEqVector) NeumannFluxes;
     typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
     typedef typename GET_PROP_TYPE(TypeTag, TimeManager) TimeManager;
 
@@ -122,6 +123,9 @@ class ColumnProblem : public ImplicitPorousMediaProblem<TypeTag>
 
     typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
     typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+    using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
+    using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
+    using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
 
     typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
 
@@ -156,11 +160,16 @@ public:
     const std::string name() const
     { return name_; }
 
-
-    void sourceAtPos(PrimaryVariables &values,
-                     const GlobalPosition &globalPos) const
+    /*!
+     * \brief Returns the source term
+     *
+     * \param values Stores the source values for the conservation equations in
+     *               \f$ [ \textnormal{unit of primary variable} / (m^\textrm{dim} \cdot s )] \f$
+     * \param globalPos The global position
+     */
+    PrimaryVariables sourceAtPos(const GlobalPosition &globalPos) const
     {
-        values = 0;
+        return PrimaryVariables(0.0);
     }
 
     // \}
@@ -177,13 +186,14 @@ public:
      * \param values The boundary types for the conservation equations
      * \param globalPos The position for which the bc type should be evaluated
      */
-    void boundaryTypesAtPos(BoundaryTypes &values,
-                            const GlobalPosition &globalPos) const
+    BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
     {
+        BoundaryTypes bcTypes;
         if (globalPos[1] < eps_)
-            values.setAllDirichlet();
+            bcTypes.setAllDirichlet();
         else
-            values.setAllNeumann();
+            bcTypes.setAllNeumann();
+        return bcTypes;
     }
 
     /*!
@@ -195,17 +205,17 @@ public:
      *
      * For this method, the \a values parameter stores primary variables.
      */
-    void dirichletAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const
+    PrimaryVariables dirichletAtPos( const GlobalPosition &globalPos) const
     {
-        initial_(values, globalPos);
+       return initial_(globalPos);
+
     }
 
     /*!
      * \brief Evaluate the boundary conditions for a neumann
      *        boundary segment.
      *
-     * \param values The neumann values for the conservation equations
-     * \param element The finite element
+      * \param element The finite element
      * \param fvGeomtry The finite-volume geometry in the box scheme
      * \param intersection The intersection between element and boundary
      * \param scvIdx The local vertex index
@@ -214,20 +224,13 @@ public:
      * For this method, the \a values parameter stores the mass flux
      * in normal direction of each phase. Negative values mean influx.
      */
-    void neumann(PrimaryVariables &values,
-                 const Element &element,
-                 const FVElementGeometry &fvGeomtry,
-                 const Intersection &intersection,
-                 const 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();
+        PrimaryVariables values(0.0);
+        const auto& globalPos = scvf.ipGlobal();
 
         // negative values for injection
         if (globalPos[1] > 1.2 - eps_)
@@ -237,14 +240,11 @@ public:
             values[Indices::contiNEqIdx] = -0.00;
             values[Indices::energyEqIdx] = -17452.97;
         }
+        return values;
     }
 
     // \}
 
-    /*!
-     * \name Volume terms
-     */
-    // \{
 
     /*!
      * \brief Evaluate the initial value for a control volume.
@@ -255,19 +255,11 @@ public:
      * For this method, the \a values parameter stores primary
      * variables.
      */
-    void initialAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const
+    PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
     {
-        initial_(values, globalPos);
+        return initial_(globalPos);
     }
 
-    /*!
-     * \brief Evaluate the initial phase state at a given position
-     *
-     * \param globalPos The global position
-     */
-    int initialPhasePresenceAtPos(const GlobalPosition &globalPos)
-    { return threePhases; }
-
 
        /*!
      * \brief Append all quantities of interest which can be derived
@@ -302,9 +294,10 @@ public:
 private:
     // internal method for the initial condition (reused for the
     // dirichlet conditions!)
-    void initial_(PrimaryVariables &values,
-                  const GlobalPosition &globalPos) const
+    PrimaryVariables initial_(const GlobalPosition &globalPos) const
     {
+        PrimaryVariables values;
+        values.setState(threePhases);
         Scalar y = globalPos[1];
 
         values[temperatureIdx] = 296.15;
@@ -355,6 +348,7 @@ private:
             values[switch2Idx] = 0.988;
         else
             values[switch2Idx] = 1.e-4;
+        return values;
     }
 
     static constexpr Scalar eps_ = 1e-6;
diff --git a/test/porousmediumflow/3p3c/implicit/columnxylolspatialparams.hh b/test/porousmediumflow/3p3c/implicit/columnxylolspatialparams.hh
index 3aa624a6f830df7e5b350d42431a878235880824..d494c61fe812610566078d779b8f5c9bd4873ed2 100644
--- a/test/porousmediumflow/3p3c/implicit/columnxylolspatialparams.hh
+++ b/test/porousmediumflow/3p3c/implicit/columnxylolspatialparams.hh
@@ -95,19 +95,22 @@ class ColumnSpatialParams : public ImplicitSpatialParams<TypeTag>
 
     typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
     typedef typename GridView::template Codim<0>::Entity Element;
-
+    using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
+    using ElementSolutionVector = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
+    using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
 
 public:
     typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
     typedef typename MaterialLaw::Params MaterialLawParams;
+    using PermeabilityType = Scalar;
 
     /*!
      * \brief The constructor
      *
      * \param gridView The grid view
      */
-    ColumnSpatialParams(const GridView &gridView)
-        : ParentType(gridView)
+    ColumnSpatialParams(const Problem& problem, const GridView &gridView)
+        : ParentType(problem, gridView)
     {
         // intrinsic permeabilities
         fineK_ = 1.4e-11;
@@ -166,15 +169,10 @@ public:
      * \brief Apply the intrinsic permeability tensor to a pressure
      *        potential gradient.
      *
-     * \param element The current finite element
-     * \param fvGeometry The current finite volume geometry of the element
-     * \param scvIdx The index of the sub-control volume
+     * \param globalPos The global position
      */
-    const Scalar intrinsicPermeability(const Element &element,
-                                       const FVElementGeometry &fvGeometry,
-                                       int scvIdx) const
+   Scalar permeabilityAtPos(const GlobalPosition& globalPos) const
     {
-        const GlobalPosition &globalPos = fvGeometry.subContVol[scvIdx].global;
         if (isFineMaterial_(globalPos))
             return fineK_;
         return coarseK_;
@@ -188,11 +186,9 @@ public:
      * \param scvIdx The local index of the sub-control volume where
      *                    the porosity needs to be defined
      */
-    double porosity(const Element &element,
-                    const FVElementGeometry &fvGeometry,
-                    const int scvIdx) const
+    Scalar porosityAtPos(const GlobalPosition& globalPos) const
     {
-        const GlobalPosition &globalPos = fvGeometry.subContVol[scvIdx].global;
+
         if (isFineMaterial_(globalPos))
             return finePorosity_;
         else
@@ -203,15 +199,11 @@ public:
     /*!
      * \brief return the parameter object for the Brooks-Corey material law which depends on the position
      *
-     * \param element The current finite element
-     * \param fvGeometry The current finite volume geometry of the element
-     * \param scvIdx The index of the sub-control volume
+     * \param globalPos The global position
      */
-    const MaterialLawParams& materialLawParams(const Element &element,
-                                               const FVElementGeometry &fvGeometry,
-                                               const int scvIdx) const
+    const MaterialLawParams& materialLawParamsAtPos(const GlobalPosition& globalPos) const
     {
-        const GlobalPosition &globalPos = fvGeometry.subContVol[scvIdx].global;
+
         if (isFineMaterial_(globalPos))
             return fineMaterialParams_;
         else
@@ -227,11 +219,8 @@ public:
      * \param fvGeometry The finite volume geometry
      * \param scvIdx The local index of the sub-control volume
      */
-    Scalar solidHeatCapacity(const Element &element,
-                             const FVElementGeometry &fvGeometry,
-                             const int scvIdx) const
+    Scalar solidHeatCapacityAtPos(const GlobalPosition& globalPos) const
     {
-        const GlobalPosition &globalPos = fvGeometry.subContVol[scvIdx].global;
         if (isFineMaterial_(globalPos))
             return fineHeatCap_;
         else
@@ -248,8 +237,8 @@ public:
      * \param scvIdx The local index of the sub-control volume
      */
     Scalar solidDensity(const Element &element,
-                        const FVElementGeometry &fvGeometry,
-                        const int scvIdx) const
+                        const SubControlVolume& scv,
+                        const ElementSolutionVector& elemSol) const
     {
         return 2650; // density of sand [kg/m^3]
     }
@@ -263,8 +252,8 @@ public:
      *                    the heat capacity needs to be defined
      */
     Scalar solidThermalConductivity(const Element &element,
-                                    const FVElementGeometry &fvGeometry,
-                                    const int scvIdx) const
+                                    const SubControlVolume& scv,
+                                    const ElementSolutionVector& elemSol) const
     {
         return lambdaSolid_;
     }
diff --git a/test/porousmediumflow/3p3c/implicit/infiltrationspatialparameters.hh b/test/porousmediumflow/3p3c/implicit/infiltrationspatialparameters.hh
index ba52bc0497a7984990188014d3a9c16c0b39c8b5..9853ffc28ffbc73ff0477b0c4260b74bc791df04 100644
--- a/test/porousmediumflow/3p3c/implicit/infiltrationspatialparameters.hh
+++ b/test/porousmediumflow/3p3c/implicit/infiltrationspatialparameters.hh
@@ -101,7 +101,8 @@ class InfiltrationSpatialParams : public ImplicitSpatialParams<TypeTag>
 
 
 public:
-    using PermeabilityType = Scalar;
+
+   using PermeabilityType = Scalar;
 
     /*!
      * \brief The constructor
@@ -138,6 +139,7 @@ public:
      *
      * \param globalPos The global position
      */
+
     PermeabilityType permeabilityAtPos(const GlobalPosition& globalPos) const
     {
         if (isFineMaterial_(globalPos))
@@ -161,9 +163,7 @@ public:
     /*!
      * \brief return the parameter object for the material law which depends on the position
      *
-     * \param element The current finite element
-     * \param fvGeometry The current finite volume geometry of the element
-     * \param scvIdx The index of the sub-control volume
+     *  \param globalPos The global position
      */
     const MaterialLawParams& materialLawParamsAtPos(const GlobalPosition& globalPos) const
     {
diff --git a/test/porousmediumflow/3p3c/implicit/kuevetteproblem.hh b/test/porousmediumflow/3p3c/implicit/kuevetteproblem.hh
index 82a3b8a190f464701d8db7b1f5fcc7202058bf0c..1a4458334d3a6b0a524db3dbd563ba53e4ae6755 100644
--- a/test/porousmediumflow/3p3c/implicit/kuevetteproblem.hh
+++ b/test/porousmediumflow/3p3c/implicit/kuevetteproblem.hh
@@ -29,7 +29,7 @@
 #include <dune/common/float_cmp.hh>
 
 #include <dumux/material/fluidsystems/h2oairmesitylene.hh>
-
+#include <dumux/implicit/cellcentered/tpfa/properties.hh>
 #include <dumux/porousmediumflow/3p3c/implicit/model.hh>
 #include <dumux/porousmediumflow/implicit/problem.hh>
 
@@ -46,7 +46,7 @@ namespace Properties
 {
 NEW_TYPE_TAG(KuevetteProblem, INHERITS_FROM(ThreePThreeCNI, KuevetteSpatialParams));
 NEW_TYPE_TAG(KuevetteBoxProblem, INHERITS_FROM(BoxModel, KuevetteProblem));
-NEW_TYPE_TAG(KuevetteCCProblem, INHERITS_FROM(CCModel, KuevetteProblem));
+NEW_TYPE_TAG(KuevetteCCProblem, INHERITS_FROM(CCTpfaModel, KuevetteProblem));
 
 // Set the grid type
 SET_TYPE_PROP(KuevetteProblem, Grid, Dune::YaspGrid<2>);
@@ -127,6 +127,7 @@ class KuevetteProblem : public ImplicitPorousMediaProblem<TypeTag>
 
 
     typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
+    typedef typename GET_PROP_TYPE(TypeTag, NumEqVector) NeumannFluxes;
     typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
     typedef typename GET_PROP_TYPE(TypeTag, TimeManager) TimeManager;
 
@@ -137,6 +138,10 @@ class KuevetteProblem : public ImplicitPorousMediaProblem<TypeTag>
     typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
     typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
 
+    using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
+    using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
+    using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+
     typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
 
     enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
@@ -178,13 +183,13 @@ public:
      * \param values The source values for the primary variables
      * \param globalPos The position of the center of the finite volume
      */
-    void sourceAtPos(PrimaryVariables &values,
-                     const GlobalPosition &globalPos) const
+    PrimaryVariables sourceAtPos(const GlobalPosition &globalPos) const
     {
-        values = 0;
+        return PrimaryVariables(0.0);
     }
 
 
+
     // \}
 
     /*!
@@ -199,13 +204,14 @@ public:
      * \param values The boundary types for the conservation equations
      * \param globalPos The position for which the bc type should be evaluated
      */
-    void boundaryTypesAtPos(BoundaryTypes &values,
-                            const GlobalPosition &globalPos) const
+    BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
     {
+        BoundaryTypes bcTypes;
         if(globalPos[0] > this->bBoxMax()[0] - eps_)
-            values.setAllDirichlet();
+            bcTypes.setAllDirichlet();
         else
-            values.setAllNeumann();
+            bcTypes.setAllNeumann();
+         return bcTypes;
     }
 
     /*!
@@ -217,24 +223,31 @@ public:
      *
      * For this method, the \a values parameter stores primary variables.
      */
-    void dirichletAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const
+    PrimaryVariables dirichletAtPos( const GlobalPosition &globalPos) const
     {
-        initial_(values, globalPos);
+       return initial_(globalPos);
     }
 
     /*!
      * \brief Evaluate the boundary conditions for a neumann
      *        boundary segment.
      *
-     * \param values The neumann values for the primary variables
-     * \param globalPos The position for which the bc type should be evaluated
+     * \param element The finite element
+     * \param fvGeomtry The finite-volume geometry in the box scheme
+     * \param intersection The intersection between element and boundary
+     * \param scvIdx The local vertex index
+     * \param boundaryFaceIdx The index of the boundary face
      *
      * For this method, the \a values parameter stores the mass flux
      * in normal direction of each phase. Negative values mean influx.
      */
-    void neumannAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const
+    NeumannFluxes neumann(const Element& element,
+                             const FVElementGeometry& fvGeometry,
+                             const ElementVolumeVariables& elemVolVars,
+                             const SubControlVolumeFace& scvf) const
     {
-        values = 0;
+        PrimaryVariables values(0.0);
+        const auto& globalPos = scvf.ipGlobal();
 
         // negative values for injection
         if (globalPos[0] < eps_)
@@ -244,6 +257,7 @@ public:
             values[Indices::contiNEqIdx] =  0.0;
             values[Indices::energyEqIdx] = -6929.;
         }
+        return values;
     }
 
     // \}
@@ -262,22 +276,9 @@ public:
      * For this method, the \a values parameter stores primary
      * variables.
      */
-    void initialAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const
+    PrimaryVariables initialAtPos( const GlobalPosition &globalPos) const
     {
-        initial_(values, globalPos);
-    }
-
-    /*!
-     * \brief Evaluate the initial phase state at a given position
-     *
-     * \param globalPos The global position
-     */
-    int initialPhasePresenceAtPos(const GlobalPosition &globalPos)
-    {
-        if (isInContaminationZone(globalPos))
-            return threePhases;
-        else
-            return wgPhaseOnly;
+       return initial_(globalPos);
     }
 
     /*!
@@ -339,9 +340,14 @@ private:
 
     // internal method for the initial condition (reused for the
     // dirichlet conditions!)
-    void initial_(PrimaryVariables &values,
-                  const GlobalPosition &globalPos) const
+    PrimaryVariables initial_(const GlobalPosition &globalPos) const
     {
+        PrimaryVariables values;
+        if (isInContaminationZone(globalPos))
+            values.setState(threePhases);
+        else
+         values.setState(wgPhaseOnly);
+
         values[pressureIdx] = 1e5 ;
         values[switch1Idx] = 0.12;
         values[switch2Idx] = 1.e-6;
@@ -351,6 +357,7 @@ private:
         {
             values[switch2Idx] = 0.07;
         }
+         return values;
     }
 
     static constexpr Scalar eps_ = 1e-6;
diff --git a/test/porousmediumflow/3p3c/implicit/kuevettespatialparams.hh b/test/porousmediumflow/3p3c/implicit/kuevettespatialparams.hh
index 8e758f4688abd9f8143e5057cb560d924baf96d4..946ee1a21958ea0349e9ef982cdd8c78fe849070 100644
--- a/test/porousmediumflow/3p3c/implicit/kuevettespatialparams.hh
+++ b/test/porousmediumflow/3p3c/implicit/kuevettespatialparams.hh
@@ -97,19 +97,23 @@ class KuevetteSpatialParams : public ImplicitSpatialParams<TypeTag>
 
     typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
     typedef typename GridView::template Codim<0>::Entity Element;
+    using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
+    using ElementSolutionVector = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
+    using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
 
 
 public:
     typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
     typedef typename MaterialLaw::Params MaterialLawParams;
+    using PermeabilityType = Scalar;
 
     /*!
      * \brief The constructor
      *
      * \param gridView The grid view
      */
-    KuevetteSpatialParams(const GridView &gridView)
-        : ParentType(gridView)
+    KuevetteSpatialParams(const Problem& problem, const GridView &gridView)
+        : ParentType(problem, gridView)
     {
         // intrinsic permeabilities
         fineK_ = 6.28e-12;
@@ -164,15 +168,10 @@ public:
      * \brief Apply the intrinsic permeability tensor to a pressure
      *        potential gradient.
      *
-     * \param element The current finite element
-     * \param fvGeometry The current finite volume geometry of the element
-     * \param scvIdx The index of the sub-control volume
+     * \param globalPos The global position
      */
-    const Scalar intrinsicPermeability(const Element &element,
-                                       const FVElementGeometry &fvGeometry,
-                                       const int scvIdx) const
+   Scalar permeabilityAtPos(const GlobalPosition& globalPos) const
     {
-        const GlobalPosition &globalPos = fvGeometry.subContVol[scvIdx].global;
         if (isFineMaterial_(globalPos))
             return fineK_;
         return coarseK_;
@@ -186,11 +185,8 @@ public:
      * \param scvIdx The local index of the sub-control volume where
      *                    the porosity needs to be defined
      */
-    double porosity(const Element &element,
-                    const FVElementGeometry &fvGeometry,
-                    const int scvIdx) const
+    Scalar porosityAtPos(const GlobalPosition& globalPos) const
     {
-        const GlobalPosition &globalPos = fvGeometry.subContVol[scvIdx].global;
         if (isFineMaterial_(globalPos))
             return finePorosity_;
         else
@@ -205,11 +201,8 @@ public:
      * \param fvGeometry The current finite volume geometry of the element
      * \param scvIdx The index of the sub-control volume
      */
-    const MaterialLawParams& materialLawParams(const Element &element,
-                                               const FVElementGeometry &fvGeometry,
-                                               const int scvIdx) const
+    const MaterialLawParams& materialLawParamsAtPos(const GlobalPosition& globalPos) const
     {
-        const GlobalPosition &globalPos = fvGeometry.subContVol[scvIdx].global;
         if (isFineMaterial_(globalPos))
             return fineMaterialParams_;
         else
@@ -225,9 +218,7 @@ public:
      * \param fvGeometry The finite volume geometry
      * \param scvIdx The local index of the sub-control volume
      */
-    Scalar solidHeatCapacity(const Element &element,
-                             const FVElementGeometry &fvGeometry,
-                             const int scvIdx) const
+    Scalar solidHeatCapacityAtPos(const GlobalPosition& globalPos) const
     {
         return 850; // specific heat capacity of sand [J / (kg K)]
     }
@@ -242,8 +233,8 @@ public:
      * \param scvIdx The local index of the sub-control volume
      */
     Scalar solidDensity(const Element &element,
-                        const FVElementGeometry &fvGeometry,
-                        const int scvIdx) const
+                        const SubControlVolume& scv,
+                        const ElementSolutionVector& elemSol) const
     {
         return 2650; // density of sand [kg/m^3]
     }
@@ -257,8 +248,8 @@ public:
      *                    the heat capacity needs to be defined
      */
     Scalar solidThermalConductivity(const Element &element,
-                                    const FVElementGeometry &fvGeometry,
-                                    const int scvIdx) const
+                                    const SubControlVolume& scv,
+                                    const ElementSolutionVector& elemSol) const
     {
         return lambdaSolid_;
     }