diff --git a/dumux/material/spatialparams/sequentialfv.hh b/dumux/material/spatialparams/sequentialfv.hh
index d491465768e7b5eb21571afec331a747a15a9a1e..4dedf42b25ee408d3a8e40af21634372bcc84a39 100644
--- a/dumux/material/spatialparams/sequentialfv.hh
+++ b/dumux/material/spatialparams/sequentialfv.hh
@@ -53,9 +53,6 @@ class SequentialFVSpatialParams: public SequentialFVSpatialParamsOneP<TypeTag>
 
     using Element = typename GridView::template Codim<0>::Entity;
     using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
-    /// @cond false
-    using MaterialLawParams = typename GetPropType<TypeTag, Properties::MaterialLaw>::Params;
-    /// @endcond
 
 public:
     SequentialFVSpatialParams(const Problem& problem)
@@ -63,15 +60,26 @@ public:
     {
     }
 
+
+    // make sure we get a deprecation warning (remove this after release 3.3)
+    template<class ElementSolution, class SubControlVolume>
+    [[deprecated("Use the new style material laws. Old material laws and this interface will no longer be supported after release 3.3")]]
+    decltype(auto) materialLawParamsDeprecated(const Element& element,
+                                               const SubControlVolume& scv,
+                                               const ElementSolution& elemSol) const
+    {
+        return this->asImp_().materialLawParams(element);
+    }
+
     /*!
      * \brief Function for defining the parameters needed by constitutive relationships (kr-sw, pc-sw, etc.).
      *
      * \return the material parameters object
      * \param element The element
      */
-    const MaterialLawParams& materialLawParams(const Element &element) const
+    const auto& materialLawParams(const Element &element) const
     {
-            return asImp_().materialLawParamsAtPos(element.geometry().center());
+        return asImp_().materialLawParamsAtPos(element.geometry().center());
     }
 
     /*!
@@ -80,7 +88,7 @@ public:
      * \return the material parameters object
      * \param globalPos The position of the center of the element
      */
-    const MaterialLawParams& materialLawParamsAtPos(const GlobalPosition& globalPos) const
+    const auto& materialLawParamsAtPos(const GlobalPosition& globalPos) const
     {
         DUNE_THROW(Dune::InvalidStateException,
                    "The spatial parameters do not provide "
diff --git a/dumux/porousmediumflow/2p/sequential/diffusion/cellcentered/pressure.hh b/dumux/porousmediumflow/2p/sequential/diffusion/cellcentered/pressure.hh
index b3c867d583c0c7d89d84a99be1e052885c4360fc..d68266be0dd9f4f526314c3563c1c2179a337f93 100644
--- a/dumux/porousmediumflow/2p/sequential/diffusion/cellcentered/pressure.hh
+++ b/dumux/porousmediumflow/2p/sequential/diffusion/cellcentered/pressure.hh
@@ -30,6 +30,8 @@
 #include <dumux/porousmediumflow/sequential/cellcentered/pressure.hh>
 #include <dumux/porousmediumflow/2p/sequential/diffusion/properties.hh>
 
+#include <dumux/common/deprecated.hh>
+
 namespace Dumux {
 /*!
  * \ingroup SequentialTwoPModel
@@ -121,7 +123,6 @@ template<class TypeTag> class FVPressure2P: public FVPressure<TypeTag>
     using Problem = GetPropType<TypeTag, Properties::Problem>;
 
     using SpatialParams = GetPropType<TypeTag, Properties::SpatialParams>;
-    using MaterialLaw = typename SpatialParams::MaterialLaw;
 
     using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
 
@@ -920,13 +921,20 @@ const Intersection& intersection, const CellData& cellData, const bool first)
             satW = cellData.saturation(wPhaseIdx);
             satNw = cellData.saturation(nPhaseIdx);
         }
-        Scalar temperature = problem_.temperature(element);
+
+        const Scalar temperature = problem_.temperature(element);
 
         // get Dirichlet pressure boundary condition
-        Scalar pressBound = boundValues[pressureIdx];
+        const Scalar pressBound = boundValues[pressureIdx];
 
         //calculate constitutive relations depending on the kind of saturation used
-        Scalar pcBound = MaterialLaw::pc(problem_.spatialParams().materialLawParams(element), satW);
+
+        // old material law interface is deprecated: Replace this by
+        // const auto& fluidMatrixInteraction = spatialParams.fluidMatrixInteractionAtPos(element.geometry().center());
+        // after the release of 3.3, when the deprecated interface is no longer supported
+        const auto fluidMatrixInteraction = Deprecated::makePcKrSw(Scalar{}, problem_.spatialParams(), element);
+
+        const Scalar pcBound = fluidMatrixInteraction.pc(satW);
 
         // determine phase pressures from primary pressure variable
         Scalar pressW = 0;
@@ -967,9 +975,9 @@ const Intersection& intersection, const CellData& cellData, const bool first)
             rhoMeanNw = 0.5 * (cellData.density(nPhaseIdx) + densityNwBound);
         }
 
-        Scalar lambdaWBound = MaterialLaw::krw(problem_.spatialParams().materialLawParams(element), satW)
+        Scalar lambdaWBound = fluidMatrixInteraction.krw(satW)
                 / viscosityWBound;
-        Scalar lambdaNwBound = MaterialLaw::krn(problem_.spatialParams().materialLawParams(element), satW)
+        Scalar lambdaNwBound = fluidMatrixInteraction.krn(satW)
                 / viscosityNwBound;
 
         Scalar fractionalWBound = lambdaWBound / (lambdaWBound + lambdaNwBound);
@@ -1095,13 +1103,17 @@ void FVPressure2P<TypeTag>::updateMaterialLaws()
 
         CellData& cellData = problem_.variables().cellData(eIdxGlobal);
 
-        Scalar temperature = problem_.temperature(element);
+        const Scalar temperature = problem_.temperature(element);
 
         // determine phase saturations from primary saturation variable
+        const Scalar satW = cellData.saturation(wPhaseIdx);
 
-        Scalar satW = cellData.saturation(wPhaseIdx);
+        // old material law interface is deprecated: Replace this by
+        // const auto& fluidMatrixInteraction = spatialParams.fluidMatrixInteractionAtPos(element.geometry().center());
+        // after the release of 3.3, when the deprecated interface is no longer supported
+        const auto fluidMatrixInteraction = Deprecated::makePcKrSw(Scalar{}, problem_.spatialParams(), element);
 
-        Scalar pc = MaterialLaw::pc(problem_.spatialParams().materialLawParams(element), satW);
+        const Scalar pc = fluidMatrixInteraction.pc(satW);
 
         // determine phase pressures from primary pressure variable
         Scalar pressW = 0;
@@ -1151,8 +1163,8 @@ void FVPressure2P<TypeTag>::updateMaterialLaws()
         }
 
         // initialize mobilities
-        Scalar mobilityW = MaterialLaw::krw(problem_.spatialParams().materialLawParams(element), satW) / viscosity_[wPhaseIdx];
-        Scalar mobilityNw = MaterialLaw::krn(problem_.spatialParams().materialLawParams(element), satW) / viscosity_[nPhaseIdx];
+        Scalar mobilityW = fluidMatrixInteraction.krw(satW) / viscosity_[wPhaseIdx];
+        Scalar mobilityNw = fluidMatrixInteraction.krn(satW) / viscosity_[nPhaseIdx];
 
         if (compressibility_)
         {
@@ -1168,7 +1180,7 @@ void FVPressure2P<TypeTag>::updateMaterialLaws()
         cellData.setFracFlowFunc(wPhaseIdx, mobilityW / (mobilityW + mobilityNw));
         cellData.setFracFlowFunc(nPhaseIdx, mobilityNw / (mobilityW + mobilityNw));
 
-        Scalar gravityDiff = (problem_.bBoxMax() - element.geometry().center()) * gravity_;
+        const Scalar gravityDiff = (problem_.bBoxMax() - element.geometry().center()) * gravity_;
 
         Scalar potW = pressW + gravityDiff * density_[wPhaseIdx];
         Scalar potNw = pressNw + gravityDiff * density_[nPhaseIdx];
diff --git a/dumux/porousmediumflow/2p/sequential/diffusion/cellcentered/velocity.hh b/dumux/porousmediumflow/2p/sequential/diffusion/cellcentered/velocity.hh
index 7bd5b3f8a91f16168cdb849d798b47d1e6d6e00b..337358c6c6f33278bf7ff94b098bf2f28bcf0234 100644
--- a/dumux/porousmediumflow/2p/sequential/diffusion/cellcentered/velocity.hh
+++ b/dumux/porousmediumflow/2p/sequential/diffusion/cellcentered/velocity.hh
@@ -28,6 +28,8 @@
 #include <dune/grid/common/gridenums.hh>
 #include <dumux/porousmediumflow/2p/sequential/diffusion/properties.hh>
 
+#include <dumux/common/deprecated.hh>
+
 namespace Dumux {
 
 /*!
@@ -64,7 +66,6 @@ class FVVelocity2P
     using Problem = GetPropType<TypeTag, Properties::Problem>;
 
     using SpatialParams = GetPropType<TypeTag, Properties::SpatialParams>;
-    using MaterialLaw = typename SpatialParams::MaterialLaw;
 
     using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
 
@@ -564,8 +565,14 @@ void FVVelocity2P<TypeTag>::calculateVelocityOnBoundary(const Intersection& inte
             satNw = cellData.saturation(nPhaseIdx);
         }
 
-        Scalar pressBound = boundValues[pressureIdx];
-        Scalar pcBound = MaterialLaw::pc(problem_.spatialParams().materialLawParams(element), satW);
+        const Scalar pressBound = boundValues[pressureIdx];
+
+        // old material law interface is deprecated: Replace this by
+        // const auto& fluidMatrixInteraction = spatialParams.fluidMatrixInteractionAtPos(element.geometry().center());
+        // after the release of 3.3, when the deprecated interface is no longer supported
+        const auto fluidMatrixInteraction = Deprecated::makePcKrSw(Scalar{}, problem_.spatialParams(), element);
+
+        const Scalar pcBound = fluidMatrixInteraction.pc(satW);
 
         // determine phase pressures from primary pressure variable
         Scalar pressWBound = 0;
@@ -582,7 +589,7 @@ void FVVelocity2P<TypeTag>::calculateVelocityOnBoundary(const Intersection& inte
         }
 
         // get temperature at current position
-        Scalar temperature = problem_.temperature(element);
+        const Scalar temperature = problem_.temperature(element);
 
         Scalar densityWBound = density_[wPhaseIdx];
         Scalar densityNwBound = density_[nPhaseIdx];
@@ -604,9 +611,9 @@ void FVVelocity2P<TypeTag>::calculateVelocityOnBoundary(const Intersection& inte
             viscosityNwBound = FluidSystem::viscosity(fluidState, nPhaseIdx) / densityNwBound;
         }
 
-        Scalar lambdaWBound = MaterialLaw::krw(problem_.spatialParams().materialLawParams(element), satW)
+        Scalar lambdaWBound = fluidMatrixInteraction.krw(satW)
                 / viscosityWBound;
-        Scalar lambdaNwBound = MaterialLaw::krn(problem_.spatialParams().materialLawParams(element), satW)
+        Scalar lambdaNwBound = fluidMatrixInteraction.krn(satW)
                 / viscosityNwBound;
 
         Scalar potentialDiffW = cellData.fluxData().upwindPotential(wPhaseIdx, isIndex);
@@ -657,7 +664,7 @@ void FVVelocity2P<TypeTag>::calculateVelocityOnBoundary(const Intersection& inte
                     (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffNw, 0.0, 1.0e-30)) ? 0.5 * (cellData.density(nPhaseIdx) + densityNwBound) : density_[nPhaseIdx];
         }
 
-        Scalar scalarPerm = permeability.two_norm();
+        const Scalar scalarPerm = permeability.two_norm();
 
         // calculate the gravity term
         Dune::FieldVector<Scalar, dimWorld> velocityW(unitOuterNormal);
diff --git a/dumux/porousmediumflow/2p/sequential/diffusion/mimetic/mimetic.hh b/dumux/porousmediumflow/2p/sequential/diffusion/mimetic/mimetic.hh
index 4056d8f7037d01fdc636e8861de5cb6f1de78edf..c457759256f595f77e1e09749bd0a8cabe030d70 100644
--- a/dumux/porousmediumflow/2p/sequential/diffusion/mimetic/mimetic.hh
+++ b/dumux/porousmediumflow/2p/sequential/diffusion/mimetic/mimetic.hh
@@ -42,6 +42,8 @@
 
 #include <dune/common/dynvector.hh>
 
+#include <dumux/common/deprecated.hh>
+
 namespace Dumux {
 
 /*!
@@ -96,7 +98,6 @@ class MimeticTwoPLocalStiffness: public LocalStiffness<TypeTag, 1>
     using CellData = GetPropType<TypeTag, Properties::CellData>;
     using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
     using FluidState = GetPropType<TypeTag, Properties::FluidState>;
-    using MaterialLaw = GetPropType<TypeTag, Properties::MaterialLaw>;
 
 public:
     // define the number of components of your system, this is used outside
@@ -682,10 +683,13 @@ void MimeticTwoPLocalStiffness<TypeTag>::assembleElementMatrices(const Element&
                     PrimaryVariables boundValues(0.0);
                     problem_.dirichlet(boundValues, intersection);
 
-                    Scalar krw = MaterialLaw::krw(problem_.spatialParams().materialLawParams(element),
-                            boundValues[saturationIdx]);
-                    Scalar krn = MaterialLaw::krn(problem_.spatialParams().materialLawParams(element),
-                            boundValues[saturationIdx]);
+                    // old material law interface is deprecated: Replace this by
+                    // const auto& fluidMatrixInteraction = spatialParams.fluidMatrixInteractionAtPos(element.geometry().center());
+                    // after the release of 3.3, when the deprecated interface is no longer supported
+                    const auto fluidMatrixInteraction = Deprecated::makePcKrSw(Scalar{}, problem_.spatialParams(), element);
+
+                    const Scalar krw = fluidMatrixInteraction.krw(boundValues[saturationIdx]);
+                    const Scalar krn = fluidMatrixInteraction.krn(boundValues[saturationIdx]);
 
                     switch (pressureType)
                     {
diff --git a/dumux/porousmediumflow/2p/sequential/diffusion/mimetic/mimeticadaptive.hh b/dumux/porousmediumflow/2p/sequential/diffusion/mimetic/mimeticadaptive.hh
index 60585c6ccc2ff05c084c6718f1c64e829d5d1bb4..3de77fd91c467ff7daa0b5336997fbe5f46d15d7 100644
--- a/dumux/porousmediumflow/2p/sequential/diffusion/mimetic/mimeticadaptive.hh
+++ b/dumux/porousmediumflow/2p/sequential/diffusion/mimetic/mimeticadaptive.hh
@@ -44,6 +44,8 @@
 #include <dune/common/dynvector.hh>
 #include <dune/common/dynmatrix.hh>
 
+#include <dumux/common/deprecated.hh>
+
 namespace Dumux {
 
 /*!
@@ -98,7 +100,6 @@ class MimeticTwoPLocalStiffnessAdaptive: public LocalStiffness<TypeTag, 1>
     using CellData = GetPropType<TypeTag, Properties::CellData>;
     using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
     using FluidState = GetPropType<TypeTag, Properties::FluidState>;
-    using MaterialLaw = GetPropType<TypeTag, Properties::MaterialLaw>;
 
     using IntersectionMapper = Dumux::IntersectionMapper<GridView>;
 
@@ -701,10 +702,13 @@ void MimeticTwoPLocalStiffnessAdaptive<TypeTag>::assembleElementMatrices(const E
                 PrimaryVariables boundValues(0.0);
                 problem_.dirichlet(boundValues, intersection);
 
-                Scalar krw = MaterialLaw::krw(problem_.spatialParams().materialLawParams(element),
-                        boundValues[saturationIdx]);
-                Scalar krn = MaterialLaw::krn(problem_.spatialParams().materialLawParams(element),
-                        boundValues[saturationIdx]);
+                // old material law interface is deprecated: Replace this by
+                // const auto& fluidMatrixInteraction = spatialParams.fluidMatrixInteractionAtPos(element.geometry().center());
+                // after the release of 3.3, when the deprecated interface is no longer supported
+                const auto fluidMatrixInteraction = Deprecated::makePcKrSw(Scalar{}, problem_.spatialParams(), element);
+
+                const Scalar krw = fluidMatrixInteraction.krw(boundValues[saturationIdx]);
+                const Scalar krn = fluidMatrixInteraction.krn(boundValues[saturationIdx]);
 
                 switch (pressureType)
                 {
diff --git a/dumux/porousmediumflow/2p/sequential/diffusion/mimetic/operator.hh b/dumux/porousmediumflow/2p/sequential/diffusion/mimetic/operator.hh
index ca84c4a800ac7d295890a0b0a7fcdc0d53bd6804..d4d10fbc0abd087601c6f664a809394d4d3a42e1 100644
--- a/dumux/porousmediumflow/2p/sequential/diffusion/mimetic/operator.hh
+++ b/dumux/porousmediumflow/2p/sequential/diffusion/mimetic/operator.hh
@@ -28,6 +28,8 @@
 #include <dumux/porousmediumflow/2p/sequential/diffusion/properties.hh>
 #include <dumux/porousmediumflow/sequential/mimetic/properties.hh>
 
+#include <dumux/common/deprecated.hh>
+
 namespace Dumux {
 
 /*!
@@ -57,7 +59,6 @@ class MimeticOperatorAssemblerTwoP: public CROperatorAssemblerTwoP<TypeTag>
     using Element = typename GridView::template Codim<0>::Entity;
 
     using CellData = GetPropType<TypeTag, Properties::CellData>;
-    using MaterialLaw = GetPropType<TypeTag, Properties::MaterialLaw>;
     using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
     using FluidState = GetPropType<TypeTag, Properties::FluidState>;
     using BoundaryTypes = GetPropType<TypeTag, Properties::SequentialBoundaryTypes>;
@@ -247,14 +248,18 @@ public:
                         PrimaryVariables boundValues(0.0);
                         problem.dirichlet(boundValues, intersection);
 
+                        // old material law interface is deprecated: Replace this by
+                        // const auto& fluidMatrixInteraction = spatialParams.fluidMatrixInteractionAtPos(element.geometry().center());
+                        // after the release of 3.3, when the deprecated interface is no longer supported
+                        const auto fluidMatrixInteraction = Deprecated::makePcKrSw(Scalar{}, problem.spatialParams(), element);
+
                         if (velocityW[idxInInside] >= 0.)
                         {
                             mobilityW = cellData.mobility(wPhaseIdx);
                         }
                         else
                         {
-                            mobilityW = MaterialLaw::krw(problem.spatialParams().materialLawParams(element),
-                                    boundValues[saturationIdx]) / viscosityW;
+                            mobilityW = fluidMatrixInteraction.krw(boundValues[saturationIdx]) / viscosityW;
                         }
 
                         if (velocityNw[idxInInside] >= 0.)
@@ -263,8 +268,7 @@ public:
                         }
                         else
                         {
-                            mobilityNw = MaterialLaw::krn(problem.spatialParams().materialLawParams(element),
-                                    boundValues[saturationIdx]) / viscosityNw;
+                            mobilityNw = fluidMatrixInteraction.krn(boundValues[saturationIdx]) / viscosityNw;
                         }
                     }
                     else
diff --git a/dumux/porousmediumflow/2p/sequential/diffusion/mimetic/operatoradaptive.hh b/dumux/porousmediumflow/2p/sequential/diffusion/mimetic/operatoradaptive.hh
index d204e8af45dc7dc305dac41f949f733f8e5e715d..32135a2b63be61cdcd051e56b6a0fadb10e038bd 100644
--- a/dumux/porousmediumflow/2p/sequential/diffusion/mimetic/operatoradaptive.hh
+++ b/dumux/porousmediumflow/2p/sequential/diffusion/mimetic/operatoradaptive.hh
@@ -28,6 +28,8 @@
 #include <dumux/porousmediumflow/2p/sequential/diffusion/properties.hh>
 #include <dumux/porousmediumflow/sequential/mimetic/properties.hh>
 
+#include <dumux/common/deprecated.hh>
+
 namespace Dumux {
 
 /*!
@@ -59,7 +61,6 @@ class MimeticOperatorAssemblerTwoPAdaptive : public CROperatorAssemblerTwoPAdapt
     using Element = typename GridView::template Codim<0>::Entity;
 
     using CellData = GetPropType<TypeTag, Properties::CellData>;
-    using MaterialLaw = GetPropType<TypeTag, Properties::MaterialLaw>;
     using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
     using FluidState = GetPropType<TypeTag, Properties::FluidState>;
     using BoundaryTypes = GetPropType<TypeTag, Properties::SequentialBoundaryTypes>;
@@ -255,14 +256,18 @@ public:
                         PrimaryVariables boundValues(0.0);
                         problem.dirichlet(boundValues, intersection);
 
+                        // old material law interface is deprecated: Replace this by
+                        // const auto& fluidMatrixInteraction = spatialParams.fluidMatrixInteractionAtPos(element.geometry().center());
+                        // after the release of 3.3, when the deprecated interface is no longer supported
+                        const auto fluidMatrixInteraction = Deprecated::makePcKrSw(Scalar{}, problem.spatialParams(), element);
+
                         if (velocityW[intersectionIdx] >= 0.)
                         {
                             mobilityW = cellData.mobility(wPhaseIdx);
                         }
                         else
                         {
-                            mobilityW = MaterialLaw::krw(problem.spatialParams().materialLawParams(element),
-                                boundValues[saturationIdx]) / viscosityW;
+                            mobilityW = fluidMatrixInteraction.krw(boundValues[saturationIdx]) / viscosityW;
                         }
 
                         if (velocityNw[intersectionIdx] >= 0.)
@@ -271,8 +276,7 @@ public:
                         }
                         else
                         {
-                            mobilityNw = MaterialLaw::krn(problem.spatialParams().materialLawParams(element),
-                                boundValues[saturationIdx]) / viscosityNw;
+                            mobilityNw = fluidMatrixInteraction.krn(boundValues[saturationIdx]) / viscosityNw;
                         }
                     }
                     else
diff --git a/dumux/porousmediumflow/2p/sequential/diffusion/mimetic/pressure.hh b/dumux/porousmediumflow/2p/sequential/diffusion/mimetic/pressure.hh
index 92fa2c1259aa26fb13c6bbb4b0d7b746482aa480..9da440fb954c4c52bbb8503da52d7454819fb1d2 100644
--- a/dumux/porousmediumflow/2p/sequential/diffusion/mimetic/pressure.hh
+++ b/dumux/porousmediumflow/2p/sequential/diffusion/mimetic/pressure.hh
@@ -32,6 +32,8 @@
 #include <dumux/porousmediumflow/2p/sequential/diffusion/mimetic/operator.hh>
 #include <dumux/porousmediumflow/2p/sequential/diffusion/mimetic/mimetic.hh>
 
+#include <dumux/common/deprecated.hh>
+
 namespace Dumux {
 
 /*!
@@ -67,7 +69,6 @@ template<class TypeTag> class MimeticPressure2P
     using Problem = GetPropType<TypeTag, Properties::Problem>;
 
     using SpatialParams = GetPropType<TypeTag, Properties::SpatialParams>;
-    using MaterialLaw = typename SpatialParams::MaterialLaw;
 
     using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
 
@@ -515,13 +516,16 @@ void MimeticPressure2P<TypeTag>::updateMaterialLaws()
 
             CellData& cellData = problem_.variables().cellData(eIdxGlobal);
 
-            Scalar satW = cellData.saturation(wPhaseIdx);
+            const Scalar satW = cellData.saturation(wPhaseIdx);
+
+            // old material law interface is deprecated: Replace this by
+            // const auto& fluidMatrixInteraction = spatialParams.fluidMatrixInteractionAtPos(element.geometry().center());
+            // after the release of 3.3, when the deprecated interface is no longer supported
+            const auto fluidMatrixInteraction = Deprecated::makePcKrSw(Scalar{}, problem_.spatialParams(), element);
 
             // initialize mobilities
-            Scalar mobilityW = MaterialLaw::krw(problem_.spatialParams().materialLawParams(element), satW)
-                    / viscosity_[wPhaseIdx];
-            Scalar mobilityNw = MaterialLaw::krn(problem_.spatialParams().materialLawParams(element), satW)
-                    / viscosity_[nPhaseIdx];
+            const Scalar mobilityW = fluidMatrixInteraction.krw(satW) / viscosity_[wPhaseIdx];
+            const Scalar mobilityNw = fluidMatrixInteraction.krn(satW) / viscosity_[nPhaseIdx];
 
             // initialize mobilities
             cellData.setMobility(wPhaseIdx, mobilityW);
diff --git a/dumux/porousmediumflow/2p/sequential/diffusion/mimetic/pressureadaptive.hh b/dumux/porousmediumflow/2p/sequential/diffusion/mimetic/pressureadaptive.hh
index fb0f271f81bc7401adb27cfe5d3feba591b44013..7ecf0b88390f2656af98dc27e73fa2d70d9ee530 100644
--- a/dumux/porousmediumflow/2p/sequential/diffusion/mimetic/pressureadaptive.hh
+++ b/dumux/porousmediumflow/2p/sequential/diffusion/mimetic/pressureadaptive.hh
@@ -32,6 +32,8 @@
 #include <dumux/porousmediumflow/2p/sequential/diffusion/mimetic/operatoradaptive.hh>
 #include <dumux/porousmediumflow/2p/sequential/diffusion/mimetic/mimeticadaptive.hh>
 
+#include <dumux/common/deprecated.hh>
+
 namespace Dumux {
 
 /*!
@@ -67,7 +69,6 @@ template<class TypeTag> class MimeticPressure2PAdaptive
     using Problem = GetPropType<TypeTag, Properties::Problem>;
 
     using SpatialParams = GetPropType<TypeTag, Properties::SpatialParams>;
-    using MaterialLaw = typename SpatialParams::MaterialLaw;
 
     using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
 
@@ -532,13 +533,16 @@ void MimeticPressure2PAdaptive<TypeTag>::updateMaterialLaws()
 
         CellData& cellData = problem_.variables().cellData(eIdxGlobal);
 
-        Scalar satW = cellData.saturation(wPhaseIdx);
+        const Scalar satW = cellData.saturation(wPhaseIdx);
+
+        // old material law interface is deprecated: Replace this by
+        // const auto& fluidMatrixInteraction = spatialParams.fluidMatrixInteractionAtPos(element.geometry().center());
+        // after the release of 3.3, when the deprecated interface is no longer supported
+        const auto fluidMatrixInteraction = Deprecated::makePcKrSw(Scalar{}, problem_.spatialParams(), element);
 
         // initialize mobilities
-        Scalar mobilityW = MaterialLaw::krw(problem_.spatialParams().materialLawParams(element), satW)
-                / viscosity_[wPhaseIdx];
-        Scalar mobilityNw = MaterialLaw::krn(problem_.spatialParams().materialLawParams(element), satW)
-                / viscosity_[nPhaseIdx];
+        const Scalar mobilityW = fluidMatrixInteraction.krw(satW) / viscosity_[wPhaseIdx];
+        const Scalar mobilityNw = fluidMatrixInteraction.krn(satW) / viscosity_[nPhaseIdx];
 
         // initialize mobilities
         cellData.setMobility(wPhaseIdx, mobilityW);
diff --git a/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/2dpressure.hh b/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/2dpressure.hh
index 34c6c1fabbee46583ad257bc8c2853498cceb422..d5491eaa97202573cbe0c36b69b581eeeceaf426 100644
--- a/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/2dpressure.hh
+++ b/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/2dpressure.hh
@@ -29,6 +29,8 @@
 #include <dumux/porousmediumflow/sequential/cellcentered/mpfa/properties.hh>
 #include "2dtransmissibilitycalculator.hh"
 
+#include <dumux/common/deprecated.hh>
+
 namespace Dumux {
 
 /*!
@@ -82,7 +84,6 @@ class FvMpfaL2dPressure2p: public FVPressure<TypeTag>
     using Problem = GetPropType<TypeTag, Properties::Problem>;
 
     using SpatialParams = GetPropType<TypeTag, Properties::SpatialParams>;
-    using MaterialLaw = typename SpatialParams::MaterialLaw;
 
     using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
 
@@ -1664,19 +1665,20 @@ void FvMpfaL2dPressure2p<TypeTag>::assemble()
 
                             }
 
-                            Scalar pcBound = MaterialLaw::pc(
-                                    problem_.spatialParams().materialLawParams(element), satWBound);
+                            // old material law interface is deprecated: Replace this by
+                            // const auto& fluidMatrixInteraction = spatialParams.fluidMatrixInteractionAtPos(element.geometry().center());
+                            // after the release of 3.3, when the deprecated interface is no longer supported
+                            const auto fluidMatrixInteraction = Deprecated::makePcKrSw(Scalar{}, problem_.spatialParams(), element);
+
+                            Scalar pcBound = fluidMatrixInteraction.pc(satWBound);
 
                             Scalar gravityDiffBound = (problem_.bBoxMax() - globalPosFace) * gravity_
                                     * (density_[nPhaseIdx] - density_[wPhaseIdx]);
 
                             pcBound += gravityDiffBound;
 
-                            Dune::FieldVector<Scalar, numPhases> lambdaBound(
-                                    MaterialLaw::krw(problem_.spatialParams().materialLawParams(element),
-                                            satWBound));
-                            lambdaBound[nPhaseIdx] = MaterialLaw::krn(
-                                    problem_.spatialParams().materialLawParams(element), satWBound);
+                            Dune::FieldVector<Scalar, numPhases> lambdaBound(fluidMatrixInteraction.krw(satWBound));
+                            lambdaBound[nPhaseIdx] = fluidMatrixInteraction.krn(satWBound);
                             lambdaBound[wPhaseIdx] /= viscosity_[wPhaseIdx];
                             lambdaBound[nPhaseIdx] /= viscosity_[nPhaseIdx];
 
@@ -1840,15 +1842,18 @@ void FvMpfaL2dPressure2p<TypeTag>::updateMaterialLaws()
 
         Scalar satW = cellData.saturation(wPhaseIdx);
 
-        Scalar pc = MaterialLaw::pc(problem_.spatialParams().materialLawParams(element), satW);
+        // old material law interface is deprecated: Replace this by
+        // const auto& fluidMatrixInteraction = spatialParams.fluidMatrixInteractionAtPos(element.geometry().center());
+        // after the release of 3.3, when the deprecated interface is no longer supported
+        const auto fluidMatrixInteraction = Deprecated::makePcKrSw(Scalar{}, problem_.spatialParams(), element);
+
+        const Scalar pc = fluidMatrixInteraction.pc(satW);
 
         cellData.setCapillaryPressure(pc);
 
         // initialize mobilities
-        Scalar mobilityW = MaterialLaw::krw(problem_.spatialParams().materialLawParams(element), satW)
-                / viscosity_[wPhaseIdx];
-        Scalar mobilityNw = MaterialLaw::krn(problem_.spatialParams().materialLawParams(element), satW)
-                / viscosity_[nPhaseIdx];
+        const Scalar mobilityW = fluidMatrixInteraction.krw(satW) / viscosity_[wPhaseIdx];
+        const Scalar mobilityNw = fluidMatrixInteraction.krn(satW) / viscosity_[nPhaseIdx];
 
         // initialize mobilities
         cellData.setMobility(wPhaseIdx, mobilityW);
diff --git a/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/2dpressureadaptive.hh b/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/2dpressureadaptive.hh
index e19b8d174061157570608baf227fc08bf82b694a..6813decae3097f5d93591889f15202b846a1804a 100644
--- a/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/2dpressureadaptive.hh
+++ b/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/2dpressureadaptive.hh
@@ -30,6 +30,8 @@
 #include <dumux/porousmediumflow/sequential/cellcentered/mpfa/properties.hh>
 #include "2dtransmissibilitycalculator.hh"
 
+#include <dumux/common/deprecated.hh>
+
 namespace Dumux {
 
 /*!
@@ -85,7 +87,6 @@ class FvMpfaL2dPressure2pAdaptive: public FVPressure<TypeTag>
     using Problem = GetPropType<TypeTag, Properties::Problem>;
 
     using SpatialParams = GetPropType<TypeTag, Properties::SpatialParams>;
-    using MaterialLaw = typename SpatialParams::MaterialLaw;
 
     using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
 
@@ -2383,19 +2384,20 @@ void FvMpfaL2dPressure2pAdaptive<TypeTag>::assemble()
 
                             }
 
-                            Scalar pcBound = MaterialLaw::pc(
-                                    problem_.spatialParams().materialLawParams(element), satWBound);
+                            // old material law interface is deprecated: Replace this by
+                            // const auto& fluidMatrixInteraction = spatialParams.fluidMatrixInteractionAtPos(element.geometry().center());
+                            // after the release of 3.3, when the deprecated interface is no longer supported
+                            const auto fluidMatrixInteraction = Deprecated::makePcKrSw(Scalar{}, problem_.spatialParams(), element);
+
+                            Scalar pcBound = fluidMatrixInteraction.pc(satWBound);
 
                             Scalar gravityDiffBound = (problem_.bBoxMax() - globalPosFace) * gravity_
                                     * (density_[nPhaseIdx] - density_[wPhaseIdx]);
 
                             pcBound += gravityDiffBound;
 
-                            Dune::FieldVector<Scalar, numPhases> lambdaBound(
-                                    MaterialLaw::krw(problem_.spatialParams().materialLawParams(element),
-                                            satWBound));
-                            lambdaBound[nPhaseIdx] = MaterialLaw::krn(
-                                    problem_.spatialParams().materialLawParams(element), satWBound);
+                            Dune::FieldVector<Scalar, numPhases> lambdaBound(fluidMatrixInteraction.krw(satWBound));
+                            lambdaBound[nPhaseIdx] = fluidMatrixInteraction.krn(satWBound);
                             lambdaBound[wPhaseIdx] /= viscosity_[wPhaseIdx];
                             lambdaBound[nPhaseIdx] /= viscosity_[nPhaseIdx];
 
@@ -2556,17 +2558,20 @@ void FvMpfaL2dPressure2pAdaptive<TypeTag>::updateMaterialLaws()
 
         CellData& cellData = problem_.variables().cellData(eIdxGlobal);
 
-        Scalar satW = cellData.saturation(wPhaseIdx);
+        const Scalar satW = cellData.saturation(wPhaseIdx);
+
+        // old material law interface is deprecated: Replace this by
+        // const auto& fluidMatrixInteraction = spatialParams.fluidMatrixInteractionAtPos(element.geometry().center());
+        // after the release of 3.3, when the deprecated interface is no longer supported
+        const auto fluidMatrixInteraction = Deprecated::makePcKrSw(Scalar{}, problem_.spatialParams(), element);
 
-        Scalar pc = MaterialLaw::pc(problem_.spatialParams().materialLawParams(element), satW);
+        const Scalar pc = fluidMatrixInteraction.pc(satW);
 
         cellData.setCapillaryPressure(pc);
 
         // initialize mobilities
-        Scalar mobilityW = MaterialLaw::krw(problem_.spatialParams().materialLawParams(element), satW)
-                / viscosity_[wPhaseIdx];
-        Scalar mobilityNw = MaterialLaw::krn(problem_.spatialParams().materialLawParams(element), satW)
-                / viscosity_[nPhaseIdx];
+        const Scalar mobilityW = fluidMatrixInteraction.krw(satW) / viscosity_[wPhaseIdx];
+        const Scalar mobilityNw = fluidMatrixInteraction.krn(satW) / viscosity_[nPhaseIdx];
 
         // initialize mobilities
         cellData.setMobility(wPhaseIdx, mobilityW);
diff --git a/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/2dpressurevelocity.hh b/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/2dpressurevelocity.hh
index 59b40398247266434f7e834a16ecd959c0638fcc..15b59fadef74300305ddf1e8ee2ef0a74a021d51 100644
--- a/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/2dpressurevelocity.hh
+++ b/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/2dpressurevelocity.hh
@@ -27,6 +27,8 @@
 #include "2dpressure.hh"
 #include "2dvelocity.hh"
 
+#include <dumux/common/deprecated.hh>
+
 namespace Dumux {
 
 /*!
@@ -61,7 +63,6 @@ template<class TypeTag> class FvMpfaL2dPressureVelocity2p: public FvMpfaL2dPress
     using PrimaryVariables = typename SolutionTypes::PrimaryVariables;
 
     using SpatialParams = GetPropType<TypeTag, Properties::SpatialParams>;
-    using MaterialLaw = typename SpatialParams::MaterialLaw;
 
     using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
 
@@ -453,8 +454,14 @@ void FvMpfaL2dPressureVelocity2p<TypeTag>::calculateVelocityOnBoundary(const Int
             satW = cellData.saturation(wPhaseIdx);
         }
 
-        Scalar pressBound = boundValues[pressureIdx];
-        Scalar pcBound = MaterialLaw::pc(problem_.spatialParams().materialLawParams(element), satW);
+        const Scalar pressBound = boundValues[pressureIdx];
+
+        // old material law interface is deprecated: Replace this by
+        // const auto& fluidMatrixInteraction = spatialParams.fluidMatrixInteractionAtPos(element.geometry().center());
+        // after the release of 3.3, when the deprecated interface is no longer supported
+        const auto fluidMatrixInteraction = Deprecated::makePcKrSw(Scalar{}, problem_.spatialParams(), element);
+
+        const Scalar pcBound = fluidMatrixInteraction.pc(satW);
 
         //determine phase pressures from primary pressure variable
         Scalar pressWBound = 0;
@@ -470,10 +477,8 @@ void FvMpfaL2dPressureVelocity2p<TypeTag>::calculateVelocityOnBoundary(const Int
             pressNwBound = pressBound;
         }
 
-        Scalar lambdaWBound = MaterialLaw::krw(problem_.spatialParams().materialLawParams(element), satW)
-                / viscosity_[wPhaseIdx];
-        Scalar lambdaNwBound = MaterialLaw::krn(problem_.spatialParams().materialLawParams(element), satW)
-                / viscosity_[nPhaseIdx];
+        const Scalar lambdaWBound = fluidMatrixInteraction.krw(satW) / viscosity_[wPhaseIdx];
+        const Scalar lambdaNwBound = fluidMatrixInteraction.krn(satW) / viscosity_[nPhaseIdx];
 
         Scalar potentialDiffW = cellData.fluxData().upwindPotential(wPhaseIdx, isIndex);
         Scalar potentialDiffNw = cellData.fluxData().upwindPotential(nPhaseIdx, isIndex);
diff --git a/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/2dpressurevelocityadaptive.hh b/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/2dpressurevelocityadaptive.hh
index d4a741869292c5ddf4ede6a4017ba094f8f2740e..87b3176ead6acf59ef8d78d999019f7bd0699b4b 100644
--- a/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/2dpressurevelocityadaptive.hh
+++ b/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/2dpressurevelocityadaptive.hh
@@ -27,6 +27,8 @@
 #include "2dpressureadaptive.hh"
 #include "2dvelocityadaptive.hh"
 
+#include <dumux/common/deprecated.hh>
+
 namespace Dumux {
 
 /*!
@@ -63,7 +65,6 @@ template<class TypeTag> class FvMpfaL2dPressureVelocity2pAdaptive: public FvMpfa
     using PrimaryVariables = typename SolutionTypes::PrimaryVariables;
 
     using SpatialParams = GetPropType<TypeTag, Properties::SpatialParams>;
-    using MaterialLaw = typename SpatialParams::MaterialLaw;
 
     using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
 
@@ -589,8 +590,14 @@ void FvMpfaL2dPressureVelocity2pAdaptive<TypeTag>::calculateVelocityOnBoundary(c
             satW = cellData.saturation(wPhaseIdx);
         }
 
-        Scalar pressBound = boundValues[pressureIdx];
-        Scalar pcBound = MaterialLaw::pc(problem_.spatialParams().materialLawParams(element), satW);
+        const Scalar pressBound = boundValues[pressureIdx];
+
+        // old material law interface is deprecated: Replace this by
+        // const auto& fluidMatrixInteraction = spatialParams.fluidMatrixInteractionAtPos(element.geometry().center());
+        // after the release of 3.3, when the deprecated interface is no longer supported
+        const auto fluidMatrixInteraction = Deprecated::makePcKrSw(Scalar{}, problem_.spatialParams(), element);
+
+        const Scalar pcBound = fluidMatrixInteraction.pc(satW);
 
         //determine phase pressures from primary pressure variable
         Scalar pressWBound = 0;
@@ -606,10 +613,8 @@ void FvMpfaL2dPressureVelocity2pAdaptive<TypeTag>::calculateVelocityOnBoundary(c
             pressNwBound = pressBound;
         }
 
-        Scalar lambdaWBound = MaterialLaw::krw(problem_.spatialParams().materialLawParams(element), satW)
-                / viscosity_[wPhaseIdx];
-        Scalar lambdaNwBound = MaterialLaw::krn(problem_.spatialParams().materialLawParams(element), satW)
-                / viscosity_[nPhaseIdx];
+        const Scalar lambdaWBound = fluidMatrixInteraction.krw(satW) / viscosity_[wPhaseIdx];
+        const Scalar lambdaNwBound = fluidMatrixInteraction.krn(satW) / viscosity_[nPhaseIdx];
 
         Scalar potentialDiffW = cellData.fluxData().upwindPotential(wPhaseIdx, isIndex);
         Scalar potentialDiffNw = cellData.fluxData().upwindPotential(nPhaseIdx, isIndex);
diff --git a/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/2dvelocity.hh b/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/2dvelocity.hh
index 269960826c4a6b674404ffa3a7f2bd0dce5e9458..0e6df329cd43c587de7c1bcfbde16aee595df24e 100644
--- a/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/2dvelocity.hh
+++ b/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/2dvelocity.hh
@@ -30,6 +30,8 @@
 #include <dumux/porousmediumflow/sequential/cellcentered/mpfa/linteractionvolume.hh>
 #include "2dtransmissibilitycalculator.hh"
 
+#include <dumux/common/deprecated.hh>
+
 namespace Dumux {
 
 /*!
@@ -66,7 +68,6 @@ template<class TypeTag> class FvMpfaL2dVelocity2p
     using Problem = GetPropType<TypeTag, Properties::Problem>;
 
     using SpatialParams = GetPropType<TypeTag, Properties::SpatialParams>;
-    using MaterialLaw = typename SpatialParams::MaterialLaw;
 
     using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
 
@@ -786,21 +787,20 @@ void FvMpfaL2dVelocity2p<TypeTag>::calculateBoundaryInteractionVolumeVelocity(In
 
                     }
 
-                    Scalar pcBound = MaterialLaw::pc(
-                            problem_.spatialParams().materialLawParams(element), satWBound);
+                    // old material law interface is deprecated: Replace this by
+                    // const auto& fluidMatrixInteraction = spatialParams.fluidMatrixInteractionAtPos(element.geometry().center());
+                    // after the release of 3.3, when the deprecated interface is no longer supported
+                    const auto fluidMatrixInteraction = Deprecated::makePcKrSw(Scalar{}, problem_.spatialParams(), element);
+
+                    Scalar pcBound = fluidMatrixInteraction.pc(satWBound);
 
                     Scalar gravityDiffBound = (problem_.bBoxMax() - globalPosFace) * gravity_
                             * (density_[nPhaseIdx] - density_[wPhaseIdx]);
 
                     pcBound += gravityDiffBound;
 
-                    Dune::FieldVector < Scalar, numPhases
-                            > lambdaBound(
-                                    MaterialLaw::krw(
-                                            problem_.spatialParams().materialLawParams(element),
-                                            satWBound));
-                    lambdaBound[nPhaseIdx] = MaterialLaw::krn(
-                            problem_.spatialParams().materialLawParams(element), satWBound);
+                    Dune::FieldVector <Scalar, numPhases> lambdaBound(fluidMatrixInteraction.krw(satWBound));
+                    lambdaBound[nPhaseIdx] = fluidMatrixInteraction.krn(satWBound);
                     lambdaBound[wPhaseIdx] /= viscosity_[wPhaseIdx];
                     lambdaBound[nPhaseIdx] /= viscosity_[nPhaseIdx];
 
diff --git a/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/2dvelocityadaptive.hh b/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/2dvelocityadaptive.hh
index 2fddf014ca55d706416d8bab54c12dd2bb8a06fb..85cbc3586bb69a5742d19fb1caf38d595068d2ac 100644
--- a/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/2dvelocityadaptive.hh
+++ b/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/2dvelocityadaptive.hh
@@ -63,7 +63,6 @@ template<class TypeTag> class FvMpfaL2dVelocity2pAdaptive : public FvMpfaL2dVelo
     using Problem = GetPropType<TypeTag, Properties::Problem>;
 
     using SpatialParams = GetPropType<TypeTag, Properties::SpatialParams>;
-    using MaterialLaw = typename SpatialParams::MaterialLaw;
 
     using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
 
diff --git a/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/3dpressure.hh b/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/3dpressure.hh
index 68a418328fb95d7ce21e5f78e3cf7f9c56e007d9..8eefe30623c7caf559eeccea2256e062271c0d69 100644
--- a/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/3dpressure.hh
+++ b/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/3dpressure.hh
@@ -32,6 +32,8 @@
 #include "3dinteractionvolumecontainer.hh"
 #include "3dtransmissibilitycalculator.hh"
 
+#include <dumux/common/deprecated.hh>
+
 namespace Dumux {
 
 /*!
@@ -85,7 +87,6 @@ class FvMpfaL3dPressure2p: public FVPressure<TypeTag>
     using Problem = GetPropType<TypeTag, Properties::Problem>;
 
     using SpatialParams = GetPropType<TypeTag, Properties::SpatialParams>;
-    using MaterialLaw = typename SpatialParams::MaterialLaw;
 
     using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
 
@@ -2351,18 +2352,20 @@ void FvMpfaL3dPressure2p<TypeTag>::assembleBoundaryInteractionVolume(Interaction
 
                     }
 
-                    Scalar pcBound = MaterialLaw::pc(
-                                                     problem_.spatialParams().materialLawParams(element), satWBound);
+                    // old material law interface is deprecated: Replace this by
+                    // const auto& fluidMatrixInteraction = spatialParams.fluidMatrixInteractionAtPos(element.geometry().center());
+                    // after the release of 3.3, when the deprecated interface is no longer supported
+                    const auto fluidMatrixInteraction = Deprecated::makePcKrSw(Scalar{}, problem_.spatialParams(), element);
+
+                    Scalar pcBound = fluidMatrixInteraction.pc(satWBound);
 
                     Scalar gravityDiffBound = (problem_.bBoxMax() - globalPosFace) * gravity_
                         * (density_[nPhaseIdx] - density_[wPhaseIdx]);
 
                     pcBound += gravityDiffBound;
 
-                    Dune::FieldVector<Scalar, numPhases>
-                      lambdaBound(MaterialLaw::krw(problem_.spatialParams().materialLawParams(element),satWBound));
-                    lambdaBound[nPhaseIdx] = MaterialLaw::krn(
-                                                              problem_.spatialParams().materialLawParams(element), satWBound);
+                    Dune::FieldVector<Scalar, numPhases> lambdaBound(fluidMatrixInteraction.krw(satWBound));
+                    lambdaBound[nPhaseIdx] = fluidMatrixInteraction.krn(satWBound);
                     lambdaBound[wPhaseIdx] /= viscosity_[wPhaseIdx];
                     lambdaBound[nPhaseIdx] /= viscosity_[nPhaseIdx];
 
@@ -2489,14 +2492,19 @@ void FvMpfaL3dPressure2p<TypeTag>::updateMaterialLaws()
 
         Scalar satW = cellData.saturation(wPhaseIdx);
 
-        Scalar pc = MaterialLaw::pc(problem_.spatialParams().materialLawParams(element), satW);
+        // old material law interface is deprecated: Replace this by
+        // const auto& fluidMatrixInteraction = spatialParams.fluidMatrixInteractionAtPos(element.geometry().center());
+        // after the release of 3.3, when the deprecated interface is no longer supported
+        const auto fluidMatrixInteraction = Deprecated::makePcKrSw(Scalar{}, problem_.spatialParams(), element);
+
+        const Scalar pc = fluidMatrixInteraction.pc(satW);
 
         cellData.setCapillaryPressure(pc);
 
         // initialize mobilities
-        Scalar mobilityW = MaterialLaw::krw(problem_.spatialParams().materialLawParams(element), satW)
+        Scalar mobilityW = fluidMatrixInteraction.krw(satW)
             / viscosity_[wPhaseIdx];
-        Scalar mobilityNw = MaterialLaw::krn(problem_.spatialParams().materialLawParams(element), satW)
+        Scalar mobilityNw = fluidMatrixInteraction.krn(satW)
             / viscosity_[nPhaseIdx];
 
         // initialize mobilities
diff --git a/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/3dpressureadaptive.hh b/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/3dpressureadaptive.hh
index a5129ed870ea096794f3e1bf2936a2b572e1677a..ef2aef01af18b077ea2085f47015db0f9901df36 100644
--- a/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/3dpressureadaptive.hh
+++ b/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/3dpressureadaptive.hh
@@ -89,7 +89,6 @@ class FvMpfaL3dPressure2pAdaptive: public FvMpfaL3dPressure2p<TypeTag>
     using Problem = GetPropType<TypeTag, Properties::Problem>;
 
     using SpatialParams = GetPropType<TypeTag, Properties::SpatialParams>;
-    using MaterialLaw = typename SpatialParams::MaterialLaw;
 
     using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
 
diff --git a/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/3dpressurevelocity.hh b/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/3dpressurevelocity.hh
index 025df282e2959caae797b596bf9dd99e97dcbfa7..595ff582d94fb2393b487115671bbbad7819f2d1 100644
--- a/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/3dpressurevelocity.hh
+++ b/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/3dpressurevelocity.hh
@@ -28,6 +28,8 @@
 #include "3dpressure.hh"
 #include "3dvelocity.hh"
 
+#include <dumux/common/deprecated.hh>
+
 namespace Dumux {
 /*!
  * \ingroup SequentialTwoPModel
@@ -77,7 +79,6 @@ template<class TypeTag> class FvMpfaL3dPressureVelocity2p: public FvMpfaL3dPress
     using SolutionTypes = GetProp<TypeTag, Properties::SolutionTypes>;
     using PrimaryVariables = typename SolutionTypes::PrimaryVariables;
     using SpatialParams = GetPropType<TypeTag, Properties::SpatialParams>;
-    using MaterialLaw = typename SpatialParams::MaterialLaw;
 
     using InteractionVolume = GetPropType<TypeTag, Properties::MPFAInteractionVolume>;
     using Intersection = typename GridView::Intersection;
@@ -452,8 +453,14 @@ void FvMpfaL3dPressureVelocity2p<TypeTag>::calculateVelocityOnBoundary(const Int
             satW = cellData.saturation(wPhaseIdx);
         }
 
-        Scalar pressBound = boundValues[pressureIdx];
-        Scalar pcBound = MaterialLaw::pc(problem_.spatialParams().materialLawParams(element), satW);
+        const Scalar pressBound = boundValues[pressureIdx];
+
+        // old material law interface is deprecated: Replace this by
+        // const auto& fluidMatrixInteraction = spatialParams.fluidMatrixInteractionAtPos(element.geometry().center());
+        // after the release of 3.3, when the deprecated interface is no longer supported
+        const auto fluidMatrixInteraction = Deprecated::makePcKrSw(Scalar{}, problem_.spatialParams(), element);
+
+        const Scalar pcBound = fluidMatrixInteraction.pc(satW);
 
         //determine phase pressures from primary pressure variable
         Scalar pressWBound = 0;
@@ -469,10 +476,8 @@ void FvMpfaL3dPressureVelocity2p<TypeTag>::calculateVelocityOnBoundary(const Int
             pressNwBound = pressBound;
         }
 
-        Scalar lambdaWBound = MaterialLaw::krw(problem_.spatialParams().materialLawParams(element), satW)
-                / viscosity_[wPhaseIdx];
-        Scalar lambdaNwBound = MaterialLaw::krn(problem_.spatialParams().materialLawParams(element), satW)
-                / viscosity_[nPhaseIdx];
+        const Scalar lambdaWBound = fluidMatrixInteraction.krw(satW) / viscosity_[wPhaseIdx];
+        const Scalar lambdaNwBound = fluidMatrixInteraction.krn(satW) / viscosity_[nPhaseIdx];
 
         Scalar potentialDiffW = cellData.fluxData().upwindPotential(wPhaseIdx, isIndex);
         Scalar potentialDiffNw = cellData.fluxData().upwindPotential(nPhaseIdx, isIndex);
diff --git a/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/3dpressurevelocityadaptive.hh b/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/3dpressurevelocityadaptive.hh
index 5782b15acf17d3b4e3bec7f2046b468e6dc51877..84a8ba579aaca865da74121371139e95e0922ec4 100644
--- a/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/3dpressurevelocityadaptive.hh
+++ b/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/3dpressurevelocityadaptive.hh
@@ -27,6 +27,8 @@
 #include "3dpressureadaptive.hh"
 #include "3dvelocityadaptive.hh"
 
+#include <dumux/common/deprecated.hh>
+
 namespace Dumux {
 
 /*!
@@ -79,7 +81,6 @@ template<class TypeTag> class FvMpfaL3dPressureVelocity2pAdaptive: public FvMpfa
     using SolutionTypes = GetProp<TypeTag, Properties::SolutionTypes>;
     using PrimaryVariables = typename SolutionTypes::PrimaryVariables;
     using SpatialParams = GetPropType<TypeTag, Properties::SpatialParams>;
-    using MaterialLaw = typename SpatialParams::MaterialLaw;
 
     using InteractionVolume = GetPropType<TypeTag, Properties::MPFAInteractionVolume>;
     using Intersection = typename GridView::Intersection;
@@ -553,8 +554,14 @@ void FvMpfaL3dPressureVelocity2pAdaptive<TypeTag>::calculateVelocityOnBoundary(c
             satW = cellData.saturation(wPhaseIdx);
         }
 
-        Scalar pressBound = boundValues[pressureIdx];
-        Scalar pcBound = MaterialLaw::pc(problem_.spatialParams().materialLawParams(element), satW);
+        const Scalar pressBound = boundValues[pressureIdx];
+
+        // old material law interface is deprecated: Replace this by
+        // const auto& fluidMatrixInteraction = spatialParams.fluidMatrixInteractionAtPos(element.geometry().center());
+        // after the release of 3.3, when the deprecated interface is no longer supported
+        const auto fluidMatrixInteraction = Deprecated::makePcKrSw(Scalar{}, problem_.spatialParams(), element);
+
+        const Scalar pcBound = fluidMatrixInteraction.pc(satW);
 
         //determine phase pressures from primary pressure variable
         Scalar pressWBound = 0;
@@ -570,10 +577,8 @@ void FvMpfaL3dPressureVelocity2pAdaptive<TypeTag>::calculateVelocityOnBoundary(c
             pressNwBound = pressBound;
         }
 
-        Scalar lambdaWBound = MaterialLaw::krw(problem_.spatialParams().materialLawParams(element), satW)
-        / viscosity_[wPhaseIdx];
-        Scalar lambdaNwBound = MaterialLaw::krn(problem_.spatialParams().materialLawParams(element), satW)
-        / viscosity_[nPhaseIdx];
+        const Scalar lambdaWBound = fluidMatrixInteraction.krw(satW) / viscosity_[wPhaseIdx];
+        const Scalar lambdaNwBound = fluidMatrixInteraction.krn(satW) / viscosity_[nPhaseIdx];
 
         Scalar potentialDiffW = cellData.fluxData().upwindPotential(wPhaseIdx, isIndex);
         Scalar potentialDiffNw = cellData.fluxData().upwindPotential(nPhaseIdx, isIndex);
diff --git a/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/3dvelocity.hh b/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/3dvelocity.hh
index 8aef1f4cb8e9b5f0a868ed53fcff349b6f268337..818f50b71bd314160b1ab206a1bfcb082ce4b9f2 100644
--- a/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/3dvelocity.hh
+++ b/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/3dvelocity.hh
@@ -29,6 +29,8 @@
 #include <dumux/porousmediumflow/sequential/cellcentered/mpfa/properties.hh>
 #include "3dtransmissibilitycalculator.hh"
 
+#include <dumux/common/deprecated.hh>
+
 namespace Dumux {
 
 /*!
@@ -63,7 +65,6 @@ template<class TypeTag> class FvMpfaL3dVelocity2p
     using Problem = GetPropType<TypeTag, Properties::Problem>;
 
     using SpatialParams = GetPropType<TypeTag, Properties::SpatialParams>;
-    using MaterialLaw = typename SpatialParams::MaterialLaw;
 
     using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
 
@@ -1972,19 +1973,20 @@ void FvMpfaL3dVelocity2p<TypeTag>::calculateBoundaryInteractionVolumeVelocity(In
 
                 }
 
-                Scalar pcBound = MaterialLaw::pc(
-                        problem_.spatialParams().materialLawParams(element), satWBound);
+                // old material law interface is deprecated: Replace this by
+                // const auto& fluidMatrixInteraction = spatialParams.fluidMatrixInteractionAtPos(element.geometry().center());
+                // after the release of 3.3, when the deprecated interface is no longer supported
+                const auto fluidMatrixInteraction = Deprecated::makePcKrSw(Scalar{}, problem_.spatialParams(), element);
+
+                Scalar pcBound = fluidMatrixInteraction.pc(satWBound);
 
                 Scalar gravityDiffBound = (problem_.bBoxMax() - globalPosFace) * gravity_
                         * (density_[nPhaseIdx] - density_[wPhaseIdx]);
 
                 pcBound += gravityDiffBound;
 
-                Dune::FieldVector<Scalar, numPhases> lambdaBound(
-                        MaterialLaw::krw(problem_.spatialParams().materialLawParams(element),
-                                satWBound));
-                lambdaBound[nPhaseIdx] = MaterialLaw::krn(
-                        problem_.spatialParams().materialLawParams(element), satWBound);
+                Dune::FieldVector<Scalar, numPhases> lambdaBound(fluidMatrixInteraction.krw(satWBound));
+                lambdaBound[nPhaseIdx] = fluidMatrixInteraction.krn(satWBound);
                 lambdaBound[wPhaseIdx] /= viscosity_[wPhaseIdx];
                 lambdaBound[nPhaseIdx] /= viscosity_[nPhaseIdx];
 
diff --git a/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/3dvelocityadaptive.hh b/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/3dvelocityadaptive.hh
index 9aa9447370a29fdfb15e209433c24f3d04073764..1afe847c7f7f7f275196c7253d25cb14402cdad3 100644
--- a/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/3dvelocityadaptive.hh
+++ b/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/3dvelocityadaptive.hh
@@ -65,7 +65,6 @@ template<class TypeTag> class FvMpfaL3dVelocity2pAdaptive: public FvMpfaL3dVeloc
     using Problem = GetPropType<TypeTag, Properties::Problem>;
 
     using SpatialParams = GetPropType<TypeTag, Properties::SpatialParams>;
-    using MaterialLaw = typename SpatialParams::MaterialLaw;
 
     using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
 
diff --git a/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/omethod/2dpressure.hh b/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/omethod/2dpressure.hh
index fe1b2e99703644bb228bdb06c0b2a259477c47f7..97f3be5fe571963888641670ad3c3898ea8337c0 100644
--- a/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/omethod/2dpressure.hh
+++ b/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/omethod/2dpressure.hh
@@ -28,6 +28,8 @@
 #include <dumux/porousmediumflow/2p/sequential/diffusion/properties.hh>
 #include <dumux/porousmediumflow/sequential/cellcentered/mpfa/properties.hh>
 
+#include <dumux/common/deprecated.hh>
+
 namespace Dumux {
 /*!
  * \ingroup SequentialTwoPModel
@@ -77,7 +79,6 @@ class FvMpfaO2dPressure2p: public FVPressure<TypeTag>
     using Problem = GetPropType<TypeTag, Properties::Problem>;
 
     using SpatialParams = GetPropType<TypeTag, Properties::SpatialParams>;
-    using MaterialLaw = typename SpatialParams::MaterialLaw;
 
     using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
 
@@ -1763,19 +1764,20 @@ void FvMpfaO2dPressure2p<TypeTag>::assemble()
 
                             }
 
-                            Scalar pcBound = MaterialLaw::pc(
-                                    problem_.spatialParams().materialLawParams(element), satWBound);
+                            // old material law interface is deprecated: Replace this by
+                            // const auto& fluidMatrixInteraction = spatialParams.fluidMatrixInteractionAtPos(element.geometry().center());
+                            // after the release of 3.3, when the deprecated interface is no longer supported
+                            const auto fluidMatrixInteraction = Deprecated::makePcKrSw(Scalar{}, problem_.spatialParams(), element);
+
+                            Scalar pcBound = fluidMatrixInteraction.pc(satWBound);
 
                             Scalar gravityDiffBound = (problem_.bBoxMax() - globalPosFace) * gravity_
                                     * (density_[nPhaseIdx] - density_[wPhaseIdx]);
 
                             pcBound += gravityDiffBound;
 
-                            Dune::FieldVector<Scalar, numPhases> lambdaBound(
-                                    MaterialLaw::krw(problem_.spatialParams().materialLawParams(element),
-                                            satWBound));
-                            lambdaBound[nPhaseIdx] = MaterialLaw::krn(
-                                    problem_.spatialParams().materialLawParams(element), satWBound);
+                            Dune::FieldVector<Scalar, numPhases> lambdaBound(fluidMatrixInteraction.krw(satWBound));
+                            lambdaBound[nPhaseIdx] = fluidMatrixInteraction.krn(satWBound);
                             lambdaBound[wPhaseIdx] /= viscosity_[wPhaseIdx];
                             lambdaBound[nPhaseIdx] /= viscosity_[nPhaseIdx];
 
@@ -1938,17 +1940,20 @@ void FvMpfaO2dPressure2p<TypeTag>::updateMaterialLaws()
 
         CellData& cellData = problem_.variables().cellData(eIdxGlobal);
 
-        Scalar satW = cellData.saturation(wPhaseIdx);
+        const Scalar satW = cellData.saturation(wPhaseIdx);
+
+        // old material law interface is deprecated: Replace this by
+        // const auto& fluidMatrixInteraction = spatialParams.fluidMatrixInteractionAtPos(element.geometry().center());
+        // after the release of 3.3, when the deprecated interface is no longer supported
+        const auto fluidMatrixInteraction = Deprecated::makePcKrSw(Scalar{}, problem_.spatialParams(), element);
 
-        Scalar pc = MaterialLaw::pc(problem_.spatialParams().materialLawParams(element), satW);
+        const Scalar pc = fluidMatrixInteraction.pc(satW);
 
         cellData.setCapillaryPressure(pc);
 
         // initialize mobilities
-        Scalar mobilityW = MaterialLaw::krw(problem_.spatialParams().materialLawParams(element), satW)
-                / viscosity_[wPhaseIdx];
-        Scalar mobilityNw = MaterialLaw::krn(problem_.spatialParams().materialLawParams(element), satW)
-                / viscosity_[nPhaseIdx];
+        const Scalar mobilityW = fluidMatrixInteraction.krw(satW) / viscosity_[wPhaseIdx];
+        const Scalar mobilityNw = fluidMatrixInteraction.krn(satW) / viscosity_[nPhaseIdx];
 
         // initialize mobilities
         cellData.setMobility(wPhaseIdx, mobilityW);
diff --git a/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/omethod/2dpressurevelocity.hh b/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/omethod/2dpressurevelocity.hh
index a38f5e39d1f612ae3e307dd04d45b0b4c1a10ffd..eb00f5a0777fd038b11ec6af18aa8747356476af 100644
--- a/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/omethod/2dpressurevelocity.hh
+++ b/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/omethod/2dpressurevelocity.hh
@@ -27,6 +27,8 @@
 #include "2dpressure.hh"
 #include "2dvelocity.hh"
 
+#include <dumux/common/deprecated.hh>
+
 namespace Dumux {
 /*!
  * \ingroup SequentialTwoPModel
@@ -60,7 +62,6 @@ template<class TypeTag> class FvMpfaO2dPressureVelocity2p: public FvMpfaO2dPress
     using PrimaryVariables = typename SolutionTypes::PrimaryVariables;
 
     using SpatialParams = GetPropType<TypeTag, Properties::SpatialParams>;
-    using MaterialLaw = typename SpatialParams::MaterialLaw;
 
     using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
     using FluidState = GetPropType<TypeTag, Properties::FluidState>;
@@ -454,8 +455,14 @@ void FvMpfaO2dPressureVelocity2p<TypeTag>::calculateVelocityOnBoundary(const Int
             satW = cellData.saturation(wPhaseIdx);
         }
 
-        Scalar pressBound = boundValues[pressureIdx];
-        Scalar pcBound = MaterialLaw::pc(problem_.spatialParams().materialLawParams(element), satW);
+        const Scalar pressBound = boundValues[pressureIdx];
+
+        // old material law interface is deprecated: Replace this by
+        // const auto& fluidMatrixInteraction = spatialParams.fluidMatrixInteractionAtPos(element.geometry().center());
+        // after the release of 3.3, when the deprecated interface is no longer supported
+        const auto fluidMatrixInteraction = Deprecated::makePcKrSw(Scalar{}, problem_.spatialParams(), element);
+
+        const Scalar pcBound = fluidMatrixInteraction.pc(satW);
 
         //determine phase pressures from primary pressure variable
         Scalar pressWBound = 0;
@@ -471,10 +478,8 @@ void FvMpfaO2dPressureVelocity2p<TypeTag>::calculateVelocityOnBoundary(const Int
             pressNwBound = pressBound;
         }
 
-        Scalar lambdaWBound = MaterialLaw::krw(problem_.spatialParams().materialLawParams(element), satW)
-                / viscosity_[wPhaseIdx];
-        Scalar lambdaNwBound = MaterialLaw::krn(problem_.spatialParams().materialLawParams(element), satW)
-                / viscosity_[nPhaseIdx];
+        const Scalar lambdaWBound = fluidMatrixInteraction.krw(satW) / viscosity_[wPhaseIdx];
+        const Scalar lambdaNwBound = fluidMatrixInteraction.krn(satW) / viscosity_[nPhaseIdx];
 
         Scalar potentialDiffW = cellData.fluxData().upwindPotential(wPhaseIdx, isIndex);
         Scalar potentialDiffNw = cellData.fluxData().upwindPotential(nPhaseIdx, isIndex);
diff --git a/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/omethod/2dvelocity.hh b/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/omethod/2dvelocity.hh
index b9ede8717aa6b397cd3b9d74621f2bb994007ee4..36a08b2bd7e74cb075c15fbfbdc338d1478a8a28 100644
--- a/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/omethod/2dvelocity.hh
+++ b/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/omethod/2dvelocity.hh
@@ -30,6 +30,8 @@
 #include <dumux/porousmediumflow/sequential/cellcentered/mpfa/properties.hh>
 #include <dumux/porousmediumflow/sequential/cellcentered/mpfa/ointeractionvolume.hh>
 
+#include <dumux/common/deprecated.hh>
+
 namespace Dumux {
 /*!
  * \ingroup SequentialTwoPModel
@@ -65,7 +67,6 @@ template<class TypeTag> class FvMpfaO2dVelocity2P
     using Problem = GetPropType<TypeTag, Properties::Problem>;
 
     using SpatialParams = GetPropType<TypeTag, Properties::SpatialParams>;
-    using MaterialLaw = typename SpatialParams::MaterialLaw;
 
     using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
 
@@ -619,21 +620,20 @@ void FvMpfaO2dVelocity2P<TypeTag>::calculateBoundaryInteractionVolumeVelocity(In
 
                     }
 
-                    Scalar pcBound = MaterialLaw::pc(
-                            problem_.spatialParams().materialLawParams(element), satWBound);
+                    // old material law interface is deprecated: Replace this by
+                    // const auto& fluidMatrixInteraction = spatialParams.fluidMatrixInteractionAtPos(element.geometry().center());
+                    // after the release of 3.3, when the deprecated interface is no longer supported
+                    const auto fluidMatrixInteraction = Deprecated::makePcKrSw(Scalar{}, problem_.spatialParams(), element);
+
+                    Scalar pcBound = fluidMatrixInteraction.pc(satWBound);
 
                     Scalar gravityDiffBound = (problem_.bBoxMax() - globalPosFace) * gravity_
                             * (density_[nPhaseIdx] - density_[wPhaseIdx]);
 
                     pcBound += gravityDiffBound;
 
-                    Dune::FieldVector < Scalar, numPhases
-                            > lambdaBound(
-                                    MaterialLaw::krw(
-                                            problem_.spatialParams().materialLawParams(element),
-                                            satWBound));
-                    lambdaBound[nPhaseIdx] = MaterialLaw::krn(
-                            problem_.spatialParams().materialLawParams(element), satWBound);
+                    Dune::FieldVector < Scalar, numPhases> lambdaBound(fluidMatrixInteraction.krw(satWBound));
+                    lambdaBound[nPhaseIdx] = fluidMatrixInteraction.krn(satWBound);
                     lambdaBound[wPhaseIdx] /= viscosity_[wPhaseIdx];
                     lambdaBound[nPhaseIdx] /= viscosity_[nPhaseIdx];
 
diff --git a/dumux/porousmediumflow/2p/sequential/transport/cellcentered/capillarydiffusion.hh b/dumux/porousmediumflow/2p/sequential/transport/cellcentered/capillarydiffusion.hh
index 4e50446756bcdb6b3aad9495fca7b06f0b62a08e..dfdc2b3ae04831af4b78662d2a6f4d99e936bb98 100644
--- a/dumux/porousmediumflow/2p/sequential/transport/cellcentered/capillarydiffusion.hh
+++ b/dumux/porousmediumflow/2p/sequential/transport/cellcentered/capillarydiffusion.hh
@@ -27,6 +27,8 @@
 #include <dumux/porousmediumflow/2p/sequential/transport/cellcentered/diffusivepart.hh>
 #include "properties.hh"
 
+#include <dumux/common/deprecated.hh>
+
 namespace Dumux {
 /*!
  * \ingroup SequentialTwoPModel
@@ -54,7 +56,6 @@ private:
       using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
 
       using SpatialParams = GetPropType<TypeTag, Properties::SpatialParams>;
-      using MaterialLaw = typename SpatialParams::MaterialLaw;
 
       using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
       using FluidState = GetPropType<TypeTag, Properties::FluidState>;
@@ -112,6 +113,11 @@ public:
         Scalar mobilityWI = 0;
         Scalar mobilityNwI = 0;
 
+        // old material law interface is deprecated: Replace this by
+        // const auto& fluidMatrixInteraction = spatialParams.fluidMatrixInteractionAtPos(element.geometry().center());
+        // after the release of 3.3, when the deprecated interface is no longer supported
+        const auto fluidMatrixInteraction = Deprecated::makePcKrSw(Scalar{}, problem_.spatialParams(), element);
+
         if (preComput_)
         {
             mobilityWI = CellDataI.mobility(wPhaseIdx);
@@ -123,9 +129,9 @@ public:
             fluidState.setPressure(wPhaseIdx, referencePressure);
             fluidState.setPressure(nPhaseIdx, referencePressure);
             fluidState.setTemperature(temperature);
-            mobilityWI = MaterialLaw::krw(problem_.spatialParams().materialLawParams(element), satI);
+            mobilityWI = fluidMatrixInteraction.krw(satI);
             mobilityWI /= FluidSystem::viscosity(fluidState, wPhaseIdx);
-            mobilityNwI = MaterialLaw::krn(problem_.spatialParams().materialLawParams(element), satI);
+            mobilityNwI = fluidMatrixInteraction.krn(satI);
             mobilityNwI /= FluidSystem::viscosity(fluidState, nPhaseIdx);
         }
 
@@ -172,9 +178,14 @@ public:
                 fluidState.setPressure(nPhaseIdx, referencePressure);
                 fluidState.setTemperature(temperature);
 
-                mobilityWJ = MaterialLaw::krw(problem_.spatialParams().materialLawParams(neighbor), satJ);
+                // old material law interface is deprecated: Replace this by
+                // const auto& fluidMatrixInteractionNeighbor = spatialParams.fluidMatrixInteractionAtPos(neighbor.geometry().center());
+                // after the release of 3.3, when the deprecated interface is no longer supported
+                const auto fluidMatrixInteractionNeighbor = Deprecated::makePcKrSw(Scalar{}, problem_.spatialParams(), neighbor);
+
+                mobilityWJ = fluidMatrixInteractionNeighbor.krw(satJ);
                 mobilityWJ /= FluidSystem::viscosity(fluidState, wPhaseIdx);
-                mobilityNwJ = MaterialLaw::krn(problem_.spatialParams().materialLawParams(neighbor), satJ);
+                mobilityNwJ = fluidMatrixInteractionNeighbor.krn(satJ);
                 mobilityNwJ /= FluidSystem::viscosity(fluidState, nPhaseIdx);
             }
             Scalar mobilityWMean = 0.5*(mobilityWI + mobilityWJ);
@@ -207,9 +218,9 @@ public:
             fluidState.setPressure(wPhaseIdx, referencePressure);
             fluidState.setPressure(nPhaseIdx, referencePressure);
             fluidState.setTemperature(temperature);
-            mobilityWJ = MaterialLaw::krw(problem_.spatialParams().materialLawParams(element), satJ);
+            mobilityWJ = fluidMatrixInteraction.krw(satJ);
             mobilityWJ /= FluidSystem::viscosity(fluidState, wPhaseIdx);
-            mobilityNwJ = MaterialLaw::krn(problem_.spatialParams().materialLawParams(element), satJ);
+            mobilityNwJ = fluidMatrixInteraction.krn(satJ);
             mobilityNwJ /= FluidSystem::viscosity(fluidState, nPhaseIdx);
 
             Scalar mobWMean = 0.5 * (mobilityWI + mobilityWJ);
diff --git a/dumux/porousmediumflow/2p/sequential/transport/cellcentered/evalcflfluxcoats.hh b/dumux/porousmediumflow/2p/sequential/transport/cellcentered/evalcflfluxcoats.hh
index 0f43ac4153a351d4c83048a9cc62f62a903cefec..4750fe67abebe790db361ed4a4bc5cb83ff5e9b4 100644
--- a/dumux/porousmediumflow/2p/sequential/transport/cellcentered/evalcflfluxcoats.hh
+++ b/dumux/porousmediumflow/2p/sequential/transport/cellcentered/evalcflfluxcoats.hh
@@ -28,6 +28,8 @@
 #include <dumux/porousmediumflow/sequential/impetproperties.hh>
 #include "evalcflflux.hh"
 
+#include <dumux/common/deprecated.hh>
+
 namespace Dumux {
 /*!
  * \ingroup SequentialTwoPModel
@@ -44,7 +46,6 @@ private:
     using Problem = GetPropType<TypeTag, Properties::Problem>;
 
     using SpatialParams = GetPropType<TypeTag, Properties::SpatialParams>;
-    using MaterialLaw = typename SpatialParams::MaterialLaw;
     using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
     using FluidState = GetPropType<TypeTag, Properties::FluidState>;
     using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
@@ -357,7 +358,12 @@ void EvalCflFluxCoats<TypeTag>::addCoatsFlux(Scalar& lambdaW, Scalar& lambdaNw,
     Scalar lambdaWI = cellDataI.mobility(wPhaseIdx);
     Scalar lambdaNwI = cellDataI.mobility(nPhaseIdx);
 
-    Scalar dpc_dsI = MaterialLaw::dpc_dsw(problem_.spatialParams().materialLawParams(element), satI);
+    // old material law interface is deprecated: Replace this by
+    // const auto& fluidMatrixInteraction = spatialParams.fluidMatrixInteractionAtPos(element.geometry().center());
+    // after the release of 3.3, when the deprecated interface is no longer supported
+    const auto fluidMatrixInteraction = Deprecated::makePcKrSw(Scalar{}, problem_.spatialParams(), element);
+
+    const Scalar dpc_dsI = fluidMatrixInteraction.dpc_dsw(satI);
 
     const GlobalPosition& unitOuterNormal = intersection.centerUnitOuterNormal();
 
@@ -407,7 +413,12 @@ void EvalCflFluxCoats<TypeTag>::addCoatsFlux(Scalar& lambdaW, Scalar& lambdaNw,
             Scalar lambdaWJ = cellDataI.mobility(wPhaseIdx);
             Scalar lambdaNwJ = cellDataI.mobility(nPhaseIdx);
 
-            Scalar dpc_dsJ = MaterialLaw::dpc_dsw(problem_.spatialParams().materialLawParams(neighbor), satJ);
+            // old material law interface is deprecated: Replace this by
+            // const auto& fluidMatrixInteraction = spatialParams.fluidMatrixInteractionAtPos(neighbor.geometry().center());
+            // after the release of 3.3, when the deprecated interface is no longer supported
+            const auto fluidMatrixInteractionNeighbor = Deprecated::makePcKrSw(Scalar{}, problem_.spatialParams(), neighbor);
+
+            const Scalar dpc_dsJ = fluidMatrixInteractionNeighbor.dpc_dsw(satJ);
 
             // compute vectorized permeabilities
             DimVector permeability(0);
@@ -448,8 +459,8 @@ void EvalCflFluxCoats<TypeTag>::addCoatsFlux(Scalar& lambdaW, Scalar& lambdaNw,
                 ds += epsDerivative_;
             }
 
-            Scalar dLambdaWDs = MaterialLaw::krw(problem_.spatialParams().materialLawParams(neighbor), abs(satPlus)) / viscosityW;
-            dLambdaWDs -= MaterialLaw::krw(problem_.spatialParams().materialLawParams(neighbor), abs(satMinus)) / viscosityW;
+            Scalar dLambdaWDs = fluidMatrixInteractionNeighbor.krw(abs(satPlus)) / viscosityW;
+            dLambdaWDs -= fluidMatrixInteractionNeighbor.krw(abs(satMinus)) / viscosityW;
             dLambdaWDs /= (ds);
 
             if (upwindNwI)
@@ -471,8 +482,8 @@ void EvalCflFluxCoats<TypeTag>::addCoatsFlux(Scalar& lambdaW, Scalar& lambdaNw,
                 ds += epsDerivative_;
             }
 
-            Scalar dLambdaNwDs = MaterialLaw::krn(problem_.spatialParams().materialLawParams(neighbor), satPlus) / viscosityNw;
-            dLambdaNwDs -= MaterialLaw::krn(problem_.spatialParams().materialLawParams(neighbor), satMinus) / viscosityNw;
+            Scalar dLambdaNwDs = fluidMatrixInteractionNeighbor.krn(satPlus) / viscosityNw;
+            dLambdaNwDs -= fluidMatrixInteractionNeighbor.krn(satMinus) / viscosityNw;
             dLambdaNwDs /= (ds);
 
             Scalar lambdaWCap = 0.5 * (lambdaWI + lambdaWJ);
@@ -559,13 +570,13 @@ void EvalCflFluxCoats<TypeTag>::addCoatsFlux(Scalar& lambdaW, Scalar& lambdaNw,
                     case pw:
                     {
                         potWBound = bcValues[eqIdxPress] + density_[wPhaseIdx] * gdeltaZ;
-                        potNwBound = bcValues[eqIdxPress] + MaterialLaw::pc(problem_.spatialParams().materialLawParams(element), satWBound)
+                        potNwBound = bcValues[eqIdxPress] + fluidMatrixInteraction.pc(satWBound)
                                                           + density_[nPhaseIdx] * gdeltaZ;
                         break;
                     }
                     case pn:
                     {
-                        potWBound = bcValues[eqIdxPress] - MaterialLaw::pc(problem_.spatialParams().materialLawParams(element),satWBound)
+                        potWBound = bcValues[eqIdxPress] - fluidMatrixInteraction.pc(satWBound)
                                                          + density_[wPhaseIdx] * gdeltaZ;
                         potNwBound = bcValues[eqIdxPress] + density_[nPhaseIdx] * gdeltaZ;
                         break;
@@ -602,12 +613,12 @@ void EvalCflFluxCoats<TypeTag>::addCoatsFlux(Scalar& lambdaW, Scalar& lambdaNw,
 
                 if (hasPotWBound && !hasPotNwBound)
                 {
-                    potNwBound = potWBound + MaterialLaw::pc(problem_.spatialParams().materialLawParams(element),satWBound)
+                    potNwBound = potWBound + fluidMatrixInteraction.pc(satWBound)
                                            + (density_[nPhaseIdx] - density_[wPhaseIdx]) * gdeltaZ;
                 }
                 else if (!hasPotWBound && hasPotNwBound)
                 {
-                   potWBound = potNwBound - MaterialLaw::pc(problem_.spatialParams().materialLawParams(element),satWBound)
+                   potWBound = potNwBound - fluidMatrixInteraction.pc(satWBound)
                                           + (density_[nPhaseIdx] - density_[wPhaseIdx]) * gdeltaZ;
                 }
             }
@@ -651,7 +662,7 @@ void EvalCflFluxCoats<TypeTag>::addCoatsFlux(Scalar& lambdaW, Scalar& lambdaNw,
                 return;
             }
 
-            Scalar dpc_dsBound = MaterialLaw::dpc_dsw(problem_.spatialParams().materialLawParams(element), satWBound);
+            const Scalar dpc_dsBound = fluidMatrixInteraction.dpc_dsw(satWBound);
 
             Scalar lambdaWBound = 0;
             Scalar lambdaNwBound = 0;
@@ -666,8 +677,8 @@ void EvalCflFluxCoats<TypeTag>::addCoatsFlux(Scalar& lambdaW, Scalar& lambdaNw,
             Scalar viscosityWBound = FluidSystem::viscosity(fluidState, wPhaseIdx);
             Scalar viscosityNwBound =
                 FluidSystem::viscosity(fluidState, nPhaseIdx);
-            lambdaWBound = MaterialLaw::krw(problem_.spatialParams().materialLawParams(element), satWBound) / viscosityWBound;
-            lambdaNwBound = MaterialLaw::krn(problem_.spatialParams().materialLawParams(element), satWBound) / viscosityNwBound;
+            lambdaWBound = fluidMatrixInteraction.krw(satWBound) / viscosityWBound;
+            lambdaNwBound = fluidMatrixInteraction.krn(satWBound) / viscosityNwBound;
 
             Scalar satUpw = 0;
             using std::max;
@@ -690,8 +701,8 @@ void EvalCflFluxCoats<TypeTag>::addCoatsFlux(Scalar& lambdaW, Scalar& lambdaNw,
                 ds += epsDerivative_;
             }
 
-            Scalar dLambdaWDs = MaterialLaw::krw(problem_.spatialParams().materialLawParams(element), satPlus) / viscosityW;
-            dLambdaWDs -= MaterialLaw::krw(problem_.spatialParams().materialLawParams(element), satMinus) / viscosityW;
+            Scalar dLambdaWDs = fluidMatrixInteraction.krw(satPlus) / viscosityW;
+            dLambdaWDs -= fluidMatrixInteraction.krw(satMinus) / viscosityW;
             dLambdaWDs /= (ds);
 
             if (cellDataI.fluxData().isUpwindCell(nPhaseIdx, indexInInside))
@@ -713,8 +724,8 @@ void EvalCflFluxCoats<TypeTag>::addCoatsFlux(Scalar& lambdaW, Scalar& lambdaNw,
                 ds += epsDerivative_;
             }
 
-            Scalar dLambdaNwDs = MaterialLaw::krn(problem_.spatialParams().materialLawParams(element), satPlus) / viscosityNw;
-            dLambdaNwDs -= MaterialLaw::krn(problem_.spatialParams().materialLawParams(element), satMinus) / viscosityNw;
+            Scalar dLambdaNwDs = fluidMatrixInteraction.krn(satPlus) / viscosityNw;
+            dLambdaNwDs -= fluidMatrixInteraction.krn(satMinus) / viscosityNw;
             dLambdaNwDs /= (ds);
 
             Scalar lambdaWCap = 0.5 * (lambdaWI + lambdaWBound);
diff --git a/dumux/porousmediumflow/2p/sequential/transport/cellcentered/evalcflfluxdefault.hh b/dumux/porousmediumflow/2p/sequential/transport/cellcentered/evalcflfluxdefault.hh
index ae2e2ac4ec4378b80aa7d3f1ecf13de490ef2247..d249715a254cfda8d4c7699b151444df743afb3a 100644
--- a/dumux/porousmediumflow/2p/sequential/transport/cellcentered/evalcflfluxdefault.hh
+++ b/dumux/porousmediumflow/2p/sequential/transport/cellcentered/evalcflfluxdefault.hh
@@ -27,6 +27,8 @@
 #include <dumux/porousmediumflow/sequential/impetproperties.hh>
 #include "evalcflflux.hh"
 
+#include <dumux/common/deprecated.hh>
+
 namespace Dumux {
 /*!
  * \ingroup SequentialTwoPModel
@@ -234,11 +236,16 @@ private:
 template<class TypeTag>
 typename EvalCflFluxDefault<TypeTag>::Scalar EvalCflFluxDefault<TypeTag>::getCflFluxFunction(const Element& element)
 {
-    Scalar residualSatW = problem_.spatialParams().materialLawParams(element).swr();
-    Scalar residualSatNw = problem_.spatialParams().materialLawParams(element).snr();
+    // old material law interface is deprecated: Replace this by
+    // const auto& fluidMatrixInteraction = spatialParams.fluidMatrixInteractionAtPos(element.geometry().center());
+    // after the release of 3.3, when the deprecated interface is no longer supported
+    const auto fluidMatrixInteraction = Deprecated::makePcKrSw(Scalar{}, problem_.spatialParams(), element);
+
+    const Scalar residualSatW = fluidMatrixInteraction.pcSwCurve().effToAbsParams().swr();
+    const Scalar residualSatNw = fluidMatrixInteraction.pcSwCurve().effToAbsParams().snr();
 
     // compute dt restriction
-    Scalar volumeCorrectionFactor = 1 - residualSatW - residualSatNw;
+    const Scalar volumeCorrectionFactor = 1 - residualSatW - residualSatNw;
     Scalar volumeCorrectionFactorOutW = 0;
     Scalar volumeCorrectionFactorOutNw = 0;
 
@@ -273,7 +280,7 @@ typename EvalCflFluxDefault<TypeTag>::Scalar EvalCflFluxDefault<TypeTag>::getCfl
     }
 
     //determine timestep
-    Scalar cFLFluxFunction = min(cFLFluxIn, cFLFluxOut);
+    const Scalar cFLFluxFunction = min(cFLFluxIn, cFLFluxOut);
 
     return cFLFluxFunction;
 }
diff --git a/dumux/porousmediumflow/2p/sequential/transport/cellcentered/gravitypart.hh b/dumux/porousmediumflow/2p/sequential/transport/cellcentered/gravitypart.hh
index 59cb77ff6a6da2c8af82a0e1b8351df8f849b061..2e2eab73755495943755449ec9b73319abdd6bf7 100644
--- a/dumux/porousmediumflow/2p/sequential/transport/cellcentered/gravitypart.hh
+++ b/dumux/porousmediumflow/2p/sequential/transport/cellcentered/gravitypart.hh
@@ -27,6 +27,8 @@
 #include <dumux/porousmediumflow/2p/sequential/transport/cellcentered/convectivepart.hh>
 #include "properties.hh"
 
+#include <dumux/common/deprecated.hh>
+
 namespace Dumux {
 /*!
  * \ingroup SequentialTwoPModel
@@ -54,7 +56,6 @@ private:
     using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
 
     using SpatialParams = GetPropType<TypeTag, Properties::SpatialParams>;
-    using MaterialLaw = typename SpatialParams::MaterialLaw;
 
     using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
     using FluidState = GetPropType<TypeTag, Properties::FluidState>;
@@ -100,6 +101,11 @@ public:
         Scalar lambdaWJ = 0;
         Scalar lambdaNwJ = 0;
 
+        // old material law interface is deprecated: Replace this by
+        // const auto& fluidMatrixInteraction = spatialParams.fluidMatrixInteractionAtPos(element.geometry().center());
+        // after the release of 3.3, when the deprecated interface is no longer supported
+        const auto fluidMatrixInteraction = Deprecated::makePcKrSw(Scalar{}, problem_.spatialParams(), element);
+
         if (preComput_)
         {
             lambdaWI=cellDataI.mobility(wPhaseIdx);
@@ -107,9 +113,9 @@ public:
         }
         else
         {
-            lambdaWI = MaterialLaw::krw(problem_.spatialParams().materialLawParams(element), satI);
+            lambdaWI = fluidMatrixInteraction.krw(satI);
             lambdaWI /= viscosity_[wPhaseIdx];
-            lambdaNwI = MaterialLaw::krn(problem_.spatialParams().materialLawParams(element), satI);
+            lambdaNwI = fluidMatrixInteraction.krn(satI);
             lambdaNwI /= viscosity_[nPhaseIdx];
         }
 
@@ -144,9 +150,14 @@ public:
             }
             else
             {
-                lambdaWJ = MaterialLaw::krw(problem_.spatialParams().materialLawParams(neighbor), satJ);
+                // old material law interface is deprecated: Replace this by
+                // const auto& fluidMatrixInteraction = spatialParams.fluidMatrixInteractionAtPos(neighbor.geometry().center());
+                // after the release of 3.3, when the deprecated interface is no longer supported
+                const auto fluidMatrixInteractionNeighbor = Deprecated::makePcKrSw(Scalar{}, problem_.spatialParams(), neighbor);
+
+                lambdaWJ = fluidMatrixInteractionNeighbor.krw(satJ);
                 lambdaWJ /= viscosity_[wPhaseIdx];
-                lambdaNwJ = MaterialLaw::krn(problem_.spatialParams().materialLawParams(neighbor), satJ);
+                lambdaNwJ = fluidMatrixInteractionNeighbor.krn(satJ);
                 lambdaNwJ /= viscosity_[nPhaseIdx];
             }
 
@@ -164,9 +175,9 @@ public:
             distVec = intersection.geometry().center() - element.geometry().center();
 
             //calculate lambda_n*f_w at the boundary
-            lambdaWJ = MaterialLaw::krw(problem_.spatialParams().materialLawParams(element), satJ);
+            lambdaWJ = fluidMatrixInteraction.krw(satJ);
             lambdaWJ /= viscosity_[wPhaseIdx];
-            lambdaNwJ = MaterialLaw::krn(problem_.spatialParams().materialLawParams(element), satJ);
+            lambdaNwJ = fluidMatrixInteraction.krn(satJ);
             lambdaNwJ /= viscosity_[nPhaseIdx];
 
             //If potential is zero always take value from the boundary!
diff --git a/dumux/porousmediumflow/2p/sequential/transport/cellcentered/saturation.hh b/dumux/porousmediumflow/2p/sequential/transport/cellcentered/saturation.hh
index d8331db50676270f8242d0e37dfeb673e0575f52..da29f8af87cde5f5c87d0698867619ea9b6a47f5 100644
--- a/dumux/porousmediumflow/2p/sequential/transport/cellcentered/saturation.hh
+++ b/dumux/porousmediumflow/2p/sequential/transport/cellcentered/saturation.hh
@@ -27,6 +27,8 @@
 #include <dumux/porousmediumflow/2p/sequential/transport/properties.hh>
 #include <dumux/porousmediumflow/sequential/cellcentered/transport.hh>
 
+#include <dumux/common/deprecated.hh>
+
 namespace Dumux {
 
 /*!
@@ -92,7 +94,6 @@ class FVSaturation2P: public FVTransport<TypeTag>
     using GravityFlux = GetPropType<TypeTag, Properties::GravityFlux>;
 
     using SpatialParams = GetPropType<TypeTag, Properties::SpatialParams>;
-    using MaterialLaw = typename SpatialParams::MaterialLaw;
 
     using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
 
@@ -803,7 +804,12 @@ void FVSaturation2P<TypeTag>::getFluxOnBoundary(Scalar& update, const Intersecti
         }
         }
 
-        Scalar pcBound = MaterialLaw::pc(problem_.spatialParams().materialLawParams(elementI), satWBound);
+        // old material law interface is deprecated: Replace this by
+        // const auto& fluidMatrixInteraction = spatialParams.fluidMatrixInteractionAtPos(elementI.geometry().center());
+        // after the release of 3.3, when the deprecated interface is no longer supported
+        const auto fluidMatrixInteraction = Deprecated::makePcKrSw(Scalar{}, problem_.spatialParams(), elementI);
+
+        const Scalar pcBound = fluidMatrixInteraction.pc(satWBound);
 
         Scalar lambdaW = 0;
         Scalar lambdaNw = 0;
@@ -821,12 +827,12 @@ void FVSaturation2P<TypeTag>::getFluxOnBoundary(Scalar& update, const Intersecti
         {
             if (compressibility_)
             {
-                lambdaW = MaterialLaw::krw(problem_.spatialParams().materialLawParams(elementI), satWBound)
+                lambdaW = fluidMatrixInteraction.krw(satWBound)
                         / FluidSystem::viscosity(cellDataI.fluidState(), wPhaseIdx);
             }
             else
             {
-                lambdaW = MaterialLaw::krw(problem_.spatialParams().materialLawParams(elementI), satWBound)
+                lambdaW = fluidMatrixInteraction.krw(satWBound)
                         / viscosity_[wPhaseIdx];
             }
         }
@@ -843,12 +849,12 @@ void FVSaturation2P<TypeTag>::getFluxOnBoundary(Scalar& update, const Intersecti
         {
             if (compressibility_)
             {
-                lambdaNw = MaterialLaw::krn(problem_.spatialParams().materialLawParams(elementI), satWBound)
+                lambdaNw = fluidMatrixInteraction.krn(satWBound)
                         / FluidSystem::viscosity(cellDataI.fluidState(), nPhaseIdx);
             }
             else
             {
-                lambdaNw = MaterialLaw::krn(problem_.spatialParams().materialLawParams(elementI), satWBound)
+                lambdaNw = fluidMatrixInteraction.krn(satWBound)
                         / viscosity_[nPhaseIdx];
             }
         }
@@ -1196,13 +1202,18 @@ void FVSaturation2P<TypeTag>::updateMaterialLaws()
         //determine phase saturations from primary saturation variable
         Scalar satW = cellData.saturation(wPhaseIdx);
 
-        Scalar pc = MaterialLaw::pc(problem_.spatialParams().materialLawParams(element), satW);
+        // old material law interface is deprecated: Replace this by
+        // const auto& fluidMatrixInteraction = spatialParams.fluidMatrixInteractionAtPos(elementI.geometry().center());
+        // after the release of 3.3, when the deprecated interface is no longer supported
+        const auto fluidMatrixInteraction = Deprecated::makePcKrSw(Scalar{}, problem_.spatialParams(), element);
+
+        const Scalar pc = fluidMatrixInteraction.pc(satW);
 
         cellData.setCapillaryPressure(pc);
 
         // initialize mobilities
-        Scalar mobilityW = MaterialLaw::krw(problem_.spatialParams().materialLawParams(element), satW) / viscosity_[wPhaseIdx];
-        Scalar mobilityNw = MaterialLaw::krn(problem_.spatialParams().materialLawParams(element), satW) / viscosity_[nPhaseIdx];
+        const Scalar mobilityW = fluidMatrixInteraction.krw(satW) / viscosity_[wPhaseIdx];
+        const Scalar mobilityNw = fluidMatrixInteraction.krn(satW) / viscosity_[nPhaseIdx];
 
         // initialize mobilities
         cellData.setMobility(wPhaseIdx, mobilityW);