diff --git a/dumux/porousmediumflow/2p/incompressiblelocalresidual.hh b/dumux/porousmediumflow/2p/incompressiblelocalresidual.hh
index 5cd473eaf8bd5c6a6b9d924208d00932be7e3191..06779218e9ddad106957529fa7c294d7d32a55d6 100644
--- a/dumux/porousmediumflow/2p/incompressiblelocalresidual.hh
+++ b/dumux/porousmediumflow/2p/incompressiblelocalresidual.hh
@@ -38,6 +38,8 @@
 #include <dumux/porousmediumflow/immiscible/localresidual.hh>
 #include <dumux/porousmediumflow/2p/formulation.hh>
 
+#include <dumux/common/deprecated.hh>
+
 namespace Dumux {
 
 /*!
@@ -169,7 +171,6 @@ public:
         static_assert(ModelTraits::priVarFormulation() == TwoPFormulation::p0s1,
                       "2p/incompressiblelocalresidual.hh: Analytic differentiation has to be checked for p1-s0 formulation!");
 
-        using MaterialLaw = typename Problem::SpatialParams::MaterialLaw;
         using AdvectionType = GetPropType<TypeTag, Properties::AdvectionType>;
 
         // evaluate the current wetting phase Darcy flux and resulting upwind weights
@@ -189,12 +190,23 @@ public:
         const auto& outsideScv = fvGeometry.scv(outsideScvIdx);
         const auto& insideVolVars = curElemVolVars[insideScvIdx];
         const auto& outsideVolVars = curElemVolVars[outsideScvIdx];
-        const auto& insideMaterialParams = problem.spatialParams().materialLawParams(element,
-                                                                                     insideScv,
-                                                                                     elementSolution<FVElementGeometry>(insideVolVars.priVars()));
-        const auto& outsideMaterialParams = problem.spatialParams().materialLawParams(outsideElement,
-                                                                                      outsideScv,
-                                                                                      elementSolution<FVElementGeometry>(outsideVolVars.priVars()));
+
+        // old material law interface is deprecated: Replace this by
+        // const auto& insidefluidMatrixInteraction = problem.spatialParams().fluidMatrixInteraction(element,
+        //                                                                                           insideScv,
+        //                                                                                           elementSolution<FVElementGeometry>(insideVolVars.priVars()));
+        // const auto& outsidefluidMatrixInteraction = problem.spatialParams().fluidMatrixInteraction(outsideElement,
+        //                                                                                            outsideScv,
+        //                                                                                            elementSolution<FVElementGeometry>(outsideVolVars.priVars()));
+        // after the release of 3.3, when the deprecated interface is no longer supported
+        const auto& insidefluidMatrixInteraction = Deprecated::makePcKrSw(Scalar{}, problem.spatialParams(),
+                                                                          element,
+                                                                          insideScv,
+                                                                          elementSolution<FVElementGeometry>(insideVolVars.priVars()));
+        const auto& outsidefluidMatrixInteraction = Deprecated::makePcKrSw(Scalar{}, problem.spatialParams(),
+                                                                           outsideElement,
+                                                                           outsideScv,
+                                                                           elementSolution<FVElementGeometry>(outsideVolVars.priVars()));
 
         // get references to the two participating derivative matrices
         auto& dI_dI = derivativeMatrices[insideScvIdx];
@@ -213,12 +225,12 @@ public:
         // derivative w.r.t. to Sn is the negative of the one w.r.t. Sw
         const auto insideSw = insideVolVars.saturation(0);
         const auto outsideSw = outsideVolVars.saturation(0);
-        const auto dKrw_dSn_inside = -1.0*MaterialLaw::dkrw_dsw(insideMaterialParams, insideSw);
-        const auto dKrw_dSn_outside = -1.0*MaterialLaw::dkrw_dsw(outsideMaterialParams, outsideSw);
-        const auto dKrn_dSn_inside = -1.0*MaterialLaw::dkrn_dsw(insideMaterialParams, insideSw);
-        const auto dKrn_dSn_outside = -1.0*MaterialLaw::dkrn_dsw(outsideMaterialParams, outsideSw);
-        const auto dpc_dSn_inside = -1.0*MaterialLaw::dpc_dsw(insideMaterialParams, insideSw);
-        const auto dpc_dSn_outside = -1.0*MaterialLaw::dpc_dsw(outsideMaterialParams, outsideSw);
+        const auto dKrw_dSn_inside = -1.0*insidefluidMatrixInteraction.dkrw_dsw(insideSw);
+        const auto dKrw_dSn_outside = -1.0*outsidefluidMatrixInteraction.dkrw_dsw(outsideSw);
+        const auto dKrn_dSn_inside = -1.0*insidefluidMatrixInteraction.dkrn_dsw(insideSw);
+        const auto dKrn_dSn_outside = -1.0*outsidefluidMatrixInteraction.dkrn_dsw(outsideSw);
+        const auto dpc_dSn_inside = -1.0*insidefluidMatrixInteraction.dpc_dsw(insideSw);
+        const auto dpc_dSn_outside = -1.0*outsidefluidMatrixInteraction.dpc_dsw(outsideSw);
 
         const auto tij = elemFluxVarsCache[scvf].advectionTij();
 
@@ -288,7 +300,6 @@ public:
         static_assert(ModelTraits::priVarFormulation() == TwoPFormulation::p0s1,
                       "2p/incompressiblelocalresidual.hh: Analytic differentiation has to be checked for p0-s1 formulation!");
 
-        using MaterialLaw = typename Problem::SpatialParams::MaterialLaw;
         using AdvectionType = GetPropType<TypeTag, Properties::AdvectionType>;
 
         // evaluate the current wetting phase Darcy flux and resulting upwind weights
@@ -310,8 +321,22 @@ public:
 
         const auto elemSol = elementSolution(element, curElemVolVars, fvGeometry);
 
-        const auto& insideMaterialParams = problem.spatialParams().materialLawParams(element, insideScv, elemSol);
-        const auto& outsideMaterialParams = problem.spatialParams().materialLawParams(element, outsideScv, elemSol);
+        // old material law interface is deprecated: Replace this by
+        // const auto& insidefluidMatrixInteraction = problem.spatialParams().fluidMatrixInteraction(element,
+        //                                                                                           insideScv,
+        //                                                                                           elemSol);
+        // const auto& outsidefluidMatrixInteraction = problem.spatialParams().fluidMatrixInteraction(element,
+        //                                                                                            outsideScv,
+        //                                                                                            elemSol);
+        // after the release of 3.3, when the deprecated interface is no longer supported
+        const auto& insidefluidMatrixInteraction = Deprecated::makePcKrSw(Scalar{}, problem.spatialParams(),
+                                                                          element,
+                                                                          insideScv,
+                                                                          elemSol);
+        const auto& outsidefluidMatrixInteraction = Deprecated::makePcKrSw(Scalar{}, problem.spatialParams(),
+                                                                           element,
+                                                                           outsideScv,
+                                                                           elemSol);
 
         // some quantities to be reused (rho & mu are constant and thus equal for all cells)
         static const auto rho_w = insideVolVars.density(0);
@@ -366,19 +391,19 @@ public:
             {
                 // partial derivative of the wetting phase flux w.r.t. S_n
                 const auto insideSw = insideVolVars.saturation(0);
-                const auto dKrw_dSn_inside = -1.0*MaterialLaw::dkrw_dsw(insideMaterialParams, insideSw);
+                const auto dKrw_dSn_inside = -1.0*insidefluidMatrixInteraction.dkrw_dsw(insideSw);
                 const auto dFluxW_dSnJ = rho_mu_flux_w*dKrw_dSn_inside*insideWeight_w;
                 dI_dJ_inside[globalJ][conti0EqIdx+0][saturationIdx] += dFluxW_dSnJ;
                 dI_dJ_outside[globalJ][conti0EqIdx+0][saturationIdx] -= dFluxW_dSnJ;
 
                 // partial derivative of the nonwetting phase flux w.r.t. S_n (k_rn contribution)
-                const auto dKrn_dSn_inside = -1.0*MaterialLaw::dkrn_dsw(insideMaterialParams, insideSw);
+                const auto dKrn_dSn_inside = -1.0*insidefluidMatrixInteraction.dkrn_dsw(insideSw);
                 const auto dFluxN_dSnJ_krn = rho_mu_flux_n*dKrn_dSn_inside*insideWeight_n;
                 dI_dJ_inside[globalJ][conti0EqIdx+1][saturationIdx] += dFluxN_dSnJ_krn;
                 dI_dJ_outside[globalJ][conti0EqIdx+1][saturationIdx] -= dFluxN_dSnJ_krn;
 
                 // partial derivative of the nonwetting phase flux w.r.t. S_n (p_c contribution)
-                const auto dFluxN_dSnJ_pc = -1.0*tj_up_n*MaterialLaw::dpc_dsw(insideMaterialParams, insideSw);
+                const auto dFluxN_dSnJ_pc = -1.0*tj_up_n*insidefluidMatrixInteraction.dpc_dsw(insideSw);
                 dI_dJ_inside[globalJ][conti0EqIdx+1][saturationIdx] += dFluxN_dSnJ_pc;
                 dI_dJ_outside[globalJ][conti0EqIdx+1][saturationIdx] -= dFluxN_dSnJ_pc;
             }
@@ -386,25 +411,34 @@ public:
             {
                 // see comments for (globalJ == insideScvIdx)
                 const auto outsideSw = outsideVolVars.saturation(0);
-                const auto dKrw_dSn_outside = -1.0*MaterialLaw::dkrw_dsw(outsideMaterialParams, outsideSw);
+                const auto dKrw_dSn_outside = -1.0*outsidefluidMatrixInteraction.dkrw_dsw(outsideSw);
                 const auto dFluxW_dSnJ = rho_mu_flux_w*dKrw_dSn_outside*outsideWeight_w;
                 dI_dJ_inside[globalJ][conti0EqIdx+0][saturationIdx] += dFluxW_dSnJ;
                 dI_dJ_outside[globalJ][conti0EqIdx+0][saturationIdx] -= dFluxW_dSnJ;
 
-                const auto dKrn_dSn_outside = -1.0*MaterialLaw::dkrn_dsw(outsideMaterialParams, outsideSw);
+                const auto dKrn_dSn_outside = -1.0*outsidefluidMatrixInteraction.dkrn_dsw(outsideSw);
                 const auto dFluxN_dSnJ_krn = rho_mu_flux_n*dKrn_dSn_outside*outsideWeight_n;
                 dI_dJ_inside[globalJ][conti0EqIdx+1][saturationIdx] += dFluxN_dSnJ_krn;
                 dI_dJ_outside[globalJ][conti0EqIdx+1][saturationIdx] -= dFluxN_dSnJ_krn;
 
-                const auto dFluxN_dSnJ_pc = -1.0*tj_up_n*MaterialLaw::dpc_dsw(outsideMaterialParams, outsideSw);
+                const auto dFluxN_dSnJ_pc = -1.0*tj_up_n*outsidefluidMatrixInteraction.dpc_dsw(outsideSw);
                 dI_dJ_inside[globalJ][conti0EqIdx+1][saturationIdx] += dFluxN_dSnJ_pc;
                 dI_dJ_outside[globalJ][conti0EqIdx+1][saturationIdx] -= dFluxN_dSnJ_pc;
             }
             else
             {
-                const auto& paramsJ = problem.spatialParams().materialLawParams(element, scvJ, elemSol);
+                // old material law interface is deprecated: Replace this by
+                // const auto& fluidMatrixInteractionJ = problem.spatialParams().fluidMatrixInteraction(element,
+                //                                                                                      scvJ,
+                //                                                                                      elemSol);
+
+                // after the release of 3.3, when the deprecated interface is no longer supported
+                const auto& fluidMatrixInteractionJ = Deprecated::makePcKrSw(Scalar{}, problem.spatialParams(),
+                                                                             element,
+                                                                             scvJ,
+                                                                             elemSol);
                 const auto swJ = curElemVolVars[scvJ].saturation(0);
-                const auto dFluxN_dSnJ_pc = tj_up_n*MaterialLaw::dpc_dsw(paramsJ, swJ);
+                const auto dFluxN_dSnJ_pc = tj_up_n*fluidMatrixInteractionJ.dpc_dsw(swJ);
                 dI_dJ_inside[globalJ][conti0EqIdx+1][saturationIdx] -= dFluxN_dSnJ_pc;
                 dI_dJ_outside[globalJ][conti0EqIdx+1][saturationIdx] += dFluxN_dSnJ_pc;
             }
@@ -434,7 +468,6 @@ public:
                                        const ElementFluxVariablesCache& elemFluxVarsCache,
                                        const SubControlVolumeFace& scvf) const
     {
-        using MaterialLaw = typename Problem::SpatialParams::MaterialLaw;
         using AdvectionType = GetPropType<TypeTag, Properties::AdvectionType>;
 
         // evaluate the current wetting phase Darcy flux and resulting upwind weights
@@ -451,9 +484,16 @@ public:
         const auto& insideScv = fvGeometry.scv(insideScvIdx);
         const auto& insideVolVars = curElemVolVars[insideScvIdx];
         const auto& outsideVolVars = curElemVolVars[scvf.outsideScvIdx()];
-        const auto& insideMaterialParams = problem.spatialParams().materialLawParams(element,
-                                                                                     insideScv,
-                                                                                     elementSolution<FVElementGeometry>(insideVolVars.priVars()));
+
+        // old material law interface is deprecated: Replace this by
+        // const auto& insidefluidMatrixInteraction = problem.spatialParams().fluidMatrixInteraction(element,
+        //                                                                                           insideScv,
+        //                                                                                           elementSolution<FVElementGeometry>(insideVolVars.priVars()));
+        // after the release of 3.3, when the deprecated interface is no longer supported
+        const auto& insidefluidMatrixInteraction = Deprecated::makePcKrSw(Scalar{}, problem.spatialParams(),
+                                                                          element,
+                                                                          insideScv,
+                                                                          elementSolution<FVElementGeometry>(insideVolVars.priVars()));
 
         // some quantities to be reused (rho & mu are constant and thus equal for all cells)
         static const auto rho_w = insideVolVars.density(0);
@@ -470,9 +510,9 @@ public:
 
         // derivative w.r.t. to Sn is the negative of the one w.r.t. Sw
         const auto insideSw = insideVolVars.saturation(0);
-        const auto dKrw_dSn_inside = -1.0*MaterialLaw::dkrw_dsw(insideMaterialParams, insideSw);
-        const auto dKrn_dSn_inside = -1.0*MaterialLaw::dkrn_dsw(insideMaterialParams, insideSw);
-        const auto dpc_dSn_inside = -1.0*MaterialLaw::dpc_dsw(insideMaterialParams, insideSw);
+        const auto dKrw_dSn_inside = -1.0*insidefluidMatrixInteraction.dkrw_dsw(insideSw);
+        const auto dKrn_dSn_inside = -1.0*insidefluidMatrixInteraction.dkrn_dsw(insideSw);
+        const auto dpc_dSn_inside = -1.0*insidefluidMatrixInteraction.dpc_dsw(insideSw);
 
         const auto tij = elemFluxVarsCache[scvf].advectionTij();
         // partial derivative of the wetting phase flux w.r.t. p_w