diff --git a/Applications/ApplicationsLib/ProjectData.cpp b/Applications/ApplicationsLib/ProjectData.cpp
index 56a706bacdf7d9619f7f9ed81d6357137ffc765b..7e444c30464e8170f9382d13eb28562038c90525 100644
--- a/Applications/ApplicationsLib/ProjectData.cpp
+++ b/Applications/ApplicationsLib/ProjectData.cpp
@@ -892,7 +892,7 @@ void ProjectData::parseProcesses(BaseLib::ConfigTree const& processes_config,
                 ProcessLib::TwoPhaseFlowWithPP::createTwoPhaseFlowWithPPProcess(
                     name, *_mesh_vec[0], std::move(jacobian_assembler),
                     _process_variables, _parameters, integration_order,
-                    process_config, _curves);
+                    process_config, _curves, _media);
         }
         else
 #endif
diff --git a/MaterialLib/MPL/Properties/RelPermBrooksCorey.h b/MaterialLib/MPL/Properties/RelPermBrooksCorey.h
index a387ac129a5c78a0dbc15c0fc21f807d3d0d57d3..3b36b0da6d13ecbd5151768f21c0afc1451f1b79 100644
--- a/MaterialLib/MPL/Properties/RelPermBrooksCorey.h
+++ b/MaterialLib/MPL/Properties/RelPermBrooksCorey.h
@@ -31,7 +31,7 @@ class Component;
 class RelPermBrooksCorey final : public Property
 {
 private:
-    Medium* _medium;
+    Medium* _medium = nullptr;
     const double _residual_liquid_saturation;
     const double _residual_gas_saturation;
     const double _min_relative_permeability_liquid;
diff --git a/MaterialLib/MPL/Properties/SaturationBrooksCorey.cpp b/MaterialLib/MPL/Properties/SaturationBrooksCorey.cpp
index 54f87cacfc5165eda475ba99e602f24b52971a74..d1491e45c82bcf74f6401a31e945da067baea4ef 100644
--- a/MaterialLib/MPL/Properties/SaturationBrooksCorey.cpp
+++ b/MaterialLib/MPL/Properties/SaturationBrooksCorey.cpp
@@ -64,11 +64,10 @@ PropertyDataType SaturationBrooksCorey::dValue(
     const double s_L_res = _residual_liquid_saturation;
     const double s_L_max = 1.0 - _residual_gas_saturation;
     const double p_b = _entry_pressure;
-    const double p_cap = std::get<double>(
-        variable_array[static_cast<int>(Variable::capillary_pressure)]);
-
-    if (p_cap <= p_b)
-        return 0.;
+    const double p_cap = std::max(
+        p_b,
+        std::get<double>(
+            variable_array[static_cast<int>(Variable::capillary_pressure)]));
 
     auto const s_L = _medium->property(PropertyType::saturation)
                          .template value<double>(variable_array, pos, t);
diff --git a/ProcessLib/TwoPhaseFlowWithPP/CreateTwoPhaseFlowWithPPMaterialProperties.cpp b/ProcessLib/TwoPhaseFlowWithPP/CreateTwoPhaseFlowWithPPMaterialProperties.cpp
deleted file mode 100644
index 1c8f49f1c58626798c65faa618ef7082fe2ef56f..0000000000000000000000000000000000000000
--- a/ProcessLib/TwoPhaseFlowWithPP/CreateTwoPhaseFlowWithPPMaterialProperties.cpp
+++ /dev/null
@@ -1,140 +0,0 @@
-/**
- * \file
- * \copyright
- * Copyright (c) 2012-2019, OpenGeoSys Community (http://www.opengeosys.org)
- *            Distributed under a Modified BSD License.
- *              See accompanying file LICENSE.txt or
- *              http://www.opengeosys.org/project/license
- *
- */
-
-#include "CreateTwoPhaseFlowWithPPMaterialProperties.h"
-#include <logog/include/logog.hpp>
-#include "BaseLib/Algorithm.h"
-#include "MaterialLib/Fluid/FluidProperty.h"
-#include "MaterialLib/PorousMedium/Porosity/Porosity.h"
-#include "MaterialLib/PorousMedium/Storage/Storage.h"
-#include "MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/CapillaryPressureSaturation.h"
-#include "MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/CreateCapillaryPressureModel.h"
-#include "MaterialLib/PorousMedium/UnsaturatedProperty/RelativePermeability/CreateRelativePermeabilityModel.h"
-#include "MaterialLib/PorousMedium/UnsaturatedProperty/RelativePermeability/RelativePermeability.h"
-#include "MathLib/InterpolationAlgorithms/PiecewiseLinearInterpolation.h"
-#include "MeshLib/Mesh.h"
-#include "MeshLib/PropertyVector.h"
-#include "ParameterLib/Parameter.h"
-#include "ParameterLib/SpatialPosition.h"
-#include "TwoPhaseFlowWithPPMaterialProperties.h"
-
-namespace ProcessLib
-{
-namespace TwoPhaseFlowWithPP
-{
-std::unique_ptr<TwoPhaseFlowWithPPMaterialProperties>
-createTwoPhaseFlowWithPPMaterialProperties(
-    BaseLib::ConfigTree const& config,
-    MeshLib::PropertyVector<int> const* const material_ids,
-    std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters)
-{
-    DBUG("Reading material properties of two-phase flow process.");
-
-    //! \ogs_file_param{prj__processes__process__TWOPHASE_FLOW_PP__material_property__fluid}
-    auto const& fluid_config = config.getConfigSubtree("fluid");
-
-    // Get fluid properties
-    //! \ogs_file_param{prj__processes__process__TWOPHASE_FLOW_PP__material_property__liquid_density}
-    auto const& rho_conf = fluid_config.getConfigSubtree("liquid_density");
-    auto liquid_density = MaterialLib::Fluid::createFluidDensityModel(rho_conf);
-    //! \ogs_file_param{prj__processes__process__TWOPHASE_FLOW_PP__material_property__gas_density}
-    auto const& rho_gas_conf = fluid_config.getConfigSubtree("gas_density");
-    auto gas_density =
-        MaterialLib::Fluid::createFluidDensityModel(rho_gas_conf);
-    //! \ogs_file_param{prj__processes__process__TWOPHASE_FLOW_PP__material_property__liquid_viscosity}
-    auto const& mu_conf = fluid_config.getConfigSubtree("liquid_viscosity");
-    auto liquid_viscosity = MaterialLib::Fluid::createViscosityModel(mu_conf);
-    //! \ogs_file_param{prj__processes__process__TWOPHASE_FLOW_PP__material_property__gas_viscosity}
-    auto const& mu_gas_conf = fluid_config.getConfigSubtree("gas_viscosity");
-    auto gas_viscosity = MaterialLib::Fluid::createViscosityModel(mu_gas_conf);
-
-    // Get porous properties
-    std::vector<int> mat_ids;
-    std::vector<int> mat_krel_ids;
-    std::vector<std::unique_ptr<MaterialLib::PorousMedium::Permeability>>
-        intrinsic_permeability_models;
-    std::vector<std::unique_ptr<MaterialLib::PorousMedium::Porosity>>
-        porosity_models;
-    std::vector<std::unique_ptr<MaterialLib::PorousMedium::Storage>>
-        storage_models;
-    std::vector<
-        std::unique_ptr<MaterialLib::PorousMedium::CapillaryPressureSaturation>>
-        capillary_pressure_models;
-    std::vector<
-        std::unique_ptr<MaterialLib::PorousMedium::RelativePermeability>>
-        relative_permeability_models;
-
-    //! \ogs_file_param{prj__processes__process__TWOPHASE_FLOW_PP__material_property__porous_medium}
-    auto const& poro_config = config.getConfigSubtree("porous_medium");
-    //! \ogs_file_param{prj__processes__process__TWOPHASE_FLOW_PP__material_property__porous_medium__porous_medium}
-    for (auto const& conf : poro_config.getConfigSubtreeList("porous_medium"))
-    {
-        //! \ogs_file_attr{prj__processes__process__TWOPHASE_FLOW_PP__material_property__porous_medium__porous_medium__id}
-        auto const id = conf.getConfigAttributeOptional<int>("id");
-        mat_ids.push_back(*id);
-
-        //! \ogs_file_param{prj__processes__process__TWOPHASE_FLOW_PP__material_property__porous_medium__porous_medium__permeability}
-        auto const& permeability_conf = conf.getConfigSubtree("permeability");
-        intrinsic_permeability_models.emplace_back(
-            MaterialLib::PorousMedium::createPermeabilityModel(
-                permeability_conf, parameters));
-
-        //! \ogs_file_param{prj__processes__process__TWOPHASE_FLOW_PP__material_property__porous_medium__porous_medium__porosity}
-        auto const& porosity_conf = conf.getConfigSubtree("porosity");
-        auto n = MaterialLib::PorousMedium::createPorosityModel(porosity_conf,
-                                                                parameters);
-        porosity_models.emplace_back(std::move(n));
-
-        //! \ogs_file_param{prj__processes__process__TWOPHASE_FLOW_PP__material_property__porous_medium__porous_medium__storage}
-        auto const& storage_conf = conf.getConfigSubtree("storage");
-        auto beta = MaterialLib::PorousMedium::createStorageModel(storage_conf);
-        storage_models.emplace_back(std::move(beta));
-
-        auto const& capillary_pressure_conf =
-            //! \ogs_file_param{prj__processes__process__TWOPHASE_FLOW_PP__material_property__porous_medium__porous_medium__capillary_pressure}
-            conf.getConfigSubtree("capillary_pressure");
-        auto pc = MaterialLib::PorousMedium::createCapillaryPressureModel(
-            capillary_pressure_conf);
-        capillary_pressure_models.emplace_back(std::move(pc));
-
-        auto const& krel_config =
-            //! \ogs_file_param{prj__processes__process__TWOPHASE_FLOW_PP__material_property__porous_medium__porous_medium__relative_permeability}
-            conf.getConfigSubtree("relative_permeability");
-        for (
-            auto const& krel_conf :
-            //! \ogs_file_param{prj__processes__process__TWOPHASE_FLOW_PP__material_property__porous_medium__porous_medium__relative_permeability__relative_permeability}
-            krel_config.getConfigSubtreeList("relative_permeability"))
-        {
-            auto const krel_id =
-                //! \ogs_file_attr{prj__processes__process__TWOPHASE_FLOW_PP__material_property__porous_medium__porous_medium__relative_permeability__relative_permeability__id}
-                krel_conf.getConfigAttributeOptional<int>("id");
-            mat_krel_ids.push_back(*krel_id);
-            auto krel_n =
-                MaterialLib::PorousMedium::createRelativePermeabilityModel(
-                    krel_conf);
-            relative_permeability_models.emplace_back(std::move(krel_n));
-        }
-        BaseLib::reorderVector(relative_permeability_models, mat_krel_ids);
-    }
-
-    BaseLib::reorderVector(intrinsic_permeability_models, mat_ids);
-    BaseLib::reorderVector(porosity_models, mat_ids);
-    BaseLib::reorderVector(storage_models, mat_ids);
-
-    return std::make_unique<TwoPhaseFlowWithPPMaterialProperties>(
-        material_ids, std::move(liquid_density), std::move(liquid_viscosity),
-        std::move(gas_density), std::move(gas_viscosity),
-        std::move(intrinsic_permeability_models), std::move(porosity_models),
-        std::move(storage_models), std::move(capillary_pressure_models),
-        std::move(relative_permeability_models));
-}
-
-}  // namespace TwoPhaseFlowWithPP
-}  // namespace ProcessLib
diff --git a/ProcessLib/TwoPhaseFlowWithPP/CreateTwoPhaseFlowWithPPMaterialProperties.h b/ProcessLib/TwoPhaseFlowWithPP/CreateTwoPhaseFlowWithPPMaterialProperties.h
deleted file mode 100644
index 9c134c7748606dbc719d7482f1b87b018cf8d875..0000000000000000000000000000000000000000
--- a/ProcessLib/TwoPhaseFlowWithPP/CreateTwoPhaseFlowWithPPMaterialProperties.h
+++ /dev/null
@@ -1,32 +0,0 @@
-/**
- * \file
- * \copyright
- * Copyright (c) 2012-2019, OpenGeoSys Community (http://www.opengeosys.org)
- *            Distributed under a Modified BSD License.
- *              See accompanying file LICENSE.txt or
- *              http://www.opengeosys.org/project/license
- *
- */
-
-#pragma once
-
-
-#include "TwoPhaseFlowWithPPMaterialProperties.h"
-namespace BaseLib
-{
-class ConfigTree;
-}
-
-namespace ProcessLib
-{
-namespace TwoPhaseFlowWithPP
-{
-std::unique_ptr<TwoPhaseFlowWithPPMaterialProperties>
-createTwoPhaseFlowWithPPMaterialProperties(
-    BaseLib::ConfigTree const& config,
-    MeshLib::PropertyVector<int> const* const material_ids,
-    std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const&
-        parameters);
-
-}  // end namespace
-}  // namespace ProcessLib
diff --git a/ProcessLib/TwoPhaseFlowWithPP/CreateTwoPhaseFlowWithPPProcess.cpp b/ProcessLib/TwoPhaseFlowWithPP/CreateTwoPhaseFlowWithPPProcess.cpp
index 1e8028cf156fb31c6824f7c14f80838e615c3860..d9feaf9c0284f1c276d9dcfd9dc56d35be996e1b 100644
--- a/ProcessLib/TwoPhaseFlowWithPP/CreateTwoPhaseFlowWithPPProcess.cpp
+++ b/ProcessLib/TwoPhaseFlowWithPP/CreateTwoPhaseFlowWithPPProcess.cpp
@@ -10,14 +10,15 @@
 #include "CreateTwoPhaseFlowWithPPProcess.h"
 #include <cassert>
 
+#include "MaterialLib/MPL/CreateMaterialSpatialDistributionMap.h"
+#include "MaterialLib/MPL/MaterialSpatialDistributionMap.h"
+
 #include "MeshLib/MeshGenerators/MeshGenerator.h"
 #include "ParameterLib/ConstantParameter.h"
 #include "ParameterLib/Utils.h"
 #include "ProcessLib/Output/CreateSecondaryVariables.h"
 #include "ProcessLib/Utils/ProcessUtils.h"
 
-#include "CreateTwoPhaseFlowWithPPMaterialProperties.h"
-#include "TwoPhaseFlowWithPPMaterialProperties.h"
 #include "TwoPhaseFlowWithPPProcess.h"
 #include "TwoPhaseFlowWithPPProcessData.h"
 namespace ProcessLib
@@ -34,7 +35,8 @@ std::unique_ptr<Process> createTwoPhaseFlowWithPPProcess(
     BaseLib::ConfigTree const& config,
     std::map<std::string,
              std::unique_ptr<MathLib::PiecewiseLinearInterpolation>> const&
-        curves)
+        curves,
+    std::map<int, std::unique_ptr<MaterialPropertyLib::Medium>> const& media)
 {
     //! \ogs_file_param{prj__processes__process__type}
     config.checkConfigParameter("type", "TWOPHASE_FLOW_PP");
@@ -80,9 +82,6 @@ std::unique_ptr<Process> createTwoPhaseFlowWithPPProcess(
         //! \ogs_file_param_special{prj__processes__process__TWOPHASE_FLOW_PP__temperature}
         "temperature", parameters, 1, &mesh);
 
-    //! \ogs_file_param{prj__processes__process__TWOPHASE_FLOW_PP__material_property}
-    auto const& mat_config = config.getConfigSubtree("material_property");
-
     auto const material_ids = materialIDs(mesh);
     if (material_ids)
     {
@@ -92,18 +91,19 @@ std::unique_ptr<Process> createTwoPhaseFlowWithPPProcess(
     {
         INFO("The twophase flow is in homogeneous porous media.");
     }
-    std::unique_ptr<TwoPhaseFlowWithPPMaterialProperties> material =
-        createTwoPhaseFlowWithPPMaterialProperties(mat_config, material_ids,
-                                                   parameters);
 
-    TwoPhaseFlowWithPPProcessData process_data{
-        specific_body_force, has_gravity, mass_lumping, temperature, std::move(material)};
+    auto media_map =
+        MaterialPropertyLib::createMaterialSpatialDistributionMap(media, mesh);
+
+    TwoPhaseFlowWithPPProcessData process_data{specific_body_force, has_gravity,
+                                               mass_lumping, temperature,
+                                               std::move(media_map)};
 
     return std::make_unique<TwoPhaseFlowWithPPProcess>(
         std::move(name), mesh, std::move(jacobian_assembler), parameters,
         integration_order, std::move(process_variables),
         std::move(process_data), std::move(secondary_variables),
-        std::move(named_function_caller), mat_config, curves);
+        std::move(named_function_caller), curves);
 }
 
 }  // namespace TwoPhaseFlowWithPP
diff --git a/ProcessLib/TwoPhaseFlowWithPP/CreateTwoPhaseFlowWithPPProcess.h b/ProcessLib/TwoPhaseFlowWithPP/CreateTwoPhaseFlowWithPPProcess.h
index 97d0be571fa6bd1cadbe927e99e1328d7d492879..0be4ae6016f3e3b8cc0ed2eeeea60a88f7e07dc8 100644
--- a/ProcessLib/TwoPhaseFlowWithPP/CreateTwoPhaseFlowWithPPProcess.h
+++ b/ProcessLib/TwoPhaseFlowWithPP/CreateTwoPhaseFlowWithPPProcess.h
@@ -13,6 +13,11 @@
 #include <memory>
 #include "ProcessLib/Process.h"
 
+namespace MaterialPropertyLib
+{
+class Medium;
+}
+
 namespace ProcessLib
 {
 namespace TwoPhaseFlowWithPP
@@ -27,6 +32,7 @@ std::unique_ptr<Process> createTwoPhaseFlowWithPPProcess(
     BaseLib::ConfigTree const& config,
     std::map<std::string,
              std::unique_ptr<MathLib::PiecewiseLinearInterpolation>> const&
-        curves);
+        curves,
+    std::map<int, std::unique_ptr<MaterialPropertyLib::Medium>> const& media);
 }  // namespace TwoPhaseFlowWithPP
 }  // namespace ProcessLib
diff --git a/ProcessLib/TwoPhaseFlowWithPP/Tests.cmake b/ProcessLib/TwoPhaseFlowWithPP/Tests.cmake
index d8f1d334e4d5e4e93158b3063633fd54a9a047ed..325112f5e821a56ddd50e6b24e6cc36d1e6ff1a0 100644
--- a/ProcessLib/TwoPhaseFlowWithPP/Tests.cmake
+++ b/ProcessLib/TwoPhaseFlowWithPP/Tests.cmake
@@ -1,15 +1,5 @@
 AddTest(
-    NAME 2D_TwoPhase_PP_Lia_quad_1
-    PATH Parabolic/TwoPhaseFlowPP/Liakopoulos
-    EXECUTABLE ogs
-    EXECUTABLE_ARGS TwoPhase_Lia_quad_short.prj
-    TESTER vtkdiff
-    REQUIREMENTS NOT OGS_USE_MPI
-    DIFF_DATA
-    h2_Liako_20.vtu twophaseflow_pcs_0_ts_218_t_20.000000.vtu saturation saturation 1e-2 1e-4
-)
-AddTest(
-    NAME 2D_TwoPhase_PP_Lia_quad_2
+    NAME 2D_TwoPhase_PP_Lia_quad_short
     PATH Parabolic/TwoPhaseFlowPP/Liakopoulos
     EXECUTABLE ogs
     EXECUTABLE_ARGS TwoPhase_Lia_quad_short.prj
@@ -17,32 +7,24 @@ AddTest(
     REQUIREMENTS NOT OGS_USE_MPI
     RUNTIME 11
     DIFF_DATA
+    h2_Liako_20.vtu twophaseflow_pcs_0_ts_218_t_20.000000.vtu saturation saturation 1e-2 1e-4
     h2_Liako_20.vtu twophaseflow_pcs_0_ts_218_t_20.000000.vtu capillary_pressure capillary_pressure 20 1e-3
     h2_Liako_20.vtu twophaseflow_pcs_0_ts_218_t_20.000000.vtu gas_pressure gas_pressure 20 1e-3
 )
 AddTest(
-    NAME 2D_TwoPhase_PP_Lia_quad_1_lg
-    PATH Parabolic/TwoPhaseFlowPP/Liakopoulos
-    EXECUTABLE ogs
-    EXECUTABLE_ARGS TwoPhase_Lia_quad_large.prj
-    TESTER vtkdiff
-    REQUIREMENTS NOT OGS_USE_MPI
-    DIFF_DATA
-    h2_Liako_1198.vtu twophaseflow_pcs_0_ts_1198_t_1000.000000.vtu SATURATION1 saturation 1e-2 1e-3
-)
-AddTest(
-    NAME 2D_TwoPhase_PP_Lia_quad_2_lg
+    NAME 2D_TwoPhase_PP_Lia_quad_large
     PATH Parabolic/TwoPhaseFlowPP/Liakopoulos
     EXECUTABLE ogs
     EXECUTABLE_ARGS TwoPhase_Lia_quad_large.prj
     TESTER vtkdiff
     REQUIREMENTS NOT OGS_USE_MPI
     DIFF_DATA
-    h2_Liako_1198.vtu twophaseflow_pcs_0_ts_1198_t_1000.000000.vtu PRESSURE1 capillary_pressure 20 1e-2
-    h2_Liako_1198.vtu twophaseflow_pcs_0_ts_1198_t_1000.000000.vtu PRESSURE2 gas_pressure 20 1e-2
+    h2_Liako_1198.vtu twophaseflow_pcs_0_ts_1198_t_1000.000000.vtu saturation saturation 1e-2 1e-3
+    h2_Liako_1198.vtu twophaseflow_pcs_0_ts_1198_t_1000.000000.vtu capillary_pressure capillary_pressure 20 1e-2
+    h2_Liako_1198.vtu twophaseflow_pcs_0_ts_1198_t_1000.000000.vtu gas_pressure gas_pressure 20 1e-2
 )
 AddTest(
-    NAME 1D_TwoPhase_PP_mcwt_1
+    NAME 1D_TwoPhase_PP_mcwt
     PATH Parabolic/TwoPhaseFlowPP/McWhorter
     EXECUTABLE ogs
     EXECUTABLE_ARGS TwoPhase_mcwt_line.prj
@@ -50,15 +32,6 @@ AddTest(
     REQUIREMENTS NOT OGS_USE_MPI
     DIFF_DATA
     mcwt_1000.vtu twophaseflow_pcs_0_ts_519_t_1000.000000.vtu SATURATION1 saturation 1e-2 1e-4
-)
-AddTest(
-    NAME 1D_TwoPhase_PP_mcwt_2
-    PATH Parabolic/TwoPhaseFlowPP/McWhorter
-    EXECUTABLE ogs
-    EXECUTABLE_ARGS TwoPhase_mcwt_line.prj
-    TESTER vtkdiff
-    REQUIREMENTS NOT OGS_USE_MPI
-    DIFF_DATA
     mcwt_1000.vtu twophaseflow_pcs_0_ts_519_t_1000.000000.vtu PRESSURE1 capillary_pressure 10 1e-3
     mcwt_1000.vtu twophaseflow_pcs_0_ts_519_t_1000.000000.vtu PRESSURE2 gas_pressure 10 1e-3
 )
diff --git a/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPLocalAssembler-impl.h b/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPLocalAssembler-impl.h
index bd624b3a92942ce97db9ea3ef9ad53803b12356c..da2261143df3f1856682847825fc7b44481a9dec 100644
--- a/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPLocalAssembler-impl.h
+++ b/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPLocalAssembler-impl.h
@@ -9,30 +9,29 @@
  */
 
 /**
-* common nomenclature
-* --------------primary variable----------------------
-* pn_int_pt    pressure for nonwetting phase at each integration point
-* pc_int_pt    capillary pressure at each integration point
-* --------------secondary variable--------------------
-* temperature              capillary pressure
-* Sw wetting               phase saturation
-* dSw_dpc                  derivative of wetting phase saturation with respect
-* to capillary pressure
-* rho_nonwet               density of nonwetting phase
-* drhononwet_dpn           derivative of nonwetting phase density with respect
-*to nonwetting phase pressure
-* rho_wet                  density of wetting phase
-* k_rel_nonwet             relative permeability of nonwetting phase
-* mu_nonwet                viscosity of nonwetting phase
-* lambda_nonwet            mobility of nonwetting phase
-* k_rel_wet                relative permeability of wetting phase
-* mu_wet                   viscosity of wetting phase
-* lambda_wet               mobility of wetting phase
-*/
+ * common nomenclature
+ * --------------primary variable----------------------
+ * pn_int_pt    pressure for nonwetting phase at each integration point
+ * pc_int_pt    capillary pressure at each integration point
+ * --------------secondary variable--------------------
+ * Sw wetting               phase saturation
+ * dSw_dpc                  derivative of wetting phase saturation with respect
+ * to capillary pressure
+ * rho_nonwet               density of nonwetting phase
+ * drhononwet_dpn           derivative of nonwetting phase density with respect
+ *to nonwetting phase pressure
+ * rho_wet                  density of wetting phase
+ * k_rel_nonwet             relative permeability of nonwetting phase
+ * mu_nonwet                viscosity of nonwetting phase
+ * lambda_nonwet            mobility of nonwetting phase
+ * k_rel_wet                relative permeability of wetting phase
+ * mu_wet                   viscosity of wetting phase
+ * lambda_wet               mobility of wetting phase
+ */
 #pragma once
 
-#include "TwoPhaseFlowWithPPLocalAssembler.h"
-
+#include "MaterialLib/MPL/Medium.h"
+#include "MaterialLib/MPL/Utils/FormEigenTensor.h"
 #include "MathLib/InterpolationAlgorithms/PiecewiseLinearInterpolation.h"
 #include "NumLib/Function/Interpolation.h"
 #include "TwoPhaseFlowWithPPProcessData.h"
@@ -41,6 +40,8 @@ namespace ProcessLib
 {
 namespace TwoPhaseFlowWithPP
 {
+namespace MPL = MaterialPropertyLib;
+
 template <typename ShapeFunction, typename IntegrationMethod,
           unsigned GlobalDim>
 void TwoPhaseFlowWithPPLocalAssembler<
@@ -90,90 +91,87 @@ void TwoPhaseFlowWithPPLocalAssembler<
     auto Bl =
         local_b.template segment<cap_pressure_size>(cap_pressure_matrix_index);
 
+    auto const& medium = *_process_data.media_map->getMedium(_element.getID());
+    auto const& liquid_phase = medium.phase("AqueousLiquid");
+    auto const& gas_phase = medium.phase("Gas");
+
     unsigned const n_integration_points =
         _integration_method.getNumberOfPoints();
 
     ParameterLib::SpatialPosition pos;
     pos.setElementID(_element.getID());
-    const int material_id =
-        _process_data.material->getMaterialID(pos.getElementID().get());
-
-    const Eigen::MatrixXd& perm = _process_data.material->getPermeability(
-        material_id, t, pos, _element.getDimension());
-    assert(perm.rows() == _element.getDimension() || perm.rows() == 1);
-    GlobalDimMatrixType permeability = GlobalDimMatrixType::Zero(
-        _element.getDimension(), _element.getDimension());
-    if (perm.rows() == _element.getDimension())
-    {
-        permeability = perm;
-    }
-    else if (perm.rows() == 1)
-    {
-        permeability.diagonal().setConstant(perm(0, 0));
-    }
 
     for (unsigned ip = 0; ip < n_integration_points; ip++)
     {
+        pos.setIntegrationPoint(ip);
+        auto const& w = _ip_data[ip].integration_weight;
+        auto const& dNdx = _ip_data[ip].dNdx;
+        auto const& M = _ip_data[ip].massOperator;
+
         double pc_int_pt = 0.;
         double pn_int_pt = 0.;
         NumLib::shapeFunctionInterpolate(local_x, _ip_data[ip].N, pn_int_pt,
                                          pc_int_pt);
 
         _pressure_wet[ip] = pn_int_pt - pc_int_pt;
+        MPL::VariableArray variables;
 
-        const double temperature = _process_data.temperature(t, pos)[0];
-        double const rho_nonwet =
-            _process_data.material->getGasDensity(pn_int_pt, temperature);
-        double const rho_wet = _process_data.material->getLiquidDensity(
-            _pressure_wet[ip], temperature);
-
-        double const Sw = _process_data.material->getSaturation(
-            material_id, t, pos, pn_int_pt, temperature, pc_int_pt);
-
-        _saturation[ip] = Sw;
-
-        double dSw_dpc = _process_data.material->getSaturationDerivative(
-            material_id, t, pos, pn_int_pt, temperature, Sw);
-
-        double const porosity = _process_data.material->getPorosity(
-            material_id, t, pos, pn_int_pt, temperature, 0);
-
-        // Assemble M matrix
-        // nonwetting
-        double const drhononwet_dpn =
-            _process_data.material->getGasDensityDerivative(pn_int_pt,
-                                                            temperature);
-
-        Mgp.noalias() +=
-            porosity * (1 - Sw) * drhononwet_dpn * _ip_data[ip].massOperator;
-        Mgpc.noalias() +=
-            -porosity * rho_nonwet * dSw_dpc * _ip_data[ip].massOperator;
-
-        Mlpc.noalias() +=
-            porosity * dSw_dpc * rho_wet * _ip_data[ip].massOperator;
-
-        // nonwet
-        double const k_rel_nonwet =
-            _process_data.material->getNonwetRelativePermeability(
-                t, pos, pn_int_pt, temperature, Sw);
-        double const mu_nonwet =
-            _process_data.material->getGasViscosity(pn_int_pt, temperature);
-        double const lambda_nonwet = k_rel_nonwet / mu_nonwet;
-
-        // wet
-        double const k_rel_wet =
-            _process_data.material->getWetRelativePermeability(
-                t, pos, _pressure_wet[ip], temperature, Sw);
-        double const mu_wet = _process_data.material->getLiquidViscosity(
-            _pressure_wet[ip], temperature);
-        double const lambda_wet = k_rel_wet / mu_wet;
-
-        laplace_operator.noalias() = _ip_data[ip].dNdx.transpose() *
-                                     permeability * _ip_data[ip].dNdx *
-                                     _ip_data[ip].integration_weight;
+        variables[static_cast<int>(MPL::Variable::phase_pressure)] = pn_int_pt;
+        variables[static_cast<int>(MPL::Variable::capillary_pressure)] =
+            pc_int_pt;
+        variables[static_cast<int>(MPL::Variable::temperature)] =
+            _process_data.temperature(t, pos)[0];
 
-        Kgp.noalias() += rho_nonwet * lambda_nonwet * laplace_operator;
+        auto const rho_nonwet = gas_phase.property(MPL::PropertyType::density)
+                                    .template value<double>(variables, pos, t);
+
+        auto const rho_wet = liquid_phase.property(MPL::PropertyType::density)
+                                 .template value<double>(variables, pos, t);
+        auto& Sw = _saturation[ip];
+        Sw = medium.property(MPL::PropertyType::saturation)
+                 .template value<double>(variables, pos, t);
+
+        auto const dSw_dpc =
+            medium.property(MPL::PropertyType::saturation)
+                .template dValue<double>(
+                    variables, MPL::Variable::capillary_pressure, pos, t);
 
+        auto const porosity = medium.property(MPL::PropertyType::porosity)
+                                  .template value<double>(variables, pos, t);
+
+        auto const drhononwet_dpn =
+            gas_phase.property(MPL::PropertyType::density)
+                .template dValue<double>(variables,
+                                         MPL::Variable::phase_pressure, pos, t);
+
+        auto const k_rel =
+            medium.property(MPL::PropertyType::relative_permeability)
+                .template value<MPL::Pair>(variables, pos, t);
+
+        auto const k_rel_wet = k_rel[0];
+        auto const k_rel_nonwet = k_rel[1];
+
+        auto const mu_nonwet = gas_phase.property(MPL::PropertyType::viscosity)
+                                   .template value<double>(variables, pos, t);
+
+        auto const lambda_nonwet = k_rel_nonwet / mu_nonwet;
+
+        auto const mu_wet = liquid_phase.property(MPL::PropertyType::viscosity)
+                                .template value<double>(variables, pos, t);
+
+        auto const lambda_wet = k_rel_wet / mu_wet;
+
+        auto const K = MPL::formEigenTensor<GlobalDim>(
+            medium.property(MPL::PropertyType::permeability)
+                .template value<double>(variables, pos, t));
+
+        Mgp.noalias() += porosity * (1 - Sw) * drhononwet_dpn * M;
+        Mgpc.noalias() += -porosity * rho_nonwet * dSw_dpc * M;
+        Mlpc.noalias() += porosity * dSw_dpc * rho_wet * M;
+
+        laplace_operator.noalias() = dNdx.transpose() * K * dNdx * w;
+
+        Kgp.noalias() += rho_nonwet * lambda_nonwet * laplace_operator;
         Klp.noalias() += rho_wet * lambda_wet * laplace_operator;
         Klpc.noalias() += -rho_wet * lambda_wet * laplace_operator;
 
@@ -181,9 +179,7 @@ void TwoPhaseFlowWithPPLocalAssembler<
         {
             auto const& b = _process_data.specific_body_force;
 
-            NodalVectorType gravity_operator = _ip_data[ip].dNdx.transpose() *
-                                               permeability * b *
-                                               _ip_data[ip].integration_weight;
+            NodalVectorType gravity_operator = dNdx.transpose() * K * b * w;
             Bg.noalias() +=
                 rho_nonwet * rho_nonwet * lambda_nonwet * gravity_operator;
             Bl.noalias() += rho_wet * rho_wet * lambda_wet * gravity_operator;
diff --git a/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPLocalAssembler.h b/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPLocalAssembler.h
index b51b04452ac35d5fe556284bd8c3f7c0955b0b03..36c5397b053665b471df8805fc3b9fe58d4a9fc5 100644
--- a/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPLocalAssembler.h
+++ b/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPLocalAssembler.h
@@ -21,7 +21,6 @@
 #include "ProcessLib/LocalAssemblerTraits.h"
 #include "ProcessLib/Utils/InitShapeMatrices.h"
 
-#include "TwoPhaseFlowWithPPMaterialProperties.h"
 #include "TwoPhaseFlowWithPPProcessData.h"
 
 namespace ProcessLib
@@ -32,13 +31,12 @@ template <typename NodalRowVectorType, typename GlobalDimNodalMatrixType,
           typename NodalMatrixType>
 struct IntegrationPointData final
 {
-    explicit IntegrationPointData(
-        NodalRowVectorType N_, GlobalDimNodalMatrixType dNdx_,
-        TwoPhaseFlowWithPPMaterialProperties& material_property_,
-        double const& integration_weight_, NodalMatrixType const massOperator_)
+    explicit IntegrationPointData(NodalRowVectorType N_,
+                                  GlobalDimNodalMatrixType dNdx_,
+                                  double const& integration_weight_,
+                                  NodalMatrixType const massOperator_)
         : N(std::move(N_)),
           dNdx(std::move(dNdx_)),
-          mat_property(material_property_),
           integration_weight(integration_weight_),
           massOperator(massOperator_)
 
@@ -46,7 +44,6 @@ struct IntegrationPointData final
     }
     NodalRowVectorType const N;
     GlobalDimNodalMatrixType const dNdx;
-    TwoPhaseFlowWithPPMaterialProperties const& mat_property;
     const double integration_weight;
     NodalMatrixType const massOperator;
 
@@ -119,7 +116,7 @@ public:
         {
             auto const& sm = shape_matrices[ip];
             _ip_data.emplace_back(
-                sm.N, sm.dNdx, *_process_data.material,
+                sm.N, sm.dNdx,
                 sm.integralMeasure * sm.detJ *
                     _integration_method.getWeightedPoint(ip).getWeight(),
                 sm.N.transpose() * sm.N * sm.integralMeasure * sm.detJ *
diff --git a/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPMaterialProperties.cpp b/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPMaterialProperties.cpp
deleted file mode 100644
index 1560d348db0699bbbbc8bdd24e8e65499b7b6612..0000000000000000000000000000000000000000
--- a/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPMaterialProperties.cpp
+++ /dev/null
@@ -1,171 +0,0 @@
-/**
- * \file
- * \copyright
- * Copyright (c) 2012-2019, OpenGeoSys Community (http://www.opengeosys.org)
- *            Distributed under a Modified BSD License.
- *              See accompanying file LICENSE.txt or
- *              http://www.opengeosys.org/project/license
- *
- */
-
-#include "TwoPhaseFlowWithPPMaterialProperties.h"
-#include <logog/include/logog.hpp>
-#include <utility>
-#include "MaterialLib/Fluid/FluidProperty.h"
-#include "MaterialLib/PorousMedium/Porosity/Porosity.h"
-#include "MaterialLib/PorousMedium/Storage/Storage.h"
-#include "MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/CapillaryPressureSaturation.h"
-#include "MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/CreateCapillaryPressureModel.h"
-#include "MaterialLib/PorousMedium/UnsaturatedProperty/RelativePermeability/CreateRelativePermeabilityModel.h"
-#include "MaterialLib/PorousMedium/UnsaturatedProperty/RelativePermeability/RelativePermeability.h"
-#include "MathLib/InterpolationAlgorithms/PiecewiseLinearInterpolation.h"
-#include "MeshLib/Mesh.h"
-#include "MeshLib/PropertyVector.h"
-#include "ParameterLib/Parameter.h"
-#include "ParameterLib/SpatialPosition.h"
-namespace ProcessLib
-{
-namespace TwoPhaseFlowWithPP
-{
-TwoPhaseFlowWithPPMaterialProperties::TwoPhaseFlowWithPPMaterialProperties(
-    MeshLib::PropertyVector<int> const* const material_ids,
-    std::unique_ptr<MaterialLib::Fluid::FluidProperty>
-        liquid_density,
-    std::unique_ptr<MaterialLib::Fluid::FluidProperty>
-        liquid_viscosity,
-    std::unique_ptr<MaterialLib::Fluid::FluidProperty>
-        gas_density,
-    std::unique_ptr<MaterialLib::Fluid::FluidProperty>
-        gas_viscosity,
-    std::vector<std::unique_ptr<MaterialLib::PorousMedium::Permeability>>&&
-        intrinsic_permeability_models,
-    std::vector<std::unique_ptr<MaterialLib::PorousMedium::Porosity>>&&
-        porosity_models,
-    std::vector<std::unique_ptr<MaterialLib::PorousMedium::Storage>>&&
-        storage_models,
-    std::vector<std::unique_ptr<
-        MaterialLib::PorousMedium::CapillaryPressureSaturation>>&&
-        capillary_pressure_models,
-    std::vector<
-        std::unique_ptr<MaterialLib::PorousMedium::RelativePermeability>>&&
-        relative_permeability_models)
-    : _liquid_density(std::move(liquid_density)),
-      _liquid_viscosity(std::move(liquid_viscosity)),
-      _gas_density(std::move(gas_density)),
-      _gas_viscosity(std::move(gas_viscosity)),
-      _material_ids(material_ids),
-      _intrinsic_permeability_models(std::move(intrinsic_permeability_models)),
-      _porosity_models(std::move(porosity_models)),
-      _storage_models(std::move(storage_models)),
-      _capillary_pressure_models(std::move(capillary_pressure_models)),
-      _relative_permeability_models(std::move(relative_permeability_models))
-{
-    DBUG("Create material properties for Two-Phase flow with PP model.");
-}
-
-int TwoPhaseFlowWithPPMaterialProperties::getMaterialID(
-    const std::size_t element_id)
-{
-    if (!_material_ids)
-    {
-        return 0;
-    }
-
-    assert(element_id < _material_ids->size());
-    return (*_material_ids)[element_id];
-}
-
-double TwoPhaseFlowWithPPMaterialProperties::getLiquidDensity(
-    const double p, const double T) const
-{
-    ArrayType vars;
-    vars[static_cast<int>(MaterialLib::Fluid::PropertyVariableType::T)] = T;
-    vars[static_cast<int>(MaterialLib::Fluid::PropertyVariableType::p)] = p;
-    return _liquid_density->getValue(vars);
-}
-
-double TwoPhaseFlowWithPPMaterialProperties::getGasDensity(const double p,
-                                                           const double T) const
-{
-    ArrayType vars;
-    vars[static_cast<int>(MaterialLib::Fluid::PropertyVariableType::T)] = T;
-    vars[static_cast<int>(MaterialLib::Fluid::PropertyVariableType::p)] = p;
-    return _gas_density->getValue(vars);
-}
-
-double TwoPhaseFlowWithPPMaterialProperties::getGasDensityDerivative(
-    const double p, const double T) const
-{
-    ArrayType vars;
-    vars[static_cast<int>(MaterialLib::Fluid::PropertyVariableType::T)] = T;
-    vars[static_cast<int>(MaterialLib::Fluid::PropertyVariableType::p)] = p;
-
-    return _gas_density->getdValue(vars,
-                                   MaterialLib::Fluid::PropertyVariableType::p);
-}
-double TwoPhaseFlowWithPPMaterialProperties::getLiquidViscosity(
-    const double p, const double T) const
-{
-    ArrayType vars;
-    vars[static_cast<int>(MaterialLib::Fluid::PropertyVariableType::T)] = T;
-    vars[static_cast<int>(MaterialLib::Fluid::PropertyVariableType::p)] = p;
-    return _liquid_viscosity->getValue(vars);
-}
-
-double TwoPhaseFlowWithPPMaterialProperties::getGasViscosity(
-    const double p, const double T) const
-{
-    ArrayType vars;
-    vars[static_cast<int>(MaterialLib::Fluid::PropertyVariableType::T)] = T;
-    vars[static_cast<int>(MaterialLib::Fluid::PropertyVariableType::p)] = p;
-    return _gas_viscosity->getValue(vars);
-}
-
-Eigen::MatrixXd TwoPhaseFlowWithPPMaterialProperties::getPermeability(
-    const int material_id, const double t,
-    const ParameterLib::SpatialPosition& pos, const int /*dim*/) const
-{
-    return _intrinsic_permeability_models[material_id]->getValue(t, pos, 0.0,
-                                                                 0.0);
-}
-
-double TwoPhaseFlowWithPPMaterialProperties::getPorosity(
-    const int material_id, const double t,
-    const ParameterLib::SpatialPosition& pos, const double /*p*/,
-    const double T, const double porosity_variable) const
-{
-    return _porosity_models[material_id]->getValue(t, pos, porosity_variable, T);
-}
-
-double TwoPhaseFlowWithPPMaterialProperties::getNonwetRelativePermeability(
-    const double /*t*/, const ParameterLib::SpatialPosition& /*pos*/,
-    const double /*p*/, const double /*T*/, const double saturation) const
-{
-    return _relative_permeability_models[0]->getValue(saturation);
-}
-
-double TwoPhaseFlowWithPPMaterialProperties::getWetRelativePermeability(
-    const double /*t*/, const ParameterLib::SpatialPosition& /*pos*/,
-    const double /*p*/, const double /*T*/, const double saturation) const
-{
-    return _relative_permeability_models[1]->getValue(saturation);
-}
-
-double TwoPhaseFlowWithPPMaterialProperties::getSaturation(
-    const int material_id, const double /*t*/,
-    const ParameterLib::SpatialPosition& /*pos*/, const double /*p*/,
-    const double /*T*/, const double pc) const
-{
-    return _capillary_pressure_models[material_id]->getSaturation(pc);
-}
-double TwoPhaseFlowWithPPMaterialProperties::getSaturationDerivative(
-    const int material_id, const double /*t*/,
-    const ParameterLib::SpatialPosition& /*pos*/, const double /*p*/,
-    const double /*T*/, const double saturation) const
-{
-    const double dpcdsw =
-        _capillary_pressure_models[material_id]->getdPcdS(saturation);
-    return 1 / dpcdsw;
-}
-}  // namespace TwoPhaseFlowWithPP
-}  // namespace ProcessLib
diff --git a/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPMaterialProperties.h b/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPMaterialProperties.h
deleted file mode 100644
index 61d21a5b6d8c10883aef6a5ff9e84729d93bb92b..0000000000000000000000000000000000000000
--- a/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPMaterialProperties.h
+++ /dev/null
@@ -1,125 +0,0 @@
-/**
- * \file
- * \copyright
- * Copyright (c) 2012-2019, OpenGeoSys Community (http://www.opengeosys.org)
- *            Distributed under a Modified BSD License.
- *              See accompanying file LICENSE.txt or
- *              http://www.opengeosys.org/project/license
- *
- */
-
-#pragma once
-
-#include <memory>
-#include <vector>
-#include "MaterialLib/Fluid/FluidPropertyHeaders.h"
-#include "MaterialLib/PhysicalConstant.h"
-#include "MaterialLib/PorousMedium/Porosity/Porosity.h"
-#include "MaterialLib/PorousMedium/PorousPropertyHeaders.h"
-#include "MaterialLib/PorousMedium/Storage/Storage.h"
-#include "MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/CapillaryPressureSaturation.h"
-#include "MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/CreateCapillaryPressureModel.h"
-#include "MaterialLib/PorousMedium/UnsaturatedProperty/RelativePermeability/CreateRelativePermeabilityModel.h"
-#include "MaterialLib/PorousMedium/UnsaturatedProperty/RelativePermeability/RelativePermeability.h"
-
-namespace MeshLib
-{
-template <typename PROP_VAL_TYPE>
-class PropertyVector;
-}
-
-namespace ProcessLib::TwoPhaseFlowWithPP
-{
-    /** This class has a collection of material properties for two-phase flow with PP model
-    *  and it makes description of the material properties for two-phase condition,
-    *  i.e. the gas/liquid density and viscosity models, respectively,
-    *  the relative permeability models with respect to two phases,
-    *  the capillary pressure-saturation relationships.
-    *  It generally provides the computation of the PDE coefficients for two-phase flow.
-    */
-
-class TwoPhaseFlowWithPPMaterialProperties
-{
-public:
-    using ArrayType = MaterialLib::Fluid::FluidProperty::ArrayType;
-
-    TwoPhaseFlowWithPPMaterialProperties(
-        MeshLib::PropertyVector<int> const* const material_ids,
-        std::unique_ptr<MaterialLib::Fluid::FluidProperty>
-            liquid_density,
-        std::unique_ptr<MaterialLib::Fluid::FluidProperty>
-            liquid_viscosity,
-        std::unique_ptr<MaterialLib::Fluid::FluidProperty>
-            gas_density,
-        std::unique_ptr<MaterialLib::Fluid::FluidProperty>
-            gas_viscosity,
-        std::vector<std::unique_ptr<MaterialLib::PorousMedium::Permeability>>&&
-            intrinsic_permeability_models,
-        std::vector<std::unique_ptr<MaterialLib::PorousMedium::Porosity>>&&
-            porosity_models,
-        std::vector<std::unique_ptr<MaterialLib::PorousMedium::Storage>>&&
-            storage_models,
-        std::vector<std::unique_ptr<
-            MaterialLib::PorousMedium::CapillaryPressureSaturation>>&&
-            capillary_pressure_models,
-        std::vector<
-            std::unique_ptr<MaterialLib::PorousMedium::RelativePermeability>>&&
-            relative_permeability_models);
-
-    int getMaterialID(const std::size_t element_id);
-
-    Eigen::MatrixXd getPermeability(const int material_id,
-                                    const double t,
-                                    const ParameterLib::SpatialPosition& pos,
-                                    const int dim) const;
-
-    double getPorosity(const int material_id, const double t,
-                       const ParameterLib::SpatialPosition& pos, const double p,
-                       const double T, const double porosity_variable) const;
-
-    double getNonwetRelativePermeability(
-        const double t, const ParameterLib::SpatialPosition& pos,
-        const double p, const double T, const double saturation) const;
-    double getWetRelativePermeability(const double t,
-                                      const ParameterLib::SpatialPosition& pos,
-                                      const double p, const double T,
-                                      const double saturation) const;
-    double getSaturation(const int material_id, const double t,
-                         const ParameterLib::SpatialPosition& pos,
-                         const double p, const double T, const double pc) const;
-    double getSaturationDerivative(const int material_id, const double t,
-                                   const ParameterLib::SpatialPosition& pos,
-                                   const double p, const double T,
-                                   const double saturation) const;
-    double getLiquidDensity(const double p, const double T) const;
-    double getGasDensity(const double p, const double T) const;
-    double getGasViscosity(const double p, const double T) const;
-    double getLiquidViscosity(const double p, const double T) const;
-    double getGasDensityDerivative(double const p, double const T) const;
-
-protected:
-    std::unique_ptr<MaterialLib::Fluid::FluidProperty> _liquid_density;
-    std::unique_ptr<MaterialLib::Fluid::FluidProperty> _liquid_viscosity;
-    std::unique_ptr<MaterialLib::Fluid::FluidProperty> _gas_density;
-    std::unique_ptr<MaterialLib::Fluid::FluidProperty> _gas_viscosity;
-
-    /** Use two phase models for different material zones.
-    *  Material IDs must be given as mesh element properties.
-    */
-    MeshLib::PropertyVector<int> const* const _material_ids;
-
-    std::vector<std::unique_ptr<MaterialLib::PorousMedium::Permeability>>
-        _intrinsic_permeability_models;
-    std::vector<std::unique_ptr<MaterialLib::PorousMedium::Porosity>>
-        _porosity_models;
-    std::vector<std::unique_ptr<MaterialLib::PorousMedium::Storage>>
-        _storage_models;
-    std::vector<
-        std::unique_ptr<MaterialLib::PorousMedium::CapillaryPressureSaturation>>
-        _capillary_pressure_models;
-    std::vector<
-        std::unique_ptr<MaterialLib::PorousMedium::RelativePermeability>>
-        _relative_permeability_models;
-};
-
-}  // namespace ProcessLib::TwoPhaseFlowWithPP
diff --git a/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPProcess.cpp b/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPProcess.cpp
index bf5bf0f4c20dae8cd40a4c3520c65250832e977f..c1cf2c66500b4922b46e7673923c5876cf866a61 100644
--- a/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPProcess.cpp
+++ b/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPProcess.cpp
@@ -35,7 +35,6 @@ TwoPhaseFlowWithPPProcess::TwoPhaseFlowWithPPProcess(
     TwoPhaseFlowWithPPProcessData&& process_data,
     SecondaryVariableCollection&& secondary_variables,
     NumLib::NamedFunctionCaller&& named_function_caller,
-    BaseLib::ConfigTree const& /*config*/,
     std::map<std::string,
              std::unique_ptr<MathLib::PiecewiseLinearInterpolation>> const&
     /*curves*/)
diff --git a/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPProcess.h b/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPProcess.h
index 83a7355615bd7d791c52876e7a4928fe1c6587ed..c020b37dfdea76930934c0690f37d5ee905f7d8a 100644
--- a/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPProcess.h
+++ b/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPProcess.h
@@ -21,7 +21,7 @@ class Element;
 class Mesh;
 template <typename PROP_VAL_TYPE>
 class PropertyVector;
-}
+}  // namespace MeshLib
 
 namespace ProcessLib
 {
@@ -48,7 +48,6 @@ public:
         TwoPhaseFlowWithPPProcessData&& process_data,
         SecondaryVariableCollection&& secondary_variables,
         NumLib::NamedFunctionCaller&& named_function_caller,
-        BaseLib::ConfigTree const& config,
         std::map<std::string,
                  std::unique_ptr<MathLib::PiecewiseLinearInterpolation>> const&
             curves);
diff --git a/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPProcessData.h b/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPProcessData.h
index 3aaabe714b2477c6a5a6466257a8f51fe1ec4549..525971480b43b59b0a8d764aacabcb3126c1aec2 100644
--- a/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPProcessData.h
+++ b/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPProcessData.h
@@ -9,7 +9,7 @@
  */
 
 #pragma once
-#include "TwoPhaseFlowWithPPMaterialProperties.h"
+#include "MaterialLib/MPL/MaterialSpatialDistributionMap.h"
 
 namespace ProcessLib
 {
@@ -30,7 +30,8 @@ struct TwoPhaseFlowWithPPProcessData
     //! Enables lumping of the mass matrix.
     bool const has_mass_lumping;
     ParameterLib::Parameter<double> const& temperature;
-    std::unique_ptr<TwoPhaseFlowWithPPMaterialProperties> material;
+    std::unique_ptr<MaterialPropertyLib::MaterialSpatialDistributionMap>
+        media_map;
 };
 
 }  // namespace TwoPhaseFlowWithPP
diff --git a/Tests/Data/Parabolic/TwoPhaseFlowPP/Liakopoulos/TwoPhase_Lia_quad_large.prj b/Tests/Data/Parabolic/TwoPhaseFlowPP/Liakopoulos/TwoPhase_Lia_quad_large.prj
index facc2425bac8cc5dcd1461651ab4a65e03da3609..d4851385f5d73499e3bd05ac3330aab7b3a33d80 100644
--- a/Tests/Data/Parabolic/TwoPhaseFlowPP/Liakopoulos/TwoPhase_Lia_quad_large.prj
+++ b/Tests/Data/Parabolic/TwoPhaseFlowPP/Liakopoulos/TwoPhase_Lia_quad_large.prj
@@ -16,69 +16,6 @@
                 <gas_pressure>gas_pressure</gas_pressure>
                 <capillary_pressure>capillary_pressure</capillary_pressure>
             </process_variables>
-            <material_property>
-                <fluid>
-                    <liquid_density>
-                        <type>Constant</type>
-                        <value> 1.e3 </value>
-                    </liquid_density>
-                    <gas_density>
-                        <type>IdealGasLaw</type>
-                        <molar_mass> 0.02896 </molar_mass>
-                    </gas_density>
-                    <liquid_viscosity>
-                        <type>Constant</type>
-                        <value> 1.e-3 </value>
-                    </liquid_viscosity>
-                    <gas_viscosity>
-                        <type>Constant</type>
-                        <value> 1.8e-5 </value>
-                    </gas_viscosity>
-                </fluid>
-                <porous_medium>
-                    <porous_medium id="0">
-                        <permeability>
-                            <permeability_tensor_entries>kappa1</permeability_tensor_entries>
-                            <type>Constant</type>
-                        </permeability>
-                        <porosity>
-                            <type>Constant</type>
-                            <porosity_parameter>constant_porosity_parameter</porosity_parameter>
-                        </porosity>
-                        <storage>
-                            <type>Constant</type>
-                            <value> 0.0 </value>
-                        </storage>
-                        <capillary_pressure>
-                            <type>Curve</type>
-                            <curve>
-                                <coords>0.900000007    0.910000009    0.920000003    0.929821599    0.939999996    0.950000006    0.960000005    0.97    0.980000004    0.982500005    0.984999999    0.9875    0.989999997    0.9925    0.995000002    0.997499999    0.99962099    1
-                                </coords>
-                                <values>9938.064    9516.018    9065.393    8589.271    8052.432    7469.886    6813.948    6052.562    5121.663    4847.584    4549.372    4220.252    3849.668    3419.508    2893.579    2174.942    1000    0
-                                </values>
-                            </curve>
-                        </capillary_pressure>
-                        <relative_permeability>
-                            <relative_permeability id="0">
-                                <type>NonWettingPhaseBrooksCoreyOilGas</type>
-                                <sr>  0.2 </sr>
-                                <smax> 1.0</smax>
-                                <m> 3 </m>
-                                <krel_min> 1.e-4 </krel_min>
-                            </relative_permeability>
-                            <relative_permeability id="1">
-                                <type>Curve</type>
-                                <curve>
-                                    <coords>0.575    0.6    0.625    0.65    0.675    0.7    0.725    0.75    0.775    0.8    0.825    0.85    0.875    0.9    0.925    0.95    0.975    1
-                                    </coords>
-                                    <values>0    0.12693    0.18214    0.2373    0.29242    0.34748    0.40248    0.45743    0.51231    0.56711    0.62184    0.67646    0.73098    0.78536    0.83958    0.89358    0.94723    1
-                                    </values>
-                                </curve>
-                            </relative_permeability>
-                        </relative_permeability>
-                    </porous_medium>
-                </porous_medium>
-            </material_property>
             <secondary_variables>
                 <secondary_variable type="static" internal_name="saturation" output_name="saturation"/>
             </secondary_variables>
@@ -87,6 +24,71 @@
             <temperature> temp </temperature>
         </process>
     </processes>
+    <media>
+      <medium id="0">
+         <phases>
+            <phase>
+               <type>AqueousLiquid</type>
+               <properties>
+                  <property>
+                     <name>density</name>
+                     <type>Constant</type>
+                     <value>1.0e3</value>
+                  </property>
+                  <property>
+                     <name>viscosity</name>
+                     <type>Constant</type>
+                     <value>1e-3</value>
+                  </property>
+              </properties>
+            </phase>
+            <phase>
+               <type>Gas</type>
+               <properties>
+                  <property>
+                     <name>density</name>
+                     <type>IdealGasLaw</type>
+                  </property>
+                  <property>
+                     <name>viscosity</name>
+                     <type>Constant</type>
+                     <value>1.8e-5</value>
+                  </property>
+                  <property>
+                     <name>molar_mass</name>
+                     <type>Constant</type>
+                     <value>0.02896</value>
+                  </property>
+               </properties>
+            </phase>
+         </phases>
+         <properties>
+
+            <property>
+               <name>permeability</name>
+               <type>Constant</type>
+               <value>4.5e-13</value>
+            </property>
+
+            <property>
+               <name>saturation</name>
+               <type>SaturationLiakopoulos</type>
+            </property>
+
+            <property>
+               <name>relative_permeability</name>
+               <type>RelPermLiakopoulos</type>
+            </property>
+
+            <property>
+               <name>porosity</name>
+               <type>Constant</type>
+               <value>0.2975</value>
+            </property>
+         </properties>
+      </medium>
+   </media>
+
     <parameters>
         <parameter>
             <name>pc_Dirichlet_bottom</name>
diff --git a/Tests/Data/Parabolic/TwoPhaseFlowPP/Liakopoulos/TwoPhase_Lia_quad_short.prj b/Tests/Data/Parabolic/TwoPhaseFlowPP/Liakopoulos/TwoPhase_Lia_quad_short.prj
index 013064cdd0ffab0d2fcab06ef59583e2ef829ad2..d4e4a13982cf81e0a98f7fe2ec782aabf39674fe 100644
--- a/Tests/Data/Parabolic/TwoPhaseFlowPP/Liakopoulos/TwoPhase_Lia_quad_short.prj
+++ b/Tests/Data/Parabolic/TwoPhaseFlowPP/Liakopoulos/TwoPhase_Lia_quad_short.prj
@@ -16,69 +16,6 @@
                 <gas_pressure>gas_pressure</gas_pressure>
                 <capillary_pressure>capillary_pressure</capillary_pressure>
             </process_variables>
-            <material_property>
-                <fluid>
-                    <liquid_density>
-                        <type>Constant</type>
-                        <value> 1.e3 </value>
-                    </liquid_density>
-                    <gas_density>
-                        <type>IdealGasLaw</type>
-                        <molar_mass> 0.02896 </molar_mass>
-                    </gas_density>
-                    <liquid_viscosity>
-                        <type>Constant</type>
-                        <value> 1.e-3 </value>
-                    </liquid_viscosity>
-                    <gas_viscosity>
-                        <type>Constant</type>
-                        <value> 1.8e-5 </value>
-                    </gas_viscosity>
-                </fluid>
-                <porous_medium>
-                    <porous_medium id="0">
-                        <permeability>
-                            <permeability_tensor_entries>kappa1</permeability_tensor_entries>
-                            <type>Constant</type>
-                        </permeability>
-                        <porosity>
-                            <type>Constant</type>
-                            <porosity_parameter>constant_porosity_parameter</porosity_parameter>
-                        </porosity>
-                        <storage>
-                            <type>Constant</type>
-                            <value> 0.0 </value>
-                        </storage>
-                        <capillary_pressure>
-                            <type>Curve</type>
-                            <curve>
-                                <coords>0.900000007    0.910000009    0.920000003    0.929821599    0.939999996    0.950000006    0.960000005    0.97    0.980000004    0.982500005    0.984999999    0.9875    0.989999997    0.9925    0.995000002    0.997499999    0.99962099    1
-                                </coords>
-                                <values>9938.064    9516.018    9065.393    8589.271    8052.432    7469.886    6813.948    6052.562    5121.663    4847.584    4549.372    4220.252    3849.668    3419.508    2893.579    2174.942    1000    0
-                                </values>
-                            </curve>
-                        </capillary_pressure>
-                        <relative_permeability>
-                            <relative_permeability id="0">
-                                <type>NonWettingPhaseBrooksCoreyOilGas</type>
-                                <sr>  0.2 </sr>
-                                <smax> 1.0</smax>
-                                <m> 3 </m>
-                                <krel_min> 1.e-4 </krel_min>
-                            </relative_permeability>
-                            <relative_permeability id="1">
-                                <type>Curve</type>
-                                <curve>
-                                    <coords>0.575    0.6    0.625    0.65    0.675    0.7    0.725    0.75    0.775    0.8    0.825    0.85    0.875    0.9    0.925    0.95    0.975    1
-                                    </coords>
-                                    <values>0    0.12693    0.18214    0.2373    0.29242    0.34748    0.40248    0.45743    0.51231    0.56711    0.62184    0.67646    0.73098    0.78536    0.83958    0.89358    0.94723    1
-                                    </values>
-                                </curve>
-                            </relative_permeability>
-                        </relative_permeability>
-                    </porous_medium>
-                </porous_medium>
-            </material_property>
             <secondary_variables>
                 <secondary_variable type="static" internal_name="saturation" output_name="saturation"/>
             </secondary_variables>
@@ -87,6 +24,71 @@
             <temperature> temp </temperature>
         </process>
     </processes>
+    <media>
+      <medium id="0">
+         <phases>
+            <phase>
+               <type>AqueousLiquid</type>
+               <properties>
+                  <property>
+                     <name>density</name>
+                     <type>Constant</type>
+                     <value>1.0e3</value>
+                  </property>
+                  <property>
+                     <name>viscosity</name>
+                     <type>Constant</type>
+                     <value>1e-3</value>
+                  </property>
+              </properties>
+            </phase>
+            <phase>
+               <type>Gas</type>
+               <properties>
+                  <property>
+                     <name>density</name>
+                     <type>IdealGasLaw</type>
+                  </property>
+                  <property>
+                     <name>viscosity</name>
+                     <type>Constant</type>
+                     <value>1.8e-5</value>
+                  </property>
+                  <property>
+                     <name>molar_mass</name>
+                     <type>Constant</type>
+                     <value>0.02896</value>
+                  </property>
+               </properties>
+            </phase>
+         </phases>
+         <properties>
+
+            <property>
+               <name>permeability</name>
+               <type>Constant</type>
+               <value>4.5e-13</value>
+            </property>
+
+            <property>
+               <name>saturation</name>
+               <type>SaturationLiakopoulos</type>
+            </property>
+
+            <property>
+               <name>relative_permeability</name>
+               <type>RelPermLiakopoulos</type>
+            </property>
+
+            <property>
+               <name>porosity</name>
+               <type>Constant</type>
+               <value>0.2975</value>
+            </property>
+
+         </properties>
+      </medium>
+   </media>
     <parameters>
         <parameter>
             <name>pc_Dirichlet_bottom</name>
diff --git a/Tests/Data/Parabolic/TwoPhaseFlowPP/Liakopoulos/h2_Liako_1198.vtu b/Tests/Data/Parabolic/TwoPhaseFlowPP/Liakopoulos/h2_Liako_1198.vtu
index c814fb18c2cc96976eab57c41709809ce83f3134..f62802d06dd7e05bf0505510eb5412e9eeb92148 100644
--- a/Tests/Data/Parabolic/TwoPhaseFlowPP/Liakopoulos/h2_Liako_1198.vtu
+++ b/Tests/Data/Parabolic/TwoPhaseFlowPP/Liakopoulos/h2_Liako_1198.vtu
@@ -1,3 +1,3 @@
 version https://git-lfs.github.com/spec/v1
-oid sha256:032baddc6cebf9d143d6c3f1c5cee4e1a4066e6726ddb079d8050b1248a0c3dc
-size 28294
+oid sha256:e88e8bd76e0132e45d4176abf91bd316802ca640e056a44db49c2daa332519f3
+size 5251
diff --git a/Tests/Data/Parabolic/TwoPhaseFlowPP/Liakopoulos/h2_Liako_20.vtu b/Tests/Data/Parabolic/TwoPhaseFlowPP/Liakopoulos/h2_Liako_20.vtu
index 58de12352af9ab53e6d3431be7bc23243d3c0397..6e1aa1a98a223b08cc9530859b3cebcafe3445cf 100644
--- a/Tests/Data/Parabolic/TwoPhaseFlowPP/Liakopoulos/h2_Liako_20.vtu
+++ b/Tests/Data/Parabolic/TwoPhaseFlowPP/Liakopoulos/h2_Liako_20.vtu
@@ -1,3 +1,3 @@
 version https://git-lfs.github.com/spec/v1
-oid sha256:47e2c373144b143886cb71ca6e4d35e1ce9ca625bdaeabe7302716baf73a46b7
-size 5087
+oid sha256:d1378a7a33c3734e6aea4e126d1aa1d44c270e8641abbf473e493265cab4ece1
+size 5263
diff --git a/Tests/Data/Parabolic/TwoPhaseFlowPP/McWhorter/TwoPhase_mcwt_line.prj b/Tests/Data/Parabolic/TwoPhaseFlowPP/McWhorter/TwoPhase_mcwt_line.prj
index e5aceb0f24fade95b86a086101669475b67412e1..9e05c08e158a4e023b3b938aec3eec6d247fc5cd 100644
--- a/Tests/Data/Parabolic/TwoPhaseFlowPP/McWhorter/TwoPhase_mcwt_line.prj
+++ b/Tests/Data/Parabolic/TwoPhaseFlowPP/McWhorter/TwoPhase_mcwt_line.prj
@@ -16,66 +16,6 @@
                 <capillary_pressure>capillary_pressure</capillary_pressure>
                 <gas_pressure>gas_pressure</gas_pressure>
             </process_variables>
-            <material_property>
-                <fluid>
-                    <liquid_density>
-                        <type>Constant</type>
-                        <value> 1.e3 </value>
-                    </liquid_density>
-                    <gas_density>
-                        <type>Constant</type>
-                        <value> 1.e3  </value>
-                    </gas_density>
-                    <liquid_viscosity>
-                        <type>Constant</type>
-                        <value> 1.e-3 </value>
-                    </liquid_viscosity>
-                    <gas_viscosity>
-                        <type>Constant</type>
-                        <value> 1.e-3 </value>
-                    </gas_viscosity>
-                </fluid>
-                <porous_medium>
-                    <porous_medium id="0">
-                        <permeability>
-                            <permeability_tensor_entries>kappa1</permeability_tensor_entries>
-                            <type>Constant</type>
-                        </permeability>
-                        <porosity>
-                            <type>Constant</type>
-                            <porosity_parameter>constant_porosity_parameter</porosity_parameter>
-                        </porosity>
-                        <storage>
-                            <type>Constant</type>
-                            <value> 0.0 </value>
-                        </storage>
-                        <capillary_pressure>
-                            <type>BrooksCorey</type>
-                            <pd> 5000 </pd>
-                            <sr> 0.0 </sr>
-                            <smax> 1. </smax>
-                            <m> 2 </m>
-                            <pc_max> 2e+7 </pc_max>
-                        </capillary_pressure>
-                        <relative_permeability>
-                            <relative_permeability id="0">
-                                <type>NonWettingPhaseBrooksCoreyOilGas</type>
-                                <sr>  0.0 </sr>
-                                <smax> 1.0</smax>
-                                <m> 2 </m>
-                                <krel_min> 1e-9 </krel_min>
-                            </relative_permeability>
-                            <relative_permeability id="1">
-                                <type>WettingPhaseBrooksCoreyOilGas</type>
-                                <sr>  0.0 </sr>
-                                <smax> 1.0</smax>
-                                <m> 2 </m>
-                                <krel_min> 0.0 </krel_min>
-                            </relative_permeability>
-                        </relative_permeability>
-                    </porous_medium>
-                </porous_medium>
-            </material_property>
             <secondary_variables>
                 <secondary_variable type="static" internal_name="saturation" output_name="saturation"/>
             </secondary_variables>
@@ -84,6 +24,82 @@
             <temperature> temp </temperature>
         </process>
     </processes>
+
+ <media>
+      <medium id="0">
+         <phases>
+            <phase>
+               <type>AqueousLiquid</type>
+               <properties>
+                  <property>
+                     <name>density</name>
+                     <type>Constant</type>
+                     <value>1.0e3</value>
+                  </property>
+                  <property>
+                     <name>viscosity</name>
+                     <type>Constant</type>
+                     <value>1e-3</value>
+                  </property>
+               </properties>
+            </phase>
+            <phase>
+               <type>Gas</type>
+               <properties>
+                  <property>
+                     <name>density</name>
+                     <type>Constant</type>
+                     <value>1.0e3</value>
+                  </property>
+                  <property>
+                     <name>viscosity</name>
+                     <type>Constant</type>
+                     <value>1.e-3</value>
+                  </property>
+                  <property>
+                     <name>molar_mass</name>
+                     <type>Constant</type>
+                     <value>0.0289</value>
+                  </property>
+               </properties>
+            </phase>
+         </phases>
+         <properties>
+            <property>
+               <name>permeability</name>
+               <type>Constant</type>
+               <value>1.0e-10</value>
+            </property>
+
+            <property>
+               <name>saturation</name>
+               <type>SaturationBrooksCorey</type>
+               <residual_liquid_saturation>0.0</residual_liquid_saturation>
+               <residual_gas_saturation>0.0</residual_gas_saturation>
+               <lambda>2.0</lambda>
+               <entry_pressure>5000</entry_pressure>
+            </property>
+
+            <property>
+               <name>relative_permeability</name>
+               <type>RelPermBrooksCorey</type>
+               <residual_liquid_saturation>0.0</residual_liquid_saturation>
+               <residual_gas_saturation>0.0</residual_gas_saturation>
+               <lambda>2</lambda>
+               <min_relative_permeability_liquid>0</min_relative_permeability_liquid>
+               <min_relative_permeability_gas>1e-9</min_relative_permeability_gas>
+            </property>
+
+            <property>
+               <name>porosity</name>
+               <type>Constant</type>
+               <value>0.3</value>
+            </property>
+
+         </properties>
+      </medium>
+   </media>
+
     <parameters>
         <parameter>
             <name>pc_Dirichlet_left</name>
@@ -218,7 +234,7 @@
         <nonlinear_solver>
             <name>     basic_newton </name>
             <type>     Newton       </type>
-            <max_iter> 50          </max_iter>
+            <max_iter> 60          </max_iter>
             <linear_solver>general_linear_solver</linear_solver>
         </nonlinear_solver>
     </nonlinear_solvers>
diff --git a/Tests/MaterialLib/TestMPLSaturationBrooksCorey.cpp b/Tests/MaterialLib/TestMPLSaturationBrooksCorey.cpp
index 45700e5fdc8b9dccd15d6d9554a606a4121fa8b1..776403b9bcaa98046143828b2b6c661634897c2f 100644
--- a/Tests/MaterialLib/TestMPLSaturationBrooksCorey.cpp
+++ b/Tests/MaterialLib/TestMPLSaturationBrooksCorey.cpp
@@ -47,38 +47,23 @@ TEST(MaterialPropertyLib, SaturationBrooksCorey)
 
     MaterialPropertyLib::VariableArray variable_array;
     ParameterLib::SpatialPosition const pos;
-    double const time = std::numeric_limits<double>::quiet_NaN();
-
-    variable_array[static_cast<int>(
-        MaterialPropertyLib::Variable::capillary_pressure)] = 0.0;
-
-    auto s_L = medium->property(MaterialPropertyLib::PropertyType::saturation)
-                   .template value<double>(variable_array, pos, time);
-    auto ds_L_dp_cap =
-        medium->property(MaterialPropertyLib::PropertyType::saturation)
-            .template dValue<double>(
-                variable_array,
-                MaterialPropertyLib::Variable::capillary_pressure,
-                pos,
-                time);
-
-    ASSERT_EQ(s_L, max_saturation);
-    ASSERT_EQ(ds_L_dp_cap, 0.);
+    double const t = std::numeric_limits<double>::quiet_NaN();
 
     for (double p_cap = 1.0; p_cap < 1.0e10; p_cap *= 1.5)
     {
         variable_array[static_cast<int>(
             MaterialPropertyLib::Variable::capillary_pressure)] = p_cap;
 
-        s_L = medium->property(MaterialPropertyLib::PropertyType::saturation)
-                  .template value<double>(variable_array, pos, time);
-        ds_L_dp_cap =
+        auto s_L =
+            medium->property(MaterialPropertyLib::PropertyType::saturation)
+                .template value<double>(variable_array, pos, t);
+        auto ds_L_dp_cap =
             medium->property(MaterialPropertyLib::PropertyType::saturation)
                 .template dValue<double>(
                     variable_array,
                     MaterialPropertyLib::Variable::capillary_pressure,
                     pos,
-                    time);
+                    t);
 
         const double s_eff =
             std::pow(ref_entry_pressure / std::max(p_cap, ref_entry_pressure),
@@ -87,11 +72,12 @@ TEST(MaterialPropertyLib, SaturationBrooksCorey)
             s_eff * (max_saturation - ref_residual_liquid_saturation) +
             ref_residual_liquid_saturation;
         const double ds_eff_dpc =
-            (p_cap <= ref_entry_pressure) ? 0. : -ref_lambda / p_cap * s_ref;
-            const double ds_L_ds_eff = 1. / (max_saturation - ref_residual_liquid_saturation);
-            const double ds_L_dpc = ds_L_ds_eff * ds_eff_dpc;
+            -ref_lambda / std::max(p_cap, ref_entry_pressure) * s_ref;
+        const double ds_L_ds_eff =
+            1. / (max_saturation - ref_residual_liquid_saturation);
+        const double ds_L_dpc = ds_L_ds_eff * ds_eff_dpc;
 
         ASSERT_NEAR(s_L, s_ref, 1.e-10);
         ASSERT_NEAR(ds_L_dp_cap, ds_L_dpc, 1.e-10);
     }
-}
\ No newline at end of file
+}