diff --git a/Documentation/ProjectFile/prj/processes/process/LIQUID_FLOW/material_property/solid/i_solid.md b/Documentation/ProjectFile/prj/processes/process/LIQUID_FLOW/material_property/solid/i_solid.md
deleted file mode 100644
index 10e806ec9120b9a321904b4b461bd8da1468dd0c..0000000000000000000000000000000000000000
--- a/Documentation/ProjectFile/prj/processes/process/LIQUID_FLOW/material_property/solid/i_solid.md
+++ /dev/null
@@ -1 +0,0 @@
-Property of solid phase of liquid flow process.
diff --git a/Documentation/ProjectFile/prj/processes/process/LIQUID_FLOW/material_property/solid/t_biot_constant.md b/Documentation/ProjectFile/prj/processes/process/LIQUID_FLOW/material_property/solid/t_biot_constant.md
deleted file mode 100644
index 38343bfd1f58f621e53be93674819bad19886aab..0000000000000000000000000000000000000000
--- a/Documentation/ProjectFile/prj/processes/process/LIQUID_FLOW/material_property/solid/t_biot_constant.md
+++ /dev/null
@@ -1 +0,0 @@
-Biot's constant for the liquid thermal expansion computation.
diff --git a/Documentation/ProjectFile/prj/processes/process/LIQUID_FLOW/material_property/solid/t_thermal_expansion.md b/Documentation/ProjectFile/prj/processes/process/LIQUID_FLOW/material_property/solid/t_thermal_expansion.md
deleted file mode 100644
index 984c939c457d5a62efcddc947cc69541f9354d42..0000000000000000000000000000000000000000
--- a/Documentation/ProjectFile/prj/processes/process/LIQUID_FLOW/material_property/solid/t_thermal_expansion.md
+++ /dev/null
@@ -1 +0,0 @@
-Liquid thermal expansion.
diff --git a/ProcessLib/HeatConduction/HeatConductionFEM-impl.h b/ProcessLib/HeatConduction/HeatConductionFEM-impl.h
deleted file mode 100644
index 8c21f2d0364c9b1b0967a2402ca24e59bb617c81..0000000000000000000000000000000000000000
--- a/ProcessLib/HeatConduction/HeatConductionFEM-impl.h
+++ /dev/null
@@ -1,167 +0,0 @@
-/**
- * \copyright
- * Copyright (c) 2012-2017, OpenGeoSys Community (http://www.opengeosys.org)
- *            Distributed under a Modified BSD License.
- *              See accompanying file LICENSE.txt or
- *              http://www.opengeosys.org/project/license
- *
- *  \file   HeatConductionFEM-impl.h
- *  Created on January 17, 2017, 3:41 PM
- */
-
-#pragma once
-
-#include "HeatConductionFEM.h"
-#include "NumLib/Function/Interpolation.h"
-
-namespace ProcessLib
-{
-namespace HeatConduction
-{
-template <typename ShapeFunction, typename IntegrationMethod,
-          unsigned GlobalDim>
-template <typename LiquidFlowVelocityCalculator>
-void LocalAssemblerData<ShapeFunction, IntegrationMethod, GlobalDim>::
-    assembleHeatTransportLiquidFlow(
-        double const t, int const material_id, SpatialPosition& pos,
-        int const gravitational_axis_id,
-        double const gravitational_acceleration,
-        Eigen::MatrixXd const& permeability,
-        ProcessLib::LiquidFlow::LiquidFlowMaterialProperties const&
-            liquid_flow_properties,
-        std::vector<double> const& local_x, std::vector<double> const& local_p,
-        std::vector<double>& local_M_data, std::vector<double>& local_K_data)
-{
-    auto const local_matrix_size = local_x.size();
-    // This assertion is valid only if all nodal d.o.f. use the same shape
-    // matrices.
-    assert(local_matrix_size == ShapeFunction::NPOINTS * NUM_NODAL_DOF);
-
-    auto local_M = MathLib::createZeroedMatrix<NodalMatrixType>(
-        local_M_data, local_matrix_size, local_matrix_size);
-    auto local_K = MathLib::createZeroedMatrix<NodalMatrixType>(
-        local_K_data, local_matrix_size, local_matrix_size);
-    const auto local_p_vec =
-        MathLib::toVector<NodalVectorType>(local_p, local_matrix_size);
-
-    unsigned const n_integration_points =
-        _integration_method.getNumberOfPoints();
-
-    const double porosity_variable = 0.;
-    for (unsigned ip = 0; ip < n_integration_points; ip++)
-    {
-        pos.setIntegrationPoint(ip);
-        auto const& sm = _shape_matrices[ip];
-        auto const& wp = _integration_method.getWeightedPoint(ip);
-        double p = 0.;
-        NumLib::shapeFunctionInterpolate(local_p, sm.N, p);
-        double T = 0.;
-        NumLib::shapeFunctionInterpolate(local_x, sm.N, T);
-
-        // Thermal conductivity of solid phase.
-        auto const k_s = _process_data.thermal_conductivity(t, pos)[0];
-        // Specific heat capacity of solid phase.
-        auto const cp_s = _process_data.heat_capacity(t, pos)[0];
-        // Density of solid phase.
-        auto const rho_s = _process_data.density(t, pos)[0];
-
-        // Thermal conductivity of liquid.
-        double const k_f = liquid_flow_properties.getThermalConductivity(p, T);
-        // Specific heat capacity of liquid.
-        double const cp_f = liquid_flow_properties.getHeatCapacity(p, T);
-        // Density of liquid.
-        double const rho_f = liquid_flow_properties.getLiquidDensity(p, T);
-
-        // Porosity of porous media.
-        double const n = liquid_flow_properties.getPorosity(
-            material_id, t, pos, porosity_variable, T);
-
-        // Effective specific heat capacity.
-        double const effective_cp = (1.0 - n) * cp_s * rho_s + n * cp_f * rho_f;
-        // Effective thermal conductivity.
-        double const effective_K = (1.0 - n) * k_s + n * k_f;
-
-        double const integration_factor =
-            sm.detJ * wp.getWeight() * sm.integralMeasure;
-
-        local_M.noalias() +=
-            effective_cp * sm.N.transpose() * sm.N * integration_factor;
-        local_K.noalias() +=
-            sm.dNdx.transpose() * effective_K * sm.dNdx * integration_factor;
-
-        // Compute fluid flow velocity:
-        double const mu =
-            liquid_flow_properties.getViscosity(p, T);  // Viscosity
-        GlobalDimVectorType const& darcy_velocity =
-            LiquidFlowVelocityCalculator::calculateVelocity(
-                local_p_vec, sm, permeability, mu,
-                rho_f * gravitational_acceleration, gravitational_axis_id);
-        // Add the advection term
-        local_K.noalias() += cp_f * rho_f * sm.N.transpose() *
-                             (darcy_velocity.transpose() * sm.dNdx) *
-                             integration_factor;
-    }
-}
-
-template <typename ShapeFunction, typename IntegrationMethod,
-          unsigned GlobalDim>
-void LocalAssemblerData<ShapeFunction, IntegrationMethod, GlobalDim>::
-    assembleWithCoupledTerm(double const t, std::vector<double> const& local_x,
-                            std::vector<double>& local_M_data,
-                            std::vector<double>& local_K_data,
-                            std::vector<double>& /*local_b_data*/,
-                            LocalCoupledSolutions const& coupled_term)
-{
-    /*
-    for (auto const& coupled_process_pair : coupled_term.coupled_processes)
-    {
-        if (coupled_process_pair.first ==
-            std::type_index(typeid(ProcessLib::LiquidFlow::LiquidFlowProcess)))
-        {
-            assert(
-                dynamic_cast<const ProcessLib::LiquidFlow::LiquidFlowProcess*>(
-                    &(coupled_process_pair.second)) != nullptr);
-
-            auto const& pcs =
-                static_cast<ProcessLib::LiquidFlow::LiquidFlowProcess const&>(
-                    coupled_process_pair.second);
-            const auto liquid_flow_prop = pcs.getLiquidFlowMaterialProperties();
-
-            const auto local_p =
-                coupled_term.local_coupled_xs.at(coupled_process_pair.first);
-
-            SpatialPosition pos;
-            pos.setElementID(_element.getID());
-            const int material_id = liquid_flow_prop->getMaterialID(pos);
-
-            const Eigen::MatrixXd& K = liquid_flow_prop->getPermeability(
-                material_id, t, pos, _element.getDimension());
-
-            if (K.size() == 1)
-            {
-                assembleHeatTransportLiquidFlow<
-                    IsotropicLiquidFlowVelocityCalculator>(
-                    t, material_id, pos, pcs.getGravitationalAxisID(),
-                    pcs.getGravitationalAcceleration(), K, *liquid_flow_prop,
-                    local_x, local_p, local_M_data, local_K_data);
-            }
-            else
-            {
-                assembleHeatTransportLiquidFlow<
-                    AnisotropicLiquidFlowVelocityCalculator>(
-                    t, material_id, pos, pcs.getGravitationalAxisID(),
-                    pcs.getGravitationalAcceleration(), K, *liquid_flow_prop,
-                    local_x, local_p, local_M_data, local_K_data);
-            }
-        }
-        else
-        {
-            OGS_FATAL(
-                "This coupled process is not presented for "
-                "HeatConduction process");
-        }
-    }*/
-}
-
-}  // namespace HeatConduction
-}  // namespace ProcessLib
diff --git a/ProcessLib/LiquidFlow/CreateLiquidFlowMaterialProperties.cpp b/ProcessLib/LiquidFlow/CreateLiquidFlowMaterialProperties.cpp
index 5ac06dc06d76cd11cd6a2d5cfdc12bea665a6f47..5ece9fa1970086d1366d82b73741ba4c1998800a 100644
--- a/ProcessLib/LiquidFlow/CreateLiquidFlowMaterialProperties.cpp
+++ b/ProcessLib/LiquidFlow/CreateLiquidFlowMaterialProperties.cpp
@@ -18,6 +18,7 @@
 #include "MeshLib/PropertyVector.h"
 
 #include "MaterialLib/Fluid/FluidProperty.h"
+#include "MaterialLib/PorousMedium/Permeability/Permeability.h"
 #include "MaterialLib/PorousMedium/Porosity/Porosity.h"
 #include "MaterialLib/PorousMedium/Storage/Storage.h"
 #include "MaterialLib/Fluid/FluidProperties/CreateFluidProperties.h"
@@ -26,7 +27,6 @@
 #include "MaterialLib/PorousMedium/PorousPropertyHeaders.h"
 
 #include "ProcessLib/Utils/ProcessUtils.h"
-#include "ProcessLib/Parameter/ConstantParameter.h"
 
 #include "LiquidFlowMaterialProperties.h"
 
@@ -97,32 +97,10 @@ createLiquidFlowMaterialProperties(
     BaseLib::reorderVector(porosity_models, mat_ids);
     BaseLib::reorderVector(storage_models, mat_ids);
 
-    //! \ogs_file_param{prj__processes__process__LIQUID_FLOW__material_property__solid}
-    auto const solid_config = config.getConfigSubtreeOptional("solid");
-    if (solid_config)
-    {
-        auto& solid_thermal_expansion = findParameter<double>(
-            //! \ogs_file_param_special{prj__processes__process__LIQUID_FLOW__material_property__solid__thermal_expansion}
-            *solid_config, "thermal_expansion", parameters, 1);
-        DBUG("Use \'%s\' as solid thermal expansion.",
-             solid_thermal_expansion.name.c_str());
-        auto& biot_constant = findParameter<double>(
-            //! \ogs_file_param_special{prj__processes__process__LIQUID_FLOW__material_property__solid__biot_constant}
-            *solid_config, "biot_constant", parameters, 1);
-        return std::make_unique<LiquidFlowMaterialProperties>(
-            std::move(fluid_properties),
-            std::move(intrinsic_permeability_models),
-            std::move(porosity_models), std::move(storage_models),
-            has_material_ids, material_ids, solid_thermal_expansion,
-            biot_constant);
-    }
-
-    ConstantParameter<double> void_parameter("void_solid_thermal_expansion",
-                                             0.);
     return std::make_unique<LiquidFlowMaterialProperties>(
         std::move(fluid_properties), std::move(intrinsic_permeability_models),
         std::move(porosity_models), std::move(storage_models), has_material_ids,
-        material_ids, void_parameter, void_parameter);
+        material_ids);
 }
 
 }  // end of namespace
diff --git a/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h
index 2d6262194aa4981e201c11d50b9a369aaca2f5a6..4a52debc71c0c7c0f2b56c1ba1c9207cea3e9f30 100644
--- a/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h
+++ b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h
@@ -16,10 +16,6 @@
 
 #include "NumLib/Function/Interpolation.h"
 
-#include "ProcessLib/CoupledSolutionsForStaggeredScheme.h"
-
-#include "ProcessLib/HeatConduction/HeatConductionProcess.h"
-
 namespace ProcessLib
 {
 namespace LiquidFlow
@@ -53,60 +49,6 @@ void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
             pos, permeability);
 }
 
-/*
-template <typename ShapeFunction, typename IntegrationMethod,
-          unsigned GlobalDim>
-void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
-    assembleWithCoupledTerm(double const t, std::vector<double> const& local_x,
-                            std::vector<double>& local_M_data,
-                            std::vector<double>& local_K_data,
-                            std::vector<double>& local_b_data,
-                            LocalCouplingTerm const& coupled_term)
-{
-    SpatialPosition pos;
-    pos.setElementID(_element.getID());
-    const int material_id = _material_properties.getMaterialID(pos);
-
-    const Eigen::MatrixXd& permeability = _material_properties.getPermeability(
-        material_id, t, pos, _element.getDimension());
-    // Note: For Inclined 1D in 2D/3D or 2D element in 3D, the first item in
-    //  the assert must be changed to permeability.rows() ==
-    //  _element->getDimension()
-    assert(permeability.rows() == GlobalDim || permeability.rows() == 1);
-
-    const double dt = coupled_term.dt;
-    for (auto const& coupled_process_pair : coupled_term.coupled_processes)
-    {
-        if (coupled_process_pair.first ==
-            std::type_index(
-                typeid(ProcessLib::HeatConduction::HeatConductionProcess)))
-        {
-            const auto local_T0 =
-                coupled_term.local_coupled_xs0.at(coupled_process_pair.first);
-            const auto local_T1 =
-                coupled_term.local_coupled_xs.at(coupled_process_pair.first);
-
-            if (permeability.size() == 1)  // isotropic or 1D problem.
-                assembleWithCoupledWithHeatTransport<IsotropicCalculator>(
-                    material_id, t, dt, local_x, local_T0, local_T1,
-                    local_M_data, local_K_data, local_b_data, pos,
-                    permeability);
-            else
-                assembleWithCoupledWithHeatTransport<AnisotropicCalculator>(
-                    material_id, t, dt, local_x, local_T0, local_T1,
-                    local_M_data, local_K_data, local_b_data, pos,
-                    permeability);
-        }
-        else
-        {
-            OGS_FATAL(
-                "This coupled process is not presented for "
-                "LiquidFlow process");
-        }
-    }
-}
-*/
-
 template <typename ShapeFunction, typename IntegrationMethod,
           unsigned GlobalDim>
 template <typename LaplacianGravityVelocityCalculator>
@@ -168,83 +110,6 @@ void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
     }
 }
 
-/*
-template <typename ShapeFunction, typename IntegrationMethod,
-          unsigned GlobalDim>
-template <typename LaplacianGravityVelocityCalculator>
-void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
-    assembleWithCoupledWithHeatTransport(
-        const int material_id, double const t, double const dt,
-        std::vector<double> const& local_x, std::vector<double> const& local_T0,
-        std::vector<double> const& local_T1, std::vector<double>& local_M_data,
-        std::vector<double>& local_K_data, std::vector<double>& local_b_data,
-        SpatialPosition const& pos, Eigen::MatrixXd const& permeability)
-{
-    auto const local_matrix_size = local_x.size();
-    assert(local_matrix_size == ShapeFunction::NPOINTS);
-
-    auto local_M = MathLib::createZeroedMatrix<NodalMatrixType>(
-        local_M_data, local_matrix_size, local_matrix_size);
-    auto local_K = MathLib::createZeroedMatrix<NodalMatrixType>(
-        local_K_data, local_matrix_size, local_matrix_size);
-    auto local_b = MathLib::createZeroedVector<NodalVectorType>(
-        local_b_data, local_matrix_size);
-
-    unsigned const n_integration_points =
-        _integration_method.getNumberOfPoints();
-
-    // TODO: The following two variables should be calculated inside the
-    //       the integration loop for non-constant porosity and storage models.
-    double porosity_variable = 0.;
-    double storage_variable = 0.;
-    for (unsigned ip = 0; ip < n_integration_points; ip++)
-    {
-        auto const& sm = _shape_matrices[ip];
-        auto const& wp = _integration_method.getWeightedPoint(ip);
-
-        double p = 0.;
-        NumLib::shapeFunctionInterpolate(local_x, sm.N, p);
-        double T0 = 0.;
-        NumLib::shapeFunctionInterpolate(local_T0, sm.N, T0);
-        double T = 0.;
-        NumLib::shapeFunctionInterpolate(local_T1, sm.N, T);
-
-        const double integration_factor =
-            sm.integralMeasure * sm.detJ * wp.getWeight();
-
-        // Assemble mass matrix, M
-        local_M.noalias() += _material_properties.getMassCoefficient(
-                                 material_id, t, pos, porosity_variable,
-                                 storage_variable, p, T) *
-                             sm.N.transpose() * sm.N * integration_factor;
-
-        // Compute density:
-        const double rho = _material_properties.getLiquidDensity(p, T);
-        const double rho_g = rho * _gravitational_acceleration;
-
-        // Compute viscosity:
-        const double mu = _material_properties.getViscosity(p, T);
-
-        // Assemble Laplacian, K, and RHS by the gravitational term
-        LaplacianGravityVelocityCalculator::calculateLaplacianAndGravityTerm(
-            local_K, local_b, sm, permeability, integration_factor, mu, rho_g,
-            _gravitational_axis_id);
-
-        // Add the thermal expansion term
-        auto const solid_thermal_expansion =
-            _material_properties.getSolidThermalExpansion(t, pos);
-        auto const biot_constant = _material_properties.getBiotConstant(t, pos);
-        auto const porosity = _material_properties.getPorosity(
-            material_id, t, pos, porosity_variable, T);
-        const double eff_thermal_expansion =
-            3.0 * (biot_constant - porosity) * solid_thermal_expansion -
-            porosity * _material_properties.getdLiquidDensity_dT(p, T) / rho;
-        local_b.noalias() +=
-            eff_thermal_expansion * (T - T0) * integration_factor * sm.N / dt;
-    }
-}
-*/
-
 template <typename ShapeFunction, typename IntegrationMethod,
           unsigned GlobalDim>
 std::vector<double> const&
@@ -272,41 +137,13 @@ LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
     //  the assert must be changed to perm.rows() == _element->getDimension()
     assert(permeability.rows() == GlobalDim || permeability.rows() == 1);
 
-    //if (!_coupling_term)
-    {
-        computeDarcyVelocity(permeability, local_x, veloctiy_cache_vectors);
-    }
-    /*
-    else
-    {
-        auto const local_coupled_xs =
-            getCurrentLocalSolutionsOfCoupledProcesses(
-                _coupling_term->coupled_xs, indices);
-        if (local_coupled_xs.empty())
-            computeDarcyVelocity(permeability, local_x, veloctiy_cache_vectors);
-        else
-            computeDarcyVelocityWithCoupling(permeability, local_x,
-                                             local_coupled_xs,
-                                             veloctiy_cache_vectors);
-    }*/
-
-    return veloctiy_cache;
-}
-
-template <typename ShapeFunction, typename IntegrationMethod,
-          unsigned GlobalDim>
-void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
-    computeDarcyVelocity(
-        Eigen::MatrixXd const& permeability,
-        std::vector<double> const& local_x,
-        MatrixOfVelocityAtIntegrationPoints& darcy_velocity_at_ips) const
-{
     if (permeability.size() == 1)  // isotropic or 1D problem.
         computeDarcyVelocityLocal<IsotropicCalculator>(local_x, permeability,
-                                                       darcy_velocity_at_ips);
+                                                       veloctiy_cache_vectors);
     else
-        computeDarcyVelocityLocal<AnisotropicCalculator>(local_x, permeability,
-                                                         darcy_velocity_at_ips);
+        computeDarcyVelocityLocal<AnisotropicCalculator>(
+            local_x, permeability, veloctiy_cache_vectors);
+    return veloctiy_cache;
 }
 
 template <typename ShapeFunction, typename IntegrationMethod,
@@ -346,67 +183,6 @@ void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
     }
 }
 
-/*
-template <typename ShapeFunction, typename IntegrationMethod,
-          unsigned GlobalDim>
-void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
-    computeDarcyVelocityWithCoupling(
-        Eigen::MatrixXd const& permeability, std::vector<double> const& local_x,
-        std::unordered_map<std::type_index, const std::vector<double>> const&
-            coupled_local_solutions,
-        MatrixOfVelocityAtIntegrationPoints& darcy_velocity_at_ips) const
-{
-    const auto local_T = coupled_local_solutions.at(std::type_index(
-        typeid(ProcessLib::HeatConduction::HeatConductionProcess)));
-
-    if (permeability.size() == 1)  // isotropic or 1D problem.
-        computeDarcyVelocityCoupledWithHeatTransportLocal<IsotropicCalculator>(
-            local_x, local_T, permeability, darcy_velocity_at_ips);
-    else
-        computeDarcyVelocityCoupledWithHeatTransportLocal<
-            AnisotropicCalculator>(local_x, local_T, permeability,
-                                   darcy_velocity_at_ips);
-}
-
-
-template <typename ShapeFunction, typename IntegrationMethod,
-          unsigned GlobalDim>
-template <typename LaplacianGravityVelocityCalculator>
-void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
-    computeDarcyVelocityCoupledWithHeatTransportLocal(
-        std::vector<double> const& local_x,
-        std::vector<double> const& local_T,
-        Eigen::MatrixXd const& permeability,
-        MatrixOfVelocityAtIntegrationPoints& darcy_velocity_at_ips) const
-{
-    auto const local_matrix_size = local_x.size();
-    assert(local_matrix_size == ShapeFunction::NPOINTS);
-
-    const auto local_p_vec =
-        MathLib::toVector<NodalVectorType>(local_x, local_matrix_size);
-
-    unsigned const n_integration_points =
-        _integration_method.getNumberOfPoints();
-
-    for (unsigned ip = 0; ip < n_integration_points; ip++)
-    {
-        auto const& sm = _shape_matrices[ip];
-        double p = 0.;
-        NumLib::shapeFunctionInterpolate(local_x, sm.N, p);
-        double T = 0.;
-        NumLib::shapeFunctionInterpolate(local_T, sm.N, T);
-
-        const double rho_g = _material_properties.getLiquidDensity(p, T) *
-                             _gravitational_acceleration;
-        // Compute viscosity:
-        const double mu = _material_properties.getViscosity(p, T);
-
-        LaplacianGravityVelocityCalculator::calculateVelocity(
-            ip, local_p_vec, sm, permeability, mu, rho_g,
-            _gravitational_axis_id, darcy_velocity_at_ips);
-    }
-}
-*/
 template <typename ShapeFunction, typename IntegrationMethod,
           unsigned GlobalDim>
 void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
diff --git a/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h
index b7ea92d5cfd7a2c25f8d3280d4673d1bdd9e49d5..d56ff9496d9e0b9fd9572241eac3a154652d87ed 100644
--- a/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h
+++ b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h
@@ -31,7 +31,6 @@
 
 namespace ProcessLib
 {
-struct CoupledSolutionsForStaggeredScheme;
 
 namespace LiquidFlow
 {
@@ -42,24 +41,11 @@ class LiquidFlowLocalAssemblerInterface
       public NumLib::ExtrapolatableElement
 {
 public:
-    LiquidFlowLocalAssemblerInterface() = default;
-
-    void setCoupledSolutionsForStaggeredScheme(
-        std::size_t const /*mesh_item_id*/,
-        CoupledSolutionsForStaggeredScheme* const coupled_solutions)
-    {
-        _coupled_solutions = coupled_solutions;
-    }
-
     virtual std::vector<double> const& getIntPtDarcyVelocity(
         const double /*t*/,
         GlobalVector const& /*current_solution*/,
         NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
         std::vector<double>& /*cache*/) const = 0;
-
-protected:
-    /// Pointer that is set from a Process class.
-    CoupledSolutionsForStaggeredScheme* _coupled_solutions;
 };
 
 template <typename ShapeFunction, typename IntegrationMethod,
@@ -89,8 +75,7 @@ public:
         double const gravitational_acceleration,
         double const reference_temperature,
         LiquidFlowMaterialProperties const& material_propertries)
-        : LiquidFlowLocalAssemblerInterface(),
-          _element(element),
+        : _element(element),
           _integration_method(integration_order),
           _shape_matrices(initShapeMatrices<ShapeFunction, ShapeMatricesType,
                                             IntegrationMethod, GlobalDim>(
@@ -107,13 +92,6 @@ public:
                   std::vector<double>& local_K_data,
                   std::vector<double>& local_b_data) override;
 
-    /*
-    void assembleWithCoupledTerm(
-        double const t, std::vector<double> const& local_x,
-        std::vector<double>& local_M_data, std::vector<double>& local_K_data,
-        std::vector<double>& local_b_data,
-        LocalCouplingTerm const& coupled_term) override;*/
-
     Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
         const unsigned integration_point) const override
     {
@@ -187,41 +165,12 @@ private:
                                  SpatialPosition const& pos,
                                  Eigen::MatrixXd const& permeability);
 
-    /*
-    template <typename LaplacianGravityVelocityCalculator>
-    void assembleWithCoupledWithHeatTransport(
-        const int material_id, double const t, double const dt,
-        std::vector<double> const& local_x, std::vector<double> const& local_T0,
-        std::vector<double> const& local_T1, std::vector<double>& local_M_data,
-        std::vector<double>& local_K_data, std::vector<double>& local_b_data,
-        SpatialPosition const& pos, Eigen::MatrixXd const& permeability);*/
-
-    void computeDarcyVelocity(
-        Eigen::MatrixXd const& permeability,
-        std::vector<double> const& local_x,
-        MatrixOfVelocityAtIntegrationPoints& darcy_velocity_at_ips) const;
-
-    /*
-    void computeDarcyVelocityWithCoupling(
-        Eigen::MatrixXd const& permeability, std::vector<double> const& local_x,
-        std::unordered_map<std::type_index, const std::vector<double>> const&
-            coupled_local_solutions,
-        MatrixOfVelocityAtIntegrationPoints& darcy_velocity_at_ips) const;*/
-
     template <typename LaplacianGravityVelocityCalculator>
     void computeDarcyVelocityLocal(
         std::vector<double> const& local_x,
         Eigen::MatrixXd const& permeability,
         MatrixOfVelocityAtIntegrationPoints& darcy_velocity_at_ips) const;
 
-    /*
-    template <typename LaplacianGravityVelocityCalculator>
-    void computeDarcyVelocityCoupledWithHeatTransportLocal(
-        std::vector<double> const& local_x,
-        std::vector<double> const& local_T,
-        Eigen::MatrixXd const& permeability,
-        MatrixOfVelocityAtIntegrationPoints& darcy_velocity_at_ips) const;*/
-
     const int _gravitational_axis_id;
     const double _gravitational_acceleration;
     const double _reference_temperature;
diff --git a/ProcessLib/LiquidFlow/LiquidFlowMaterialProperties.cpp b/ProcessLib/LiquidFlow/LiquidFlowMaterialProperties.cpp
index a6297371edcad2104c9f09404e4565d157c77965..d6748d70adf92ca102d42129d38103592fc6f707 100644
--- a/ProcessLib/LiquidFlow/LiquidFlowMaterialProperties.cpp
+++ b/ProcessLib/LiquidFlow/LiquidFlowMaterialProperties.cpp
@@ -18,10 +18,9 @@
 
 #include "MeshLib/PropertyVector.h"
 
-#include "ProcessLib/Parameter/Parameter.h"
 #include "ProcessLib/Parameter/SpatialPosition.h"
 
-#include "MaterialLib/Fluid/FluidProperty.h"
+#include "MaterialLib/PorousMedium/Permeability/Permeability.h"
 #include "MaterialLib/PorousMedium/Porosity/Porosity.h"
 #include "MaterialLib/PorousMedium/Storage/Storage.h"
 
@@ -54,17 +53,6 @@ double LiquidFlowMaterialProperties::getLiquidDensity(const double p,
         MaterialLib::Fluid::FluidPropertyType::Density, vars);
 }
 
-double LiquidFlowMaterialProperties::getdLiquidDensity_dT(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 _fluid_properties->getdValue(
-        MaterialLib::Fluid::FluidPropertyType::Density, vars,
-        MaterialLib::Fluid::PropertyVariableType::T);
-}
-
 double LiquidFlowMaterialProperties::getViscosity(const double p,
                                                   const double T) const
 {
@@ -75,24 +63,14 @@ double LiquidFlowMaterialProperties::getViscosity(const double p,
         MaterialLib::Fluid::FluidPropertyType::Viscosity, vars);
 }
 
-double LiquidFlowMaterialProperties::getHeatCapacity(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 _fluid_properties->getValue(
-        MaterialLib::Fluid::FluidPropertyType::HeatCapacity, vars);
-}
-
-double LiquidFlowMaterialProperties::getThermalConductivity(
-    const double p, const double T) const
+double LiquidFlowMaterialProperties::getPorosity(const int material_id,
+                                                 const double t,
+                                                 const SpatialPosition& pos,
+                                                 const double porosity_variable,
+                                                 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 _fluid_properties->getValue(
-        MaterialLib::Fluid::FluidPropertyType::ThermalConductivity, vars);
+    return _porosity_models[material_id]->getValue(t, pos, porosity_variable,
+                                                   T);
 }
 
 double LiquidFlowMaterialProperties::getMassCoefficient(
@@ -125,17 +103,5 @@ Eigen::MatrixXd const& LiquidFlowMaterialProperties::getPermeability(
                                                                  0.0);
 }
 
-double LiquidFlowMaterialProperties::getSolidThermalExpansion(
-    const double t, const SpatialPosition& pos) const
-{
-    return _solid_thermal_expansion(t, pos)[0];
-}
-
-double LiquidFlowMaterialProperties::getBiotConstant(
-    const double t, const SpatialPosition& pos) const
-{
-    return _biot_constant(t, pos)[0];
-}
-
 }  // end of namespace
 }  // end of namespace
diff --git a/ProcessLib/LiquidFlow/LiquidFlowMaterialProperties.h b/ProcessLib/LiquidFlow/LiquidFlowMaterialProperties.h
index 8145d0bd743c02a0c0e572b312d1bc204ac70393..55b21322cb1c940e1f3bb4c2bf859cdb92617c3d 100644
--- a/ProcessLib/LiquidFlow/LiquidFlowMaterialProperties.h
+++ b/ProcessLib/LiquidFlow/LiquidFlowMaterialProperties.h
@@ -19,16 +19,19 @@
 #include "MaterialLib/Fluid/FluidProperty.h"
 #include "MaterialLib/Fluid/FluidProperties/FluidProperties.h"
 
-#include "MaterialLib/PorousMedium/Porosity/Porosity.h"
-#include "MaterialLib/PorousMedium/Permeability/Permeability.h"
-#include "MaterialLib/PorousMedium/Storage/Storage.h"
-
 namespace MaterialLib
 {
 namespace Fluid
 {
 class FluidProperties;
 }
+
+namespace PorousMedium
+{
+class Permeability;
+class Porosity;
+class Storage;
+}
 }
 
 namespace BaseLib
@@ -44,9 +47,6 @@ class PropertyVector;
 
 namespace ProcessLib
 {
-template <typename T>
-struct Parameter;
-
 class SpatialPosition;
 
 namespace LiquidFlow
@@ -69,18 +69,14 @@ public:
         std::vector<std::unique_ptr<MaterialLib::PorousMedium::Storage>>&&
             storage_models,
         bool const has_material_ids,
-        MeshLib::PropertyVector<int> const& material_ids,
-        Parameter<double> const& solid_thermal_expansion,
-        Parameter<double> const& biot_constant)
+        MeshLib::PropertyVector<int> const& material_ids)
         : _has_material_ids(has_material_ids),
           _material_ids(material_ids),
           _fluid_properties(std::move(fluid_properties)),
           _intrinsic_permeability_models(
               std::move(intrinsic_permeability_models)),
           _porosity_models(std::move(porosity_models)),
-          _storage_models(std::move(storage_models)),
-          _solid_thermal_expansion(solid_thermal_expansion),
-          _biot_constant(biot_constant)
+          _storage_models(std::move(storage_models))
     {
     }
 
@@ -115,26 +111,11 @@ public:
 
     double getLiquidDensity(const double p, const double T) const;
 
-    double getdLiquidDensity_dT(const double p, const double T) const;
-
     double getViscosity(const double p, const double T) const;
 
-    double getHeatCapacity(const double p, const double T) const;
-
-    double getThermalConductivity(const double p, const double T) const;
-
     double getPorosity(const int material_id, const double t,
                        const SpatialPosition& pos,
-                       const double porosity_variable, const double T) const
-    {
-        return _porosity_models[material_id]->getValue(t, pos,
-                                                       porosity_variable, T);
-    }
-
-    double getSolidThermalExpansion(const double t,
-                                    const SpatialPosition& pos) const;
-
-    double getBiotConstant(const double t, const SpatialPosition& pos) const;
+                       const double porosity_variable, const double T) const;
 
 private:
     /// A flag to indicate whether the reference member, _material_ids,
@@ -155,8 +136,6 @@ private:
     const std::vector<std::unique_ptr<MaterialLib::PorousMedium::Storage>>
         _storage_models;
 
-    Parameter<double> const& _solid_thermal_expansion;
-    Parameter<double> const& _biot_constant;
     // Note: For the statistical data of porous media, they could be read from
     // vtu files directly. This can be done by using property vectors directly.
     // Such property vectors will be added here if they are needed.
diff --git a/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp b/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp
index 87f3b55ec607369b29f0ccd2875ac1b7670a020f..b2259f5a2d0dc069d2a93d998456e180d6eb1e22 100644
--- a/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp
+++ b/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp
@@ -20,6 +20,7 @@
 #include "LiquidFlowLocalAssembler.h"
 #include "LiquidFlowMaterialProperties.h"
 
+#include "MaterialLib/PorousMedium/Permeability/Permeability.h"
 #include "MaterialLib/PorousMedium/Porosity/Porosity.h"
 #include "MaterialLib/PorousMedium/Storage/Storage.h"
 
@@ -112,14 +113,5 @@ void LiquidFlowProcess::computeSecondaryVariableConcrete(const double t,
         _coupled_solutions);
 }
 
-void LiquidFlowProcess::setStaggeredCouplingTermToLocalAssemblers()
-{
-    DBUG("Compute the velocity for LiquidFlowProcess.");
-    /*
-    GlobalExecutor::executeMemberOnDereferenced(
-        &LiquidFlowLocalAssemblerInterface::setStaggeredCouplingTerm,
-        _local_assemblers, _coupled_solutions);*/
-}
-
 }  // end of namespace
 }  // end of namespace
diff --git a/ProcessLib/LiquidFlow/LiquidFlowProcess.h b/ProcessLib/LiquidFlow/LiquidFlowProcess.h
index 0ab5f4ac69be1c2b474d958753dd313863ac5574..a925468d379cf22c7f0b3af421b5c59cd3ab0a25 100644
--- a/ProcessLib/LiquidFlow/LiquidFlowProcess.h
+++ b/ProcessLib/LiquidFlow/LiquidFlowProcess.h
@@ -83,13 +83,6 @@ public:
         return _gravitational_acceleration;
     }
 
-    LiquidFlowMaterialProperties* getLiquidFlowMaterialProperties() const
-    {
-        return _material_properties.get();
-    }
-
-    void setStaggeredCouplingTermToLocalAssemblers() override;
-
 private:
     void initializeConcreteProcess(
         NumLib::LocalToGlobalIndexMap const& dof_table,