diff --git a/Applications/ApplicationsLib/ProjectData.cpp b/Applications/ApplicationsLib/ProjectData.cpp
index b59d477fec4b78c6a9af3cbb29105fcf44057232..b48ceea6af7a6c9b9296156b171d61ef5e51fd59 100644
--- a/Applications/ApplicationsLib/ProjectData.cpp
+++ b/Applications/ApplicationsLib/ProjectData.cpp
@@ -1101,7 +1101,7 @@ void ProjectData::parseProcesses(
                 createThermalTwoPhaseFlowWithPPProcess(
                     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/CheckMaterialSpatialDistributionMap.h b/MaterialLib/MPL/CheckMaterialSpatialDistributionMap.h
index ae563862f6a3396362e5207dd9292110ffdd1898..a5d8edf00c54c1e1fcf6c2ae72efa038243291c6 100644
--- a/MaterialLib/MPL/CheckMaterialSpatialDistributionMap.h
+++ b/MaterialLib/MPL/CheckMaterialSpatialDistributionMap.h
@@ -18,13 +18,15 @@ namespace MaterialPropertyLib
 {
 template <typename ContainerMedium,
           typename ContainerSolid,
-          typename ContainerLiquid>
+          typename ContainerLiquid,
+          typename ContainerGas>
 void checkMaterialSpatialDistributionMap(
     MeshLib::Mesh const& mesh,
     MaterialPropertyLib::MaterialSpatialDistributionMap const& media_map,
     ContainerMedium const& required_properties_medium,
     ContainerSolid const& required_properties_solid_phase,
-    ContainerLiquid const& required_properties_liquid_phase)
+    ContainerLiquid const& required_properties_liquid_phase,
+    ContainerGas const& required_properties_gas_phase)
 {
     for (auto const& element : mesh.getElements())
     {
@@ -42,6 +44,11 @@ void checkMaterialSpatialDistributionMap(
                 medium.phase("AqueousLiquid"),
                 required_properties_liquid_phase);
         }
+        if (!required_properties_gas_phase.empty())
+        {
+            MaterialPropertyLib::checkRequiredProperties(
+                medium.phase("Gas"), required_properties_gas_phase);
+        }
         if (!required_properties_solid_phase.empty())
         {
             MaterialPropertyLib::checkRequiredProperties(
diff --git a/MaterialLib/TwoPhaseModels/CreateTwoPhaseFlowMaterialProperties.cpp b/MaterialLib/TwoPhaseModels/CreateTwoPhaseFlowMaterialProperties.cpp
index 857682891a01fd573e7be376f8c335a40a112c8d..1256b3e9ec3f03fe5f23fa6d830a6d6f779bd175 100644
--- a/MaterialLib/TwoPhaseModels/CreateTwoPhaseFlowMaterialProperties.cpp
+++ b/MaterialLib/TwoPhaseModels/CreateTwoPhaseFlowMaterialProperties.cpp
@@ -43,10 +43,6 @@ createTwoPhaseFlowMaterialProperties(
     //! \ogs_file_param{material__twophase_flow__material_property__fluid__liquid_density}
     auto const& rho_conf = fluid_config.getConfigSubtree("liquid_density");
     auto liquid_density = MaterialLib::Fluid::createFluidDensityModel(rho_conf);
-    //! \ogs_file_param{material__twophase_flow__material_property__fluid__gas_density}
-    auto const& rho_gas_conf = fluid_config.getConfigSubtree("gas_density");
-    auto gas_density =
-        MaterialLib::Fluid::createFluidDensityModel(rho_gas_conf);
     //! \ogs_file_param{material__twophase_flow__material_property__fluid__liquid_viscosity}
     auto const& mu_conf = fluid_config.getConfigSubtree("liquid_viscosity");
     auto liquid_viscosity = MaterialLib::Fluid::createViscosityModel(mu_conf);
@@ -130,8 +126,8 @@ createTwoPhaseFlowMaterialProperties(
     return std::forward_as_tuple(
         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(liquid_viscosity), 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)),
diff --git a/MaterialLib/TwoPhaseModels/TwoPhaseFlowWithPPMaterialProperties.cpp b/MaterialLib/TwoPhaseModels/TwoPhaseFlowWithPPMaterialProperties.cpp
index 80c0f9b79b5b9fa84811d7d4469529d7afdc1023..0b17cecd31af0161a73a4cf3cc93f31f2d6e0ab3 100644
--- a/MaterialLib/TwoPhaseModels/TwoPhaseFlowWithPPMaterialProperties.cpp
+++ b/MaterialLib/TwoPhaseModels/TwoPhaseFlowWithPPMaterialProperties.cpp
@@ -32,7 +32,6 @@ TwoPhaseFlowWithPPMaterialProperties::TwoPhaseFlowWithPPMaterialProperties(
     MeshLib::PropertyVector<int> 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,
@@ -49,7 +48,6 @@ TwoPhaseFlowWithPPMaterialProperties::TwoPhaseFlowWithPPMaterialProperties(
     : _material_ids(material_ids),
       _liquid_density(std::move(liquid_density)),
       _liquid_viscosity(std::move(liquid_viscosity)),
-      _gas_density(std::move(gas_density)),
       _gas_viscosity(std::move(gas_viscosity)),
       _intrinsic_permeability_models(std::move(intrinsic_permeability_models)),
       _porosity_models(std::move(porosity_models)),
@@ -81,15 +79,6 @@ double TwoPhaseFlowWithPPMaterialProperties::getLiquidDensity(
     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::getLiquidViscosity(
     const double p, const double T) const
 {
diff --git a/MaterialLib/TwoPhaseModels/TwoPhaseFlowWithPPMaterialProperties.h b/MaterialLib/TwoPhaseModels/TwoPhaseFlowWithPPMaterialProperties.h
index 5cce8068a52643f7e1eb2ff5f785d6aad31efce1..27f19500b124f0c62098e84508fa8c01af9a307b 100644
--- a/MaterialLib/TwoPhaseModels/TwoPhaseFlowWithPPMaterialProperties.h
+++ b/MaterialLib/TwoPhaseModels/TwoPhaseFlowWithPPMaterialProperties.h
@@ -56,8 +56,6 @@ public:
             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>>&&
diff --git a/ProcessLib/HT/CreateHTProcess.cpp b/ProcessLib/HT/CreateHTProcess.cpp
index 5eb5e983c1e3befa2fa8d7b80519226f55950976..ed8cc003b17bf9eb40e48ddaa704d358745c0bf5 100644
--- a/ProcessLib/HT/CreateHTProcess.cpp
+++ b/ProcessLib/HT/CreateHTProcess.cpp
@@ -49,9 +49,13 @@ void checkMPLProperties(
         MaterialPropertyLib::PropertyType::thermal_conductivity,
         MaterialPropertyLib::PropertyType::storage};
 
+    std::array<MaterialPropertyLib::PropertyType, 0> const
+        required_gas_properties{};
+
     MaterialPropertyLib::checkMaterialSpatialDistributionMap(
         mesh, media_map, required_property_medium,
-        required_property_solid_phase, required_property_liquid_phase);
+        required_property_solid_phase, required_property_liquid_phase,
+        required_gas_properties);
 }
 
 std::unique_ptr<Process> createHTProcess(
diff --git a/ProcessLib/HeatConduction/CreateHeatConductionProcess.cpp b/ProcessLib/HeatConduction/CreateHeatConductionProcess.cpp
index 9015dac356feb569a9d31128d65fdc977d8201fa..14c31c6f88d77d07850562834d6f087963b44d16 100644
--- a/ProcessLib/HeatConduction/CreateHeatConductionProcess.cpp
+++ b/ProcessLib/HeatConduction/CreateHeatConductionProcess.cpp
@@ -34,7 +34,7 @@ void checkMPLProperties(
     std::array<MaterialPropertyLib::PropertyType, 0> const empty{};
 
     MaterialPropertyLib::checkMaterialSpatialDistributionMap(
-        mesh, media_map, required_medium_properties, empty, empty);
+        mesh, media_map, required_medium_properties, empty, empty, empty);
 }
 
 std::unique_ptr<Process> createHeatConductionProcess(
diff --git a/ProcessLib/LiquidFlow/CreateLiquidFlowProcess.cpp b/ProcessLib/LiquidFlow/CreateLiquidFlowProcess.cpp
index f7e266ee46ac4596c473725b8e376c9a7ade67e4..61c646866cc97d3ae8876a128c8fb8dbede9af7a 100644
--- a/ProcessLib/LiquidFlow/CreateLiquidFlowProcess.cpp
+++ b/ProcessLib/LiquidFlow/CreateLiquidFlowProcess.cpp
@@ -44,9 +44,12 @@ void checkMPLProperties(
     std::array<MaterialPropertyLib::PropertyType, 0> const
         required_solid_properties{};
 
+    std::array<MaterialPropertyLib::PropertyType, 0> const
+        required_gas_properties{};
+
     MaterialPropertyLib::checkMaterialSpatialDistributionMap(
         mesh, media_map, required_medium_properties, required_solid_properties,
-        required_liquid_properties);
+        required_liquid_properties, required_gas_properties);
 }
 
 std::unique_ptr<Process> createLiquidFlowProcess(
diff --git a/ProcessLib/RichardsFlow/CreateRichardsFlowProcess.cpp b/ProcessLib/RichardsFlow/CreateRichardsFlowProcess.cpp
index bba81507a474bf7d71849ca08490a2c82321db4f..f32e7396950050e981e329df9d145703905486b1 100644
--- a/ProcessLib/RichardsFlow/CreateRichardsFlowProcess.cpp
+++ b/ProcessLib/RichardsFlow/CreateRichardsFlowProcess.cpp
@@ -42,9 +42,11 @@ void checkMPLProperties(
     std::array<MaterialPropertyLib::PropertyType, 0> const
         required_solid_properties{};
 
+    std::array<MaterialPropertyLib::PropertyType, 0> const empty{};
+
     MaterialPropertyLib::checkMaterialSpatialDistributionMap(
         mesh, media_map, required_medium_properties, required_solid_properties,
-        required_liquid_properties);
+        required_liquid_properties, empty);
 }
 
 std::unique_ptr<Process> createRichardsFlowProcess(
diff --git a/ProcessLib/SteadyStateDiffusion/CreateSteadyStateDiffusion.cpp b/ProcessLib/SteadyStateDiffusion/CreateSteadyStateDiffusion.cpp
index 2c44652c8dbc3c0577554f432d9d2006cb5aaa31..81179dce34ee92ca3d45f38f3a07f16128b74a65 100644
--- a/ProcessLib/SteadyStateDiffusion/CreateSteadyStateDiffusion.cpp
+++ b/ProcessLib/SteadyStateDiffusion/CreateSteadyStateDiffusion.cpp
@@ -35,7 +35,7 @@ void checkMPLProperties(
     std::array<MaterialPropertyLib::PropertyType, 0> const empty{};
 
     MaterialPropertyLib::checkMaterialSpatialDistributionMap(
-        mesh, media_map, required_medium_properties, empty, empty);
+        mesh, media_map, required_medium_properties, empty, empty, empty);
 }
 
 std::unique_ptr<Process> createSteadyStateDiffusion(
diff --git a/ProcessLib/ThermalTwoPhaseFlowWithPP/CreateThermalTwoPhaseFlowWithPPMaterialProperties.cpp b/ProcessLib/ThermalTwoPhaseFlowWithPP/CreateThermalTwoPhaseFlowWithPPMaterialProperties.cpp
index 7e73637fcf011042643e084edf2acf5f07e8ff76..5b931b5653ff41734b04d184aeca7e843c2942e8 100644
--- a/ProcessLib/ThermalTwoPhaseFlowWithPP/CreateThermalTwoPhaseFlowWithPPMaterialProperties.cpp
+++ b/ProcessLib/ThermalTwoPhaseFlowWithPP/CreateThermalTwoPhaseFlowWithPPMaterialProperties.cpp
@@ -49,19 +49,6 @@ createThermalTwoPhaseFlowWithPPMaterialProperties(
         std::move(std::get<0>(two_phase_model_tuple));
     auto const& fluid_config = std::get<1>(two_phase_model_tuple);
 
-    // Get fluid properties
-    auto const& spec_heat_capacity_solid_conf =
-        //! \ogs_file_param{prj__processes__process__TWOPHASE_FLOW_THERMAL__material_property__specific_heat_capacity_solid}
-        fluid_config.getConfigSubtree("specific_heat_capacity_solid");
-    auto specific_heat_capacity_solid =
-        MaterialLib::Fluid::createSpecificFluidHeatCapacityModel(
-            spec_heat_capacity_solid_conf);
-    auto const& spec_heat_capacity_water_conf =
-        //! \ogs_file_param{prj__processes__process__TWOPHASE_FLOW_THERMAL__material_property__specific_heat_capacity_water}
-        fluid_config.getConfigSubtree("specific_heat_capacity_water");
-    auto specific_heat_capacity_water =
-        MaterialLib::Fluid::createSpecificFluidHeatCapacityModel(
-            spec_heat_capacity_water_conf);
     auto const& spec_heat_capacity_air_conf =
         //! \ogs_file_param{prj__processes__process__TWOPHASE_FLOW_THERMAL__material_property__specific_heat_capacity_air}
         fluid_config.getConfigSubtree("specific_heat_capacity_air");
@@ -93,8 +80,6 @@ createThermalTwoPhaseFlowWithPPMaterialProperties(
 
     return std::make_unique<ThermalTwoPhaseFlowWithPPMaterialProperties>(
         std::move(two_phase_material_model),
-        std::move(specific_heat_capacity_solid),
-        std::move(specific_heat_capacity_water),
         std::move(specific_heat_capacity_air),
         std::move(specific_heat_capacity_vapor),
         std::move(thermal_conductivity_dry_solid),
diff --git a/ProcessLib/ThermalTwoPhaseFlowWithPP/CreateThermalTwoPhaseFlowWithPPProcess.cpp b/ProcessLib/ThermalTwoPhaseFlowWithPP/CreateThermalTwoPhaseFlowWithPPProcess.cpp
index b8fda16a42559cde27d79059bc7d2822426834fd..e5f0efdac728014fdb7e85bbc0b34a75419665cb 100644
--- a/ProcessLib/ThermalTwoPhaseFlowWithPP/CreateThermalTwoPhaseFlowWithPPProcess.cpp
+++ b/ProcessLib/ThermalTwoPhaseFlowWithPP/CreateThermalTwoPhaseFlowWithPPProcess.cpp
@@ -11,6 +11,8 @@
 
 #include <cassert>
 
+#include "MaterialLib/MPL/CheckMaterialSpatialDistributionMap.h"
+#include "MaterialLib/MPL/CreateMaterialSpatialDistributionMap.h"
 #include "ParameterLib/ConstantParameter.h"
 #include "ParameterLib/Utils.h"
 #include "ProcessLib/Output/CreateSecondaryVariables.h"
@@ -24,6 +26,49 @@ namespace ProcessLib
 {
 namespace ThermalTwoPhaseFlowWithPP
 {
+void checkMPLProperties(
+    std::map<int, std::shared_ptr<MaterialPropertyLib::Medium>> const& media)
+{
+    std::array const required_property_medium = {
+        MaterialPropertyLib::PropertyType::porosity,
+        MaterialPropertyLib::PropertyType::permeability};
+
+    std::array const required_property_solid_phase = {
+        MaterialPropertyLib::PropertyType::specific_heat_capacity,
+        MaterialPropertyLib::PropertyType::density};
+
+    std::array const required_property_liquid_phase = {
+        MaterialPropertyLib::PropertyType::viscosity,
+        MaterialPropertyLib::PropertyType::specific_heat_capacity,
+        MaterialPropertyLib::PropertyType::density};
+
+    std::array const required_property_gas_phase = {
+        MaterialPropertyLib::PropertyType::viscosity};
+
+    std::array const required_property_vapor_component = {
+        MaterialPropertyLib::specific_heat_capacity};
+
+    std::array const required_property_dry_air_component = {
+        MaterialPropertyLib::specific_heat_capacity};
+
+    for (auto const& m : media)
+    {
+        auto const& gas_phase = m.second->phase("Gas");
+        checkRequiredProperties(*m.second, required_property_medium);
+        checkRequiredProperties(gas_phase, required_property_gas_phase);
+        checkRequiredProperties(m.second->phase("AqueousLiquid"),
+                                required_property_liquid_phase);
+        checkRequiredProperties(m.second->phase("Solid"),
+                                required_property_solid_phase);
+
+        // TODO (BM): should use index to identify components (same for impl.h)
+        checkRequiredProperties(gas_phase.component("w"),
+                                required_property_vapor_component);
+        checkRequiredProperties(gas_phase.component("a"),
+                                required_property_dry_air_component);
+    }
+}
+
 std::unique_ptr<Process> createThermalTwoPhaseFlowWithPPProcess(
     std::string name,
     MeshLib::Mesh& mesh,
@@ -34,7 +79,8 @@ std::unique_ptr<Process> createThermalTwoPhaseFlowWithPPProcess(
     BaseLib::ConfigTree const& config,
     std::map<std::string,
              std::unique_ptr<MathLib::PiecewiseLinearInterpolation>> const&
-        curves)
+        curves,
+    std::map<int, std::shared_ptr<MaterialPropertyLib::Medium>> const& media)
 {
     //! \ogs_file_param{prj__processes__process__type}
     config.checkConfigParameter("type", "THERMAL_TWOPHASE_WITH_PP");
@@ -87,12 +133,6 @@ std::unique_ptr<Process> createThermalTwoPhaseFlowWithPPProcess(
 
     // Parameter for the density of the solid.
 
-    auto const& density_solid = ParameterLib::findParameter<double>(
-        config,
-        //! \ogs_file_param_special{prj__processes__process__TWOPHASE_FLOW_THERMAL__density_solid}
-        "density_solid", parameters, 1, &mesh);
-    DBUG("Use '{:s}' as density_solid parameter.", density_solid.name);
-
     // Parameter for the latent heat of evaporation.
     auto const& latent_heat_evaporation = ParameterLib::findParameter<double>(
         config,
@@ -108,12 +148,20 @@ std::unique_ptr<Process> createThermalTwoPhaseFlowWithPPProcess(
         createThermalTwoPhaseFlowWithPPMaterialProperties(
             mat_config, materialIDs(mesh), parameters);
 
-    ThermalTwoPhaseFlowWithPPProcessData process_data{specific_body_force,
+    auto media_map =
+        MaterialPropertyLib::createMaterialSpatialDistributionMap(media, mesh);
+
+    DBUG(
+        "Check the media properties of ThermalTwoPhaseFlowWithPP  process ...");
+    checkMPLProperties(media);
+    DBUG("Media properties verified.");
+
+    ThermalTwoPhaseFlowWithPPProcessData process_data{std::move(media_map),
+                                                      specific_body_force,
                                                       has_gravity,
                                                       mass_lumping,
                                                       diff_coeff_b,
                                                       diff_coeff_a,
-                                                      density_solid,
                                                       latent_heat_evaporation,
                                                       std::move(material)};
 
diff --git a/ProcessLib/ThermalTwoPhaseFlowWithPP/CreateThermalTwoPhaseFlowWithPPProcess.h b/ProcessLib/ThermalTwoPhaseFlowWithPP/CreateThermalTwoPhaseFlowWithPPProcess.h
index b5e835030e9fabeacc4838ccfa2e5cf3150a9080..d5adf1836c16ac6cec6693931330d147457e7d9d 100644
--- a/ProcessLib/ThermalTwoPhaseFlowWithPP/CreateThermalTwoPhaseFlowWithPPProcess.h
+++ b/ProcessLib/ThermalTwoPhaseFlowWithPP/CreateThermalTwoPhaseFlowWithPPProcess.h
@@ -12,6 +12,7 @@
 
 #include <memory>
 
+#include "MaterialLib/MPL/Medium.h"
 #include "ProcessLib/Process.h"
 
 namespace ProcessLib
@@ -28,6 +29,7 @@ std::unique_ptr<Process> createThermalTwoPhaseFlowWithPPProcess(
     BaseLib::ConfigTree const& config,
     std::map<std::string,
              std::unique_ptr<MathLib::PiecewiseLinearInterpolation>> const&
-        curves);
+        curves,
+    std::map<int, std::shared_ptr<MaterialPropertyLib::Medium>> const& media);
 }  // namespace ThermalTwoPhaseFlowWithPP
 }  // namespace ProcessLib
diff --git a/ProcessLib/ThermalTwoPhaseFlowWithPP/ThermalTwoPhaseFlowWithPPLocalAssembler-impl.h b/ProcessLib/ThermalTwoPhaseFlowWithPP/ThermalTwoPhaseFlowWithPPLocalAssembler-impl.h
index 8789b5a132e2cc32afbb0e92fabc8e5ef4039203..7c228943ceee8fae46e0fefd24af0756a30e2ffd 100644
--- a/ProcessLib/ThermalTwoPhaseFlowWithPP/ThermalTwoPhaseFlowWithPPLocalAssembler-impl.h
+++ b/ProcessLib/ThermalTwoPhaseFlowWithPP/ThermalTwoPhaseFlowWithPPLocalAssembler-impl.h
@@ -43,11 +43,13 @@
 */
 #pragma once
 
-#include "ThermalTwoPhaseFlowWithPPLocalAssembler.h"
-
+#include "MaterialLib/MPL/Medium.h"
+#include "MaterialLib/MPL/Property.h"
+#include "MaterialLib/MPL/Utils/FormEffectiveThermalConductivity.h"
+#include "MaterialLib/MPL/Utils/FormEigenTensor.h"
 #include "MathLib/InterpolationAlgorithms/PiecewiseLinearInterpolation.h"
 #include "NumLib/Function/Interpolation.h"
-
+#include "ThermalTwoPhaseFlowWithPPLocalAssembler.h"
 #include "ThermalTwoPhaseFlowWithPPProcessData.h"
 
 namespace ProcessLib
@@ -57,7 +59,7 @@ namespace ThermalTwoPhaseFlowWithPP
 template <typename ShapeFunction, typename IntegrationMethod, int GlobalDim>
 void ThermalTwoPhaseFlowWithPPLocalAssembler<
     ShapeFunction, IntegrationMethod,
-    GlobalDim>::assemble(double const t, double const /*dt*/,
+    GlobalDim>::assemble(double const t, double const dt,
                          std::vector<double> const& local_x,
                          std::vector<double> const& /*local_xdot*/,
                          std::vector<double>& local_M_data,
@@ -150,35 +152,48 @@ void ThermalTwoPhaseFlowWithPPLocalAssembler<
     auto const pc_nodal_values =
         Eigen::Map<const NodalVectorType>(&local_x[num_nodes], num_nodes);
 
-    const Eigen::MatrixXd& perm = two_phase_material_model.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));
-    }
-
+    MaterialPropertyLib::VariableArray vars;
     for (unsigned ip = 0; ip < n_integration_points; ip++)
     {
+        auto const& ip_data = _ip_data[ip];
+        auto const& N = ip_data.N;
+        auto const& dNdx = ip_data.dNdx;
+        auto const& w = ip_data.integration_weight;
+        auto const& mass_operator = ip_data.mass_operator;
+        auto const& diffusion_operator = ip_data.diffusion_operator;
         double pg_int_pt = 0.;
         double pc_int_pt = 0.;
         double T_int_pt = 0.0;
-        NumLib::shapeFunctionInterpolate(local_x, _ip_data[ip].N, pg_int_pt,
-                                         pc_int_pt, T_int_pt);
-
-        double const density_water =
-            two_phase_material_model.getLiquidDensity(pg_int_pt, T_int_pt);
+        NumLib::shapeFunctionInterpolate(local_x, N, pg_int_pt, pc_int_pt,
+                                         T_int_pt);
+        double const ideal_gas_constant_times_T_int_pt =
+            IdealGasConstant * T_int_pt;
+        vars[static_cast<int>(MaterialPropertyLib::Variable::temperature)] =
+            T_int_pt;
+        vars[static_cast<int>(
+            MaterialPropertyLib::Variable::capillary_pressure)] = pc_int_pt;
+        vars[static_cast<int>(MaterialPropertyLib::Variable::phase_pressure)] =
+            pg_int_pt;
+
+        auto const& medium =
+            *_process_data.media_map->getMedium(this->_element.getID());
+        auto const& liquid_phase = medium.phase("AqueousLiquid");
+        auto const& solid_phase = medium.phase("Solid");
+        auto const& gas_phase = medium.phase("Gas");
+
+        auto const& vapor_component = gas_phase.component("w");
+        auto const& dry_air_component = gas_phase.component("a");
+
+        auto const density_water =
+            liquid_phase.property(MaterialPropertyLib::PropertyType::density)
+                .template value<double>(vars, pos, t, dt);
 
         double const Sw = two_phase_material_model.getSaturation(
             material_id, t, pos, pg_int_pt, T_int_pt, pc_int_pt);
 
         _saturation[ip] = Sw;
+        vars[static_cast<int>(
+            MaterialPropertyLib::Variable::liquid_saturation)] = Sw;
 
         double dSwdpc =
             (pc_int_pt > two_phase_material_model.getCapillaryPressure(
@@ -200,10 +215,12 @@ void ThermalTwoPhaseFlowWithPPLocalAssembler<
         double const X_gas_nonwet =
             x_gas_nonwet /
             (x_gas_nonwet + x_vapor_nonwet * water_mol_mass / air_mol_mass);
-        double const mol_density_nonwet = pg_int_pt / IdealGasConstant / T_int_pt;
+        double const mol_density_nonwet =
+            pg_int_pt / ideal_gas_constant_times_T_int_pt;
         double const mol_density_water = density_water / water_mol_mass;
 
-        double const d_mol_density_nonwet_d_pg = 1 / IdealGasConstant / T_int_pt;
+        double const d_mol_density_nonwet_d_pg =
+            1 / ideal_gas_constant_times_T_int_pt;
         double const d_p_vapor_nonwet_d_T =
             _process_data.material->calculateDerivativedPgwdT(
                 pc_int_pt, T_int_pt, density_water);
@@ -211,39 +228,50 @@ void ThermalTwoPhaseFlowWithPPLocalAssembler<
             _process_data.material->calculateDerivativedPgwdPC(
                 pc_int_pt, T_int_pt, density_water);
         double const d_mol_density_nonwet_d_T =
-            -pg_int_pt / IdealGasConstant / T_int_pt / T_int_pt;
+            -pg_int_pt / ideal_gas_constant_times_T_int_pt / T_int_pt;
         double const d_x_gas_nonwet_d_pg =
             p_vapor_nonwet / pg_int_pt / pg_int_pt;
         double const d_x_gas_nonwet_d_pc = -d_p_vapor_nonwet_d_pc / pg_int_pt;
         double const d_x_gas_nonwet_d_T = -d_p_vapor_nonwet_d_T / pg_int_pt;
 
         double const density_nonwet_gas =
-            p_gas_nonwet * air_mol_mass / IdealGasConstant / T_int_pt;
+            p_gas_nonwet * air_mol_mass / ideal_gas_constant_times_T_int_pt;
         double const density_nonwet_vapor =
-            p_vapor_nonwet * water_mol_mass / IdealGasConstant / T_int_pt;
+            p_vapor_nonwet * water_mol_mass / ideal_gas_constant_times_T_int_pt;
         double const density_nonwet = density_nonwet_gas + density_nonwet_vapor;
         double const density_wet = density_water;
-        double const density_solid = _process_data.density_solid(t, pos)[0];
+        auto const density_solid =
+            solid_phase.property(MaterialPropertyLib::PropertyType::density)
+                .template value<double>(vars, pos, t, dt);
         // Derivative of nonwet phase density in terms of T
         double const d_density_nonwet_d_T =
-            _process_data.material->calculatedDensityNonwetdT (
-                p_gas_nonwet, p_vapor_nonwet, pc_int_pt, T_int_pt, density_water);
+            _process_data.material->calculatedDensityNonwetdT(
+                p_gas_nonwet, p_vapor_nonwet, pc_int_pt, T_int_pt,
+                density_water);
 
         _pressure_wetting[ip] = pg_int_pt - pc_int_pt;
         // heat capacity of nonwet phase
         double const heat_capacity_dry_gas =
-            _process_data.material->getSpecificHeatCapacityAir(pg_int_pt,
-                                                               T_int_pt);
+            dry_air_component
+                .property(
+                    MaterialPropertyLib::PropertyType::specific_heat_capacity)
+                .template value<double>(vars, pos, t, dt);
         const double heat_capacity_water_vapor =
-            _process_data.material->getSpecificHeatCapacityVapor(pg_int_pt,
-                                                                 T_int_pt);
-
-        double const heat_capacity_water =
-            _process_data.material->getSpecificHeatCapacityWater(pg_int_pt,
-                                                                 T_int_pt);
-        double const heat_capacity_solid =
-            _process_data.material->getSpecificHeatCapacitySolid(pg_int_pt,
-                                                                 T_int_pt);
+            vapor_component
+                .property(
+                    MaterialPropertyLib::PropertyType::specific_heat_capacity)
+                .template value<double>(vars, pos, t, dt);
+
+        auto const heat_capacity_water =
+            liquid_phase
+                .property(
+                    MaterialPropertyLib::PropertyType::specific_heat_capacity)
+                .template value<double>(vars, pos, t, dt);
+        auto const heat_capacity_solid =
+            solid_phase
+                .property(
+                    MaterialPropertyLib::PropertyType::specific_heat_capacity)
+                .template value<double>(vars, pos, t, dt);
         double const latent_heat_evaporation =
             _process_data.latent_heat_evaporation(t, pos)[0];
 
@@ -273,34 +301,35 @@ void ThermalTwoPhaseFlowWithPPLocalAssembler<
             d_enthalpy_gas_nonwet_d_T * X_gas_nonwet;
         // Assemble M matrix
         // nonwetting
-        double const porosity = two_phase_material_model.getPorosity(
-            material_id, t, pos, pg_int_pt, T_int_pt, 0);
+        auto const porosity =
+            medium.property(MaterialPropertyLib::PropertyType::porosity)
+                .template value<double>(vars, pos, t, dt);
 
-        Mgp.noalias() += porosity *
-                         ((1 - Sw) * (mol_density_nonwet * d_x_gas_nonwet_d_pg +
-                                      x_gas_nonwet * d_mol_density_nonwet_d_pg)) *
-                         _ip_data[ip].mass_operator;
+        Mgp.noalias() +=
+            porosity *
+            ((1 - Sw) * (mol_density_nonwet * d_x_gas_nonwet_d_pg +
+                         x_gas_nonwet * d_mol_density_nonwet_d_pg)) *
+            mass_operator;
         Mgpc.noalias() += porosity *
                           ((1 - Sw) * mol_density_nonwet * d_x_gas_nonwet_d_pc -
                            mol_density_nonwet * x_gas_nonwet * dSwdpc) *
-                          _ip_data[ip].mass_operator;
-        Mgt.noalias() += porosity *
-                         ((1 - Sw) * (mol_density_nonwet * d_x_gas_nonwet_d_T +
-                                      x_gas_nonwet * d_mol_density_nonwet_d_T)) *
-                         _ip_data[ip].mass_operator;
-
-        Mlpc.noalias() +=
+                          mass_operator;
+        Mgt.noalias() +=
             porosity *
-            ((1 - Sw) * d_p_vapor_nonwet_d_pc / IdealGasConstant / T_int_pt +
-             mol_density_nonwet * x_vapor_nonwet * (-dSwdpc) +
-             dSwdpc * mol_density_water) *
-            _ip_data[ip].mass_operator;
-        Mlt.noalias() +=
-            porosity *
-            ((1 - Sw) *
-             (d_p_vapor_nonwet_d_T / IdealGasConstant / T_int_pt -
-              p_vapor_nonwet / IdealGasConstant / T_int_pt / T_int_pt)) *
-            _ip_data[ip].mass_operator;
+            ((1 - Sw) * (mol_density_nonwet * d_x_gas_nonwet_d_T +
+                         x_gas_nonwet * d_mol_density_nonwet_d_T)) *
+            mass_operator;
+
+        Mlpc.noalias() += porosity *
+                          ((1 - Sw) * d_p_vapor_nonwet_d_pc /
+                               ideal_gas_constant_times_T_int_pt +
+                           mol_density_nonwet * x_vapor_nonwet * (-dSwdpc) +
+                           dSwdpc * mol_density_water) *
+                          mass_operator;
+        Mlt.noalias() += porosity *
+                         ((1 - Sw) / ideal_gas_constant_times_T_int_pt *
+                          (d_p_vapor_nonwet_d_T - p_vapor_nonwet / T_int_pt)) *
+                         mass_operator;
 
         Mep.noalias() +=
             porosity *
@@ -309,27 +338,30 @@ void ThermalTwoPhaseFlowWithPPLocalAssembler<
              mol_density_nonwet * (water_mol_mass - air_mol_mass) *
                  d_x_gas_nonwet_d_pg * enthalpy_nonwet -
              1) *
-            (1 - Sw) * _ip_data[ip].mass_operator;
+            (1 - Sw) * mass_operator;
         Mepc.noalias() +=
-            porosity * (density_wet * internal_energy_wet -
-                        density_nonwet * internal_energy_nonwet) *
-                dSwdpc * _ip_data[ip].mass_operator +
-            porosity * ((water_mol_mass - air_mol_mass) * enthalpy_nonwet /
-                        IdealGasConstant / T_int_pt) *
-                (1 - Sw) * d_p_vapor_nonwet_d_pc * _ip_data[ip].mass_operator;
+            porosity *
+                (density_wet * internal_energy_wet -
+                 density_nonwet * internal_energy_nonwet) *
+                dSwdpc * mass_operator +
+            porosity *
+                ((water_mol_mass - air_mol_mass) * enthalpy_nonwet /
+                 ideal_gas_constant_times_T_int_pt) *
+                (1 - Sw) * d_p_vapor_nonwet_d_pc * mass_operator;
         Met.noalias() +=
             ((1 - porosity) * density_solid * heat_capacity_solid +
              porosity * ((1 - Sw) * (d_density_nonwet_d_T * enthalpy_nonwet +
                                      density_nonwet * d_enthalpy_nonwet_d_T) +
                          Sw * density_wet * heat_capacity_water)) *
-            _ip_data[ip].mass_operator;
+            mass_operator;
 
         // nonwet
         double const k_rel_nonwet =
             two_phase_material_model.getNonwetRelativePermeability(
                 t, pos, _pressure_wetting[ip], T_int_pt, Sw);
-        double const mu_nonwet = two_phase_material_model.getGasViscosity(
-            _pressure_wetting[ip], T_int_pt);
+        auto const mu_nonwet =
+            gas_phase.property(MaterialPropertyLib::PropertyType::viscosity)
+                .template value<double>(vars, pos, t, dt);
         double const lambda_nonwet = k_rel_nonwet / mu_nonwet;
         double const diffusion_coeff_component_gas =
             _process_data.diffusion_coeff_component_b(t, pos)[0];
@@ -338,104 +370,138 @@ void ThermalTwoPhaseFlowWithPPLocalAssembler<
         double const k_rel_wet =
             two_phase_material_model.getWetRelativePermeability(
                 t, pos, pg_int_pt, T_int_pt, Sw);
-        double const mu_wet =
-            two_phase_material_model.getLiquidViscosity(pg_int_pt, T_int_pt);
+        auto const mu_wet =
+            liquid_phase.property(MaterialPropertyLib::PropertyType::viscosity)
+                .template value<double>(vars, pos, t, dt);
         double const lambda_wet = k_rel_wet / mu_wet;
 
+        auto const permeability =
+            MaterialPropertyLib::formEigenTensor<GlobalDim>(
+                medium.property(MaterialPropertyLib::PropertyType::permeability)
+                    .value(vars, pos, t, dt));
+
         GlobalDimVectorType const velocity_nonwet =
-            -lambda_nonwet * permeability *
-            (_ip_data[ip].dNdx * pg_nodal_values);
+            -lambda_nonwet * permeability * (dNdx * pg_nodal_values);
         GlobalDimVectorType const velocity_wet =
             -lambda_wet * permeability *
-            (_ip_data[ip].dNdx * (pg_nodal_values - pc_nodal_values));
-
-        laplace_operator.noalias() = _ip_data[ip].dNdx.transpose() *
-                                     permeability * _ip_data[ip].dNdx *
-                                     _ip_data[ip].integration_weight;
-
-        Ket.noalias() +=
-            _ip_data[ip].integration_weight * _ip_data[ip].N.transpose() *
-                (d_density_nonwet_d_T * enthalpy_nonwet +
-                 density_nonwet * d_enthalpy_nonwet_d_T) *
-                velocity_nonwet.transpose() * _ip_data[ip].dNdx +
-            _ip_data[ip].integration_weight * _ip_data[ip].N.transpose() *
-                heat_capacity_water * density_water * velocity_wet.transpose() *
-                _ip_data[ip].dNdx;
-
-        double const heat_conductivity_dry_solid =
-            _process_data.material->getThermalConductivityDrySolid(pg_int_pt,
-                                                                   T_int_pt);
-        double const heat_conductivity_wet_solid =
-            _process_data.material->getThermalConductivityWetSolid(pg_int_pt,
-                                                                   T_int_pt);
-        double const heat_conductivity_unsaturated =
-            _process_data.material->calculateUnsatHeatConductivity(
-                t, pos, Sw, heat_conductivity_dry_solid,
-                heat_conductivity_wet_solid);
+            (dNdx * (pg_nodal_values - pc_nodal_values));
+
+        laplace_operator.noalias() = dNdx.transpose() * permeability * dNdx * w;
+
+        Ket.noalias() += w * N.transpose() *
+                             (d_density_nonwet_d_T * enthalpy_nonwet +
+                              density_nonwet * d_enthalpy_nonwet_d_T) *
+                             velocity_nonwet.transpose() * dNdx +
+                         w * N.transpose() * heat_capacity_water *
+                             density_water * velocity_wet.transpose() * dNdx;
+
         // Laplace
-        Kgp.noalias() +=
-            (mol_density_nonwet * x_gas_nonwet * lambda_nonwet) * laplace_operator +
-            ((1 - Sw) * porosity * diffusion_coeff_component_gas *
-             mol_density_nonwet * d_x_gas_nonwet_d_pg) *
-                _ip_data[ip].diffusion_operator;
+        Kgp.noalias() += (mol_density_nonwet * x_gas_nonwet * lambda_nonwet) *
+                             laplace_operator +
+                         ((1 - Sw) * porosity * diffusion_coeff_component_gas *
+                          mol_density_nonwet * d_x_gas_nonwet_d_pg) *
+                             diffusion_operator;
         Kgpc.noalias() += ((1 - Sw) * porosity * diffusion_coeff_component_gas *
                            mol_density_nonwet * d_x_gas_nonwet_d_pc) *
-                          _ip_data[ip].diffusion_operator;
+                          diffusion_operator;
         Kgt.noalias() += ((1 - Sw) * porosity * diffusion_coeff_component_gas *
                           mol_density_nonwet * d_x_gas_nonwet_d_T) *
-                         _ip_data[ip].diffusion_operator;
+                         diffusion_operator;
 
         Klp.noalias() += (mol_density_nonwet * x_vapor_nonwet * lambda_nonwet) *
                              laplace_operator +
                          mol_density_water * lambda_wet * laplace_operator -
                          ((1 - Sw) * porosity * diffusion_coeff_component_gas *
                           mol_density_nonwet * d_x_gas_nonwet_d_pg) *
-                             _ip_data[ip].diffusion_operator;
+                             diffusion_operator;
         Klpc.noalias() += (-mol_density_water * lambda_wet * laplace_operator) -
                           ((1 - Sw) * porosity * diffusion_coeff_component_gas *
                            mol_density_nonwet * d_x_gas_nonwet_d_pc) *
-                              _ip_data[ip].diffusion_operator;
+                              diffusion_operator;
         Klt.noalias() += -((1 - Sw) * porosity * diffusion_coeff_component_gas *
                            mol_density_nonwet * d_x_gas_nonwet_d_T) *
-                         _ip_data[ip].diffusion_operator;
+                         diffusion_operator;
 
-        Kep.noalias() +=
-            (lambda_nonwet * density_nonwet * enthalpy_nonwet +
-             lambda_wet * density_wet * enthalpy_wet) *
-                laplace_operator +
-            (1 - Sw) * porosity * diffusion_coeff_component_gas *
-                mol_density_nonwet * (air_mol_mass * enthalpy_nonwet_gas -
-                                  water_mol_mass * enthalpy_nonwet_vapor) *
-                d_x_gas_nonwet_d_pg * _ip_data[ip].diffusion_operator;
+        Kep.noalias() += (lambda_nonwet * density_nonwet * enthalpy_nonwet +
+                          lambda_wet * density_wet * enthalpy_wet) *
+                             laplace_operator +
+                         (1 - Sw) * porosity * diffusion_coeff_component_gas *
+                             mol_density_nonwet *
+                             (air_mol_mass * enthalpy_nonwet_gas -
+                              water_mol_mass * enthalpy_nonwet_vapor) *
+                             d_x_gas_nonwet_d_pg * diffusion_operator;
         Kepc.noalias() +=
             -lambda_wet * enthalpy_wet * density_wet * laplace_operator +
             (1 - Sw) * porosity * diffusion_coeff_component_gas *
-                mol_density_nonwet * (air_mol_mass * enthalpy_nonwet_gas -
-                                  water_mol_mass * enthalpy_nonwet_vapor) *
-                d_x_gas_nonwet_d_pc * _ip_data[ip].diffusion_operator;
-        Ket.noalias() +=
-            _ip_data[ip].dNdx.transpose() * heat_conductivity_unsaturated *
-                _ip_data[ip].dNdx * _ip_data[ip].integration_weight +
-            (1 - Sw) * porosity * diffusion_coeff_component_gas *
-                mol_density_nonwet * (air_mol_mass * enthalpy_nonwet_gas -
-                                  water_mol_mass * enthalpy_nonwet_vapor) *
-                d_x_gas_nonwet_d_T * _ip_data[ip].diffusion_operator;
+                mol_density_nonwet *
+                (air_mol_mass * enthalpy_nonwet_gas -
+                 water_mol_mass * enthalpy_nonwet_vapor) *
+                d_x_gas_nonwet_d_pc * diffusion_operator;
+
+        if (medium.hasProperty(
+                MaterialPropertyLib::PropertyType::thermal_conductivity))
+        {
+            auto const lambda =
+                medium
+                    .property(
+                        MaterialPropertyLib::PropertyType::thermal_conductivity)
+                    .value(vars, pos, t, dt);
+
+            GlobalDimMatrixType const heat_conductivity_unsaturated =
+                MaterialPropertyLib::formEigenTensor<GlobalDim>(lambda);
+
+            Ket.noalias() +=
+                dNdx.transpose() * heat_conductivity_unsaturated * dNdx * w +
+                (1 - Sw) * porosity * diffusion_coeff_component_gas *
+                    mol_density_nonwet *
+                    (air_mol_mass * enthalpy_nonwet_gas -
+                     water_mol_mass * enthalpy_nonwet_vapor) *
+                    d_x_gas_nonwet_d_T * diffusion_operator;
+        }
+        else
+        {
+            auto const thermal_conductivity_solid =
+                solid_phase
+                    .property(
+                        MaterialPropertyLib::PropertyType::thermal_conductivity)
+                    .value(vars, pos, t, dt);
+
+            auto const thermal_conductivity_fluid =
+                liquid_phase
+                    .property(
+                        MaterialPropertyLib::PropertyType::thermal_conductivity)
+                    .template value<double>(vars, pos, t, dt) *
+                Sw;
+
+            GlobalDimMatrixType const heat_conductivity_unsaturated =
+                MaterialPropertyLib::formEffectiveThermalConductivity<
+                    GlobalDim>(thermal_conductivity_solid,
+                               thermal_conductivity_fluid, porosity);
+
+            Ket.noalias() +=
+                dNdx.transpose() * heat_conductivity_unsaturated * dNdx * w +
+                (1 - Sw) * porosity * diffusion_coeff_component_gas *
+                    mol_density_nonwet *
+                    (air_mol_mass * enthalpy_nonwet_gas -
+                     water_mol_mass * enthalpy_nonwet_vapor) *
+                    d_x_gas_nonwet_d_T * diffusion_operator;
+        }
 
         if (_process_data.has_gravity)
         {
             auto const& b = _process_data.specific_body_force;
-            NodalVectorType gravity_operator = _ip_data[ip].dNdx.transpose() *
-                                               permeability * b *
-                                               _ip_data[ip].integration_weight;
-            Bg.noalias() +=
-                (mol_density_nonwet * x_gas_nonwet * lambda_nonwet * density_nonwet) *
-                gravity_operator;
-            Bl.noalias() +=
-                (mol_density_water * lambda_wet * density_wet +
-                 mol_density_nonwet * x_vapor_nonwet * lambda_nonwet * density_nonwet) *
-                gravity_operator;
+            NodalVectorType gravity_operator =
+                dNdx.transpose() * permeability * b * w;
+            Bg.noalias() += (mol_density_nonwet * x_gas_nonwet * lambda_nonwet *
+                             density_nonwet) *
+                            gravity_operator;
+            Bl.noalias() += (mol_density_water * lambda_wet * density_wet +
+                             mol_density_nonwet * x_vapor_nonwet *
+                                 lambda_nonwet * density_nonwet) *
+                            gravity_operator;
             Be.noalias() +=
-                (lambda_nonwet * density_nonwet * density_nonwet * enthalpy_nonwet +
+                (lambda_nonwet * density_nonwet * density_nonwet *
+                     enthalpy_nonwet +
                  lambda_wet * density_wet * density_wet * enthalpy_wet) *
                 gravity_operator;
         }  // end of has gravity
diff --git a/ProcessLib/ThermalTwoPhaseFlowWithPP/ThermalTwoPhaseFlowWithPPMaterialProperties.cpp b/ProcessLib/ThermalTwoPhaseFlowWithPP/ThermalTwoPhaseFlowWithPPMaterialProperties.cpp
index 24712ae8146ff843e9cc7c3485f433ab8d924c7f..cf360536bfa8ca217dd77edb7ef71b8b54ba3fff 100644
--- a/ProcessLib/ThermalTwoPhaseFlowWithPP/ThermalTwoPhaseFlowWithPPMaterialProperties.cpp
+++ b/ProcessLib/ThermalTwoPhaseFlowWithPP/ThermalTwoPhaseFlowWithPPMaterialProperties.cpp
@@ -38,10 +38,6 @@ ThermalTwoPhaseFlowWithPPMaterialProperties::
         std::unique_ptr<MaterialLib::TwoPhaseFlowWithPP::
                             TwoPhaseFlowWithPPMaterialProperties>&&
             two_phase_material_model,
-        std::unique_ptr<MaterialLib::Fluid::FluidProperty>&&
-            specific_heat_capacity_solid,
-        std::unique_ptr<MaterialLib::Fluid::FluidProperty>&&
-            specific_heat_capacity_water,
         std::unique_ptr<MaterialLib::Fluid::FluidProperty>&&
             specific_heat_capacity_air,
         std::unique_ptr<MaterialLib::Fluid::FluidProperty>&&
@@ -53,8 +49,6 @@ ThermalTwoPhaseFlowWithPPMaterialProperties::
         std::unique_ptr<MaterialLib::Fluid::WaterVaporProperties>&&
             water_vapor_properties)
     : _two_phase_material_model(std::move(two_phase_material_model)),
-      _specific_heat_capacity_solid(std::move(specific_heat_capacity_solid)),
-      _specific_heat_capacity_water(std::move(specific_heat_capacity_water)),
       _specific_heat_capacity_air(std::move(specific_heat_capacity_air)),
       _specific_heat_capacity_vapor(std::move(specific_heat_capacity_vapor)),
       _thermal_conductivity_dry_solid(
@@ -66,26 +60,6 @@ ThermalTwoPhaseFlowWithPPMaterialProperties::
     DBUG("Create material properties for non-isothermal two-phase flow model.");
 }
 
-double
-ThermalTwoPhaseFlowWithPPMaterialProperties::getSpecificHeatCapacitySolid(
-    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 _specific_heat_capacity_solid->getValue(vars);
-}
-
-double
-ThermalTwoPhaseFlowWithPPMaterialProperties::getSpecificHeatCapacityWater(
-    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 _specific_heat_capacity_water->getValue(vars);
-}
-
 double ThermalTwoPhaseFlowWithPPMaterialProperties::getSpecificHeatCapacityAir(
     const double p, const double T) const
 {
diff --git a/ProcessLib/ThermalTwoPhaseFlowWithPP/ThermalTwoPhaseFlowWithPPMaterialProperties.h b/ProcessLib/ThermalTwoPhaseFlowWithPP/ThermalTwoPhaseFlowWithPPMaterialProperties.h
index e050951d152f2be561f142cadd1a9581ba700b40..93d40864a64a9e562231b5c635030aa917f99a9b 100644
--- a/ProcessLib/ThermalTwoPhaseFlowWithPP/ThermalTwoPhaseFlowWithPPMaterialProperties.h
+++ b/ProcessLib/ThermalTwoPhaseFlowWithPP/ThermalTwoPhaseFlowWithPPMaterialProperties.h
@@ -53,10 +53,6 @@ public:
         std::unique_ptr<MaterialLib::TwoPhaseFlowWithPP::
                             TwoPhaseFlowWithPPMaterialProperties>&&
             two_phase_material_model,
-        std::unique_ptr<MaterialLib::Fluid::FluidProperty>&&
-            specific_heat_capacity_solid,
-        std::unique_ptr<MaterialLib::Fluid::FluidProperty>&&
-            specific_heat_capacity_water,
         std::unique_ptr<MaterialLib::Fluid::FluidProperty>&&
             specific_heat_capacity_air,
         std::unique_ptr<MaterialLib::Fluid::FluidProperty>&&
@@ -68,7 +64,6 @@ public:
         std::unique_ptr<MaterialLib::Fluid::WaterVaporProperties>&&
             water_vapor_properties);
 
-    double getSpecificHeatCapacitySolid(const double p, const double T) const;
     double getSpecificHeatCapacityWater(const double p, const double T) const;
     double getSpecificHeatCapacityAir(const double p, const double T) const;
     double getSpecificHeatCapacityVapor(const double p, const double T) const;
diff --git a/ProcessLib/ThermalTwoPhaseFlowWithPP/ThermalTwoPhaseFlowWithPPProcessData.h b/ProcessLib/ThermalTwoPhaseFlowWithPP/ThermalTwoPhaseFlowWithPPProcessData.h
index a2227d529bb36d46301e74a49de758e712d86888..349fa5b603c80a27e60ef6be2986019eea0d9463 100644
--- a/ProcessLib/ThermalTwoPhaseFlowWithPP/ThermalTwoPhaseFlowWithPPProcessData.h
+++ b/ProcessLib/ThermalTwoPhaseFlowWithPP/ThermalTwoPhaseFlowWithPPProcessData.h
@@ -10,6 +10,7 @@
 
 #pragma once
 
+#include "MaterialLib/MPL/MaterialSpatialDistributionMap.h"
 #include "ThermalTwoPhaseFlowWithPPMaterialProperties.h"
 
 namespace ProcessLib
@@ -21,13 +22,14 @@ namespace ThermalTwoPhaseFlowWithPP
 {
 struct ThermalTwoPhaseFlowWithPPProcessData
 {
+    std::unique_ptr<MaterialPropertyLib::MaterialSpatialDistributionMap>
+        media_map;
     Eigen::VectorXd const specific_body_force;
 
     bool const has_gravity;
     bool const has_mass_lumping;
     ParameterLib::Parameter<double> const& diffusion_coeff_component_b;
     ParameterLib::Parameter<double> const& diffusion_coeff_component_a;
-    ParameterLib::Parameter<double> const& density_solid;
     ParameterLib::Parameter<double> const& latent_heat_evaporation;
     std::unique_ptr<ThermalTwoPhaseFlowWithPPMaterialProperties> material;
 };
diff --git a/Tests/Data/Parabolic/ThermalTwoPhaseFlowPP/HeatPipe/Twophase_HeatPipe_quad_curve_large.prj b/Tests/Data/Parabolic/ThermalTwoPhaseFlowPP/HeatPipe/Twophase_HeatPipe_quad_curve_large.prj
index 11e963d38089c2c1320a743649bb5a5ee20af108..9dccb1611af57eed3cbd8b649112f2845d910a6b 100644
--- a/Tests/Data/Parabolic/ThermalTwoPhaseFlowPP/HeatPipe/Twophase_HeatPipe_quad_curve_large.prj
+++ b/Tests/Data/Parabolic/ThermalTwoPhaseFlowPP/HeatPipe/Twophase_HeatPipe_quad_curve_large.prj
@@ -23,10 +23,6 @@
                         <type>Constant</type>
                         <value> 1.e3 </value>
                     </liquid_density>
-                    <gas_density>
-                        <type>IdealGasLaw</type>
-                        <molar_mass> 0.029 </molar_mass>
-                    </gas_density>
                     <liquid_viscosity>
                         <type>Constant</type>
                         <value> 2.938e-4 </value>
@@ -35,14 +31,6 @@
                         <type>Constant</type>
                         <value> 1.8e-5 </value>
                     </gas_viscosity>
-                    <specific_heat_capacity_solid>
-                        <type>Constant</type>
-                        <value> 700 </value>
-                    </specific_heat_capacity_solid>
-                    <specific_heat_capacity_water>
-                        <type>Constant</type>
-                        <value> 4187 </value>
-                    </specific_heat_capacity_water>
                     <specific_heat_capacity_air>
                         <type>Constant</type>
                         <value> 733 </value>
@@ -110,10 +98,104 @@
             <mass_lumping> true </mass_lumping>
             <diffusion_coeff_component_b> diff_coef_b </diffusion_coeff_component_b>
             <diffusion_coeff_component_a> diff_coef_a </diffusion_coeff_component_a>
-            <density_solid> rho_solid </density_solid>
             <latent_heat_evaporation> latent_heat_evp </latent_heat_evaporation>
         </process>
     </processes>
+    <media>
+        <medium id="0">
+            <phases>
+                <phase>
+                    <type>Gas</type>
+                    <components>
+                        <component>
+                            <name>w</name>
+                            <properties>
+                                <property>
+                                    <name>specific_heat_capacity</name>
+                                    <type>Constant</type>
+                                    <value>4187</value>
+                                </property>
+                            </properties>
+                        </component>
+                        <component>
+                            <name>a</name>
+                            <properties>
+                                <property>
+                                    <name>specific_heat_capacity</name>
+                                    <type>Constant</type>
+                                    <value>733</value>
+                                </property>
+                            </properties>
+                        </component>
+                    </components>
+                    <properties>
+                        <property>
+                            <name>viscosity</name>
+                            <type>Constant</type>
+                            <value> 1.8e-5 </value>
+                        </property>
+                        <property>
+                            <name>density</name>
+                            <type>IdealGasLaw</type>
+                        </property>
+                    </properties>
+                </phase>
+                <phase>
+                    <type>AqueousLiquid</type>
+                    <properties>
+                        <property>
+                            <name>viscosity</name>
+                            <type>Constant</type>
+                            <value> 2.938e-4 </value>
+                        </property>
+                        <property>
+                            <name>specific_heat_capacity</name>
+                            <type>Constant</type>
+                            <value> 4187 </value>
+                        </property>
+                        <property>
+                            <name>density</name>
+                            <type>Constant</type>
+                            <value> 1.e3 </value>
+                        </property>
+                    </properties>
+                </phase>
+                <phase>
+                    <type>Solid</type>
+                    <properties>
+                        <property>
+                            <name>specific_heat_capacity</name>
+                            <type>Constant</type>
+                            <value> 700 </value>
+                        </property>
+                        <property>
+                            <name>density</name>
+                            <type>Constant</type>
+                            <value>2650</value>
+                        </property>
+                    </properties>
+                </phase>
+            </phases>
+            <properties>
+                <property>
+                    <name>porosity</name>
+                    <type>Constant</type>
+                    <value>0.4</value>
+                </property>
+                <property>
+                    <name>permeability</name>
+                    <type>Constant</type>
+                    <value>1e-12</value>
+                </property>
+                <property>
+                    <name>thermal_conductivity</name>
+                    <type>SoilThermalConductivitySomerton</type>
+                    <dry_thermal_conductivity>k_T_dry</dry_thermal_conductivity>
+                    <wet_thermal_conductivity>k_T_wet</wet_thermal_conductivity>
+                </property>
+            </properties>
+        </medium>
+    </media>
     <parameters>
         <parameter>
             <name>HEN0</name>
@@ -170,11 +252,6 @@
             <type>Constant</type>
             <value>0.0</value>
         </parameter>
-        <parameter>
-            <name>rho_solid</name>
-            <type>Constant</type>
-            <value>2650</value>
-        </parameter>
         <parameter>
             <name>latent_heat_evp</name>
             <type>Constant</type>
@@ -190,6 +267,16 @@
             <type>Constant</type>
             <values>1e-12</values>
         </parameter>
+        <parameter>
+            <name>k_T_dry</name>
+            <type>Constant</type>
+            <value>0.582</value>
+        </parameter>
+        <parameter>
+            <name>k_T_wet</name>
+            <type>Constant</type>
+            <value>1.14</value>
+        </parameter>
     </parameters>
     <time_loop>
         <processes>
diff --git a/Tests/Data/Parabolic/ThermalTwoPhaseFlowPP/HeatPipe/Twophase_HeatPipe_quad_curve_small.prj b/Tests/Data/Parabolic/ThermalTwoPhaseFlowPP/HeatPipe/Twophase_HeatPipe_quad_curve_small.prj
index 8461eb5339684677f6f12c958346bd62ce5c5f02..1b91e86e8e020e5c3fc63de872e65e25efab62a1 100644
--- a/Tests/Data/Parabolic/ThermalTwoPhaseFlowPP/HeatPipe/Twophase_HeatPipe_quad_curve_small.prj
+++ b/Tests/Data/Parabolic/ThermalTwoPhaseFlowPP/HeatPipe/Twophase_HeatPipe_quad_curve_small.prj
@@ -23,10 +23,6 @@
                         <type>Constant</type>
                         <value> 1.e3 </value>
                     </liquid_density>
-                    <gas_density>
-                        <type>IdealGasLaw</type>
-                        <molar_mass> 0.029 </molar_mass>
-                    </gas_density>
                     <liquid_viscosity>
                         <type>Constant</type>
                         <value> 2.938e-4 </value>
@@ -35,14 +31,6 @@
                         <type>Constant</type>
                         <value> 1.8e-5 </value>
                     </gas_viscosity>
-                    <specific_heat_capacity_solid>
-                        <type>Constant</type>
-                        <value> 700 </value>
-                    </specific_heat_capacity_solid>
-                    <specific_heat_capacity_water>
-                        <type>Constant</type>
-                        <value> 4187 </value>
-                    </specific_heat_capacity_water>
                     <specific_heat_capacity_air>
                         <type>Constant</type>
                         <value> 733 </value>
@@ -110,10 +98,104 @@
             <mass_lumping> true </mass_lumping>
             <diffusion_coeff_component_b> diff_coef_b </diffusion_coeff_component_b>
             <diffusion_coeff_component_a> diff_coef_a </diffusion_coeff_component_a>
-            <density_solid> rho_solid </density_solid>
             <latent_heat_evaporation> latent_heat_evp </latent_heat_evaporation>
         </process>
     </processes>
+    <media>
+        <medium id="0">
+            <phases>
+                <phase>
+                    <type>Gas</type>
+                    <components>
+                        <component>
+                            <name>w</name>
+                            <properties>
+                                <property>
+                                    <name>specific_heat_capacity</name>
+                                    <type>Constant</type>
+                                    <value>4187</value>
+                                </property>
+                            </properties>
+                        </component>
+                        <component>
+                            <name>a</name>
+                            <properties>
+                                <property>
+                                    <name>specific_heat_capacity</name>
+                                    <type>Constant</type>
+                                    <value>733</value>
+                                </property>
+                            </properties>
+                        </component>
+                    </components>
+                    <properties>
+                        <property>
+                            <name>viscosity</name>
+                            <type>Constant</type>
+                            <value> 1.8e-5 </value>
+                        </property>
+                        <property>
+                            <name>density</name>
+                            <type>IdealGasLaw</type>
+                        </property>
+                    </properties>
+                </phase>
+                <phase>
+                    <type>AqueousLiquid</type>
+                    <properties>
+                        <property>
+                            <name>viscosity</name>
+                            <type>Constant</type>
+                            <value> 2.938e-4 </value>
+                        </property>
+                        <property>
+                            <name>specific_heat_capacity</name>
+                            <type>Constant</type>
+                            <value> 4187 </value>
+                        </property>
+                        <property>
+                            <name>density</name>
+                            <type>Constant</type>
+                            <value> 1.e3 </value>
+                        </property>
+                    </properties>
+                </phase>
+                <phase>
+                    <type>Solid</type>
+                    <properties>
+                        <property>
+                            <name>specific_heat_capacity</name>
+                            <type>Constant</type>
+                            <value> 700 </value>
+                        </property>
+                        <property>
+                            <name>density</name>
+                            <type>Constant</type>
+                            <value>2650</value>
+                        </property>
+                    </properties>
+                </phase>
+            </phases>
+            <properties>
+                <property>
+                    <name>porosity</name>
+                    <type>Constant</type>
+                    <value>0.4</value>
+                </property>
+                <property>
+                    <name>permeability</name>
+                    <type>Constant</type>
+                    <value>1e-12</value>
+                </property>
+                <property>
+                    <name>thermal_conductivity</name>
+                    <type>SoilThermalConductivitySomerton</type>
+                    <dry_thermal_conductivity>k_T_dry</dry_thermal_conductivity>
+                    <wet_thermal_conductivity>k_T_wet</wet_thermal_conductivity>
+                </property>
+            </properties>
+        </medium>
+    </media>
     <parameters>
         <parameter>
             <name>HEN0</name>
@@ -170,11 +252,6 @@
             <type>Constant</type>
             <value>0.0</value>
         </parameter>
-        <parameter>
-            <name>rho_solid</name>
-            <type>Constant</type>
-            <value>2650</value>
-        </parameter>
         <parameter>
             <name>latent_heat_evp</name>
             <type>Constant</type>
@@ -190,6 +267,16 @@
             <type>Constant</type>
             <values>1e-12</values>
         </parameter>
+        <parameter>
+            <name>k_T_dry</name>
+            <type>Constant</type>
+            <value>0.582</value>
+        </parameter>
+        <parameter>
+            <name>k_T_wet</name>
+            <type>Constant</type>
+            <value>1.14</value>
+        </parameter>
     </parameters>
     <time_loop>
         <processes>