diff --git a/Documentation/ProjectFile/prj/processes/process/RICHARDS_MECHANICS/t_intrinsic_permeability.md b/Documentation/ProjectFile/prj/processes/process/RICHARDS_MECHANICS/t_intrinsic_permeability.md
deleted file mode 100644
index f945fe4d68fe60bacc884d5d1504e0ad52e4f6c7..0000000000000000000000000000000000000000
--- a/Documentation/ProjectFile/prj/processes/process/RICHARDS_MECHANICS/t_intrinsic_permeability.md
+++ /dev/null
@@ -1 +0,0 @@
-\copydoc ProcessLib::RichardsMechanics::RichardsMechanicsProcessData::intrinsic_permeability
diff --git a/ProcessLib/Parameter/ConstantParameter.cpp b/ProcessLib/Parameter/ConstantParameter.cpp
index 94d7f6929f3ee40645cfc28392f38f1de3c3245b..12ac93c33fe0cfe88112ef1264fcf08a26632a92 100644
--- a/ProcessLib/Parameter/ConstantParameter.cpp
+++ b/ProcessLib/Parameter/ConstantParameter.cpp
@@ -25,12 +25,20 @@ std::unique_ptr<ParameterBase> createConstantParameter(
     // then required tag 'values'.
     {
         //! \ogs_file_param{prj__parameters__parameter__Constant__value}
-        auto const value = config.getConfigParameterOptional<double>("value");
+        auto const value =
+            config.getConfigParameterOptional<std::vector<double>>("value");
 
         if (value)
         {
-            DBUG("Using value %g for constant parameter.", *value);
-            return std::make_unique<ConstantParameter<double>>(name, *value);
+            if (value->size() != 1)
+            {
+                OGS_FATAL(
+                    "Expected to read exactly one value, but %d were given.",
+                    value->size());
+            }
+            DBUG("Using value %g for constant parameter.", (*value)[0]);
+            return std::make_unique<ConstantParameter<double>>(name,
+                                                               (*value)[0]);
         }
     }
 
diff --git a/ProcessLib/RichardsMechanics/CreateRichardsMechanicsProcess.cpp b/ProcessLib/RichardsMechanics/CreateRichardsMechanicsProcess.cpp
index 40fe159072d5b2284b36c0fc66f4add87757dbb5..44a939c5534a60e2ca3f3d7ff0a4e8be6174b451 100644
--- a/ProcessLib/RichardsMechanics/CreateRichardsMechanicsProcess.cpp
+++ b/ProcessLib/RichardsMechanics/CreateRichardsMechanicsProcess.cpp
@@ -105,15 +105,6 @@ std::unique_ptr<Process> createRichardsMechanicsProcess(
         MaterialLib::Solids::createConstitutiveRelations<DisplacementDim>(
             parameters, config);
 
-    // Intrinsic permeability
-    auto& intrinsic_permeability = findParameter<double>(
-        config,
-        //! \ogs_file_param_special{prj__processes__process__RICHARDS_MECHANICS__intrinsic_permeability}
-        "intrinsic_permeability", parameters, 1);
-
-    DBUG("Use '%s' as intrinsic conductivity parameter.",
-         intrinsic_permeability.name.c_str());
-
     // Fluid bulk modulus
     auto& fluid_bulk_modulus = findParameter<double>(
         config,
@@ -197,7 +188,6 @@ std::unique_ptr<Process> createRichardsMechanicsProcess(
         material_ids,
         std::move(flow_material),
         std::move(solid_constitutive_relations),
-        intrinsic_permeability,
         fluid_bulk_modulus,
         biot_coefficient,
         solid_density,
diff --git a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h
index 250df5ed328f776953e2d48ce266547534bf471a..dc0f58a397f7e664a4285c1fafeb6da13c9b14a1 100644
--- a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h
+++ b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h
@@ -13,6 +13,8 @@
 
 #include "RichardsMechanicsFEM.h"
 
+#include <cassert>
+
 #include "MaterialLib/SolidModels/SelectSolidConstitutiveRelation.h"
 #include "MathLib/KelvinVector.h"
 #include "NumLib/Function/Interpolation.h"
@@ -22,6 +24,29 @@ namespace ProcessLib
 {
 namespace RichardsMechanics
 {
+template <int DisplacementDim>
+Eigen::Matrix<double, DisplacementDim, DisplacementDim> intrinsicPermeability(
+    double const t, SpatialPosition const& x_position, int const material_id,
+    RichardsFlow::RichardsFlowMaterialProperties const& material)
+{
+    const Eigen::MatrixXd& permeability =
+        material.getPermeability(material_id, t, x_position, DisplacementDim);
+    if (permeability.rows() == DisplacementDim)
+    {
+        return permeability;
+    }
+    if (permeability.rows() == 1)
+    {
+        return Eigen::Matrix<double, DisplacementDim,
+                             DisplacementDim>::Identity() *
+               permeability(0, 0);
+    }
+
+    OGS_FATAL(
+        "Intrinsic permeability dimension is neither %d nor one, but %dx%d.",
+        DisplacementDim, permeability.rows(), permeability.cols());
+}
+
 template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
           typename IntegrationMethod, int DisplacementDim>
 RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
@@ -191,12 +216,14 @@ void RichardsMechanicsLocalAssembler<
         double const k_rel =
             _process_data.flow_material->getRelativePermeability(
                 t, x_position, -p_cap_ip, temperature, S_L);
-        auto const mu = _process_data.flow_material->getFluidViscosity(
+
+        double const mu = _process_data.flow_material->getFluidViscosity(
             -p_cap_ip, temperature);
 
-        double const K_intrinsic =
-            _process_data.intrinsic_permeability(t, x_position)[0];
-        double const K_over_mu = K_intrinsic / mu;
+        GlobalDimMatrixType const rho_K_over_mu =
+            intrinsicPermeability<DisplacementDim>(
+                t, x_position, material_id, *_process_data.flow_material) *
+            (rho_LR * k_rel / mu);
 
         //
         // displacement equation, displacement part
@@ -231,11 +258,10 @@ void RichardsMechanicsLocalAssembler<
 
         K.template block<pressure_size, pressure_size>(pressure_index,
                                                        pressure_index)
-            .noalias() +=
-            dNdx_p.transpose() * rho_LR * k_rel * K_over_mu * dNdx_p * w;
+            .noalias() += dNdx_p.transpose() * rho_K_over_mu * dNdx_p * w;
 
         rhs.template segment<pressure_size>(pressure_index).noalias() +=
-            dNdx_p.transpose() * rho_LR * rho_LR * k_rel * K_over_mu * b * w;
+            dNdx_p.transpose() * rho_LR * rho_K_over_mu * b * w;
 
         //
         // displacement equation, pressure part
@@ -401,12 +427,15 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
         double const k_rel =
             _process_data.flow_material->getRelativePermeability(
                 t, x_position, -p_cap_ip, temperature, S_L);
-        auto const mu = _process_data.flow_material->getFluidViscosity(
+
+        double const mu = _process_data.flow_material->getFluidViscosity(
             -p_cap_ip, temperature);
 
-        double const K_intrinsic =
-            _process_data.intrinsic_permeability(t, x_position)[0];
-        double const K_over_mu = K_intrinsic / mu;
+        GlobalDimMatrixType const rho_Ki_over_mu =
+            intrinsicPermeability<DisplacementDim>(
+                t, x_position, material_id, *_process_data.flow_material) *
+            (rho_LR / mu);
+
         //
         // displacement equation, displacement part
         //
@@ -461,7 +490,7 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
         // pressure equation, pressure part.
         //
         laplace_p.noalias() +=
-            dNdx_p.transpose() * rho_LR * k_rel * K_over_mu * dNdx_p * w;
+            dNdx_p.transpose() * k_rel * rho_Ki_over_mu * dNdx_p * w;
 
         double const a0 = (alpha - porosity) / K_SR;
         double const specific_storage =
@@ -508,17 +537,17 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
         local_Jac
             .template block<pressure_size, pressure_size>(pressure_index,
                                                           pressure_index)
-            .noalias() += dNdx_p.transpose() * rho_LR * K_over_mu * grad_p_cap *
+            .noalias() += dNdx_p.transpose() * rho_Ki_over_mu * grad_p_cap *
                           dk_rel_dS_l * dS_L_dp_cap * N_p * w;
 
         local_Jac
             .template block<pressure_size, pressure_size>(pressure_index,
                                                           pressure_index)
-            .noalias() += dNdx_p.transpose() * rho_LR * rho_LR * K_over_mu * b *
+            .noalias() += dNdx_p.transpose() * rho_LR * rho_Ki_over_mu * b *
                           dk_rel_dS_l * dS_L_dp_cap * N_p * w;
 
         local_rhs.template segment<pressure_size>(pressure_index).noalias() +=
-            dNdx_p.transpose() * rho_LR * rho_LR * k_rel * K_over_mu * b * w;
+            dNdx_p.transpose() * rho_LR * k_rel * rho_Ki_over_mu * b * w;
     }
 
     // pressure equation, pressure part.
@@ -626,6 +655,9 @@ std::vector<double> const& RichardsMechanicsLocalAssembler<
             pressure_size> const>(local_x.data() + pressure_index,
                                   pressure_size);
 
+    auto const material_id =
+        _process_data.flow_material->getMaterialID(_element.getID());
+
     unsigned const n_integration_points =
         _integration_method.getNumberOfPoints();
 
@@ -640,10 +672,11 @@ std::vector<double> const& RichardsMechanicsLocalAssembler<
         NumLib::shapeFunctionInterpolate(-p_L, N_p, p_cap_ip);
 
         auto const temperature = _process_data.temperature(t, x_position)[0];
-        auto const mu = _process_data.flow_material->getFluidViscosity(
-            -p_cap_ip, temperature);
-        double const K_over_mu =
-            _process_data.intrinsic_permeability(t, x_position)[0] / mu;
+        GlobalDimMatrixType const K_over_mu =
+            intrinsicPermeability<DisplacementDim>(
+                t, x_position, material_id, *_process_data.flow_material) /
+            _process_data.flow_material->getFluidViscosity(-p_cap_ip,
+                                                           temperature);
         auto const rho_LR = _process_data.flow_material->getFluidDensity(
             -p_cap_ip, temperature);
         auto const& b = _process_data.specific_body_force;
@@ -651,7 +684,7 @@ std::vector<double> const& RichardsMechanicsLocalAssembler<
         // Compute the velocity
         auto const& dNdx_p = _ip_data[ip].dNdx_p;
         cache_matrix.col(ip).noalias() =
-            -K_over_mu * dNdx_p * p_L - K_over_mu * rho_LR * b;
+            -K_over_mu * dNdx_p * p_L - rho_LR * K_over_mu * b;
     }
 
     return cache;
diff --git a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM.h b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM.h
index 756705788fb24e340edd595a79956d2f6afb7838..1632ab1a4831e243818819a2166d166054c87cd1 100644
--- a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM.h
+++ b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM.h
@@ -49,6 +49,9 @@ public:
     using ShapeMatricesTypePressure =
         ShapeMatrixPolicyType<ShapeFunctionPressure, DisplacementDim>;
 
+    using GlobalDimMatrixType =
+        typename ShapeMatricesTypePressure::GlobalDimMatrixType;
+
     static int const KelvinVectorSize =
         MathLib::KelvinVector::KelvinVectorDimensions<DisplacementDim>::value;
     using Invariants = MathLib::KelvinVector::Invariants<KelvinVectorSize>;
diff --git a/ProcessLib/RichardsMechanics/RichardsMechanicsProcessData.h b/ProcessLib/RichardsMechanics/RichardsMechanicsProcessData.h
index 986309cdf133f608abbb02624a89a8cf2da8687c..6b22b03b1de2a6dc1845f7407a3df33e7a148327 100644
--- a/ProcessLib/RichardsMechanics/RichardsMechanicsProcessData.h
+++ b/ProcessLib/RichardsMechanics/RichardsMechanicsProcessData.h
@@ -42,7 +42,6 @@ struct RichardsMechanicsProcessData
                  std::unique_ptr<
                      MaterialLib::Solids::MechanicsBase<DisplacementDim>>>&&
             solid_materials_,
-        Parameter<double> const& intrinsic_permeability_,
         Parameter<double> const& fluid_bulk_modulus_,
         Parameter<double> const& biot_coefficient_,
         Parameter<double> const& solid_density_,
@@ -53,7 +52,6 @@ struct RichardsMechanicsProcessData
         : material_ids(material_ids_),
           flow_material{std::move(flow_material_)},
           solid_materials{std::move(solid_materials_)},
-          intrinsic_permeability(intrinsic_permeability_),
           fluid_bulk_modulus(fluid_bulk_modulus_),
           biot_coefficient(biot_coefficient_),
           solid_density(solid_density_),
@@ -83,8 +81,6 @@ struct RichardsMechanicsProcessData
         int,
         std::unique_ptr<MaterialLib::Solids::MechanicsBase<DisplacementDim>>>
         solid_materials;
-    /// Permeability of the solid. A scalar quantity, Parameter<double>.
-    Parameter<double> const& intrinsic_permeability;
     /// Fluid's bulk modulus. A scalar quantity, Parameter<double>.
     Parameter<double> const& fluid_bulk_modulus;
     /// Biot coefficient. A scalar quantity, Parameter<double>.
diff --git a/ProcessLib/RichardsMechanics/Tests.cmake b/ProcessLib/RichardsMechanics/Tests.cmake
index e121a9815d8e281eb40934ddab7fde792950d603..c8de028c6f9d20c8b9b8c1aa41dfc5a961694841 100644
--- a/ProcessLib/RichardsMechanics/Tests.cmake
+++ b/ProcessLib/RichardsMechanics/Tests.cmake
@@ -72,6 +72,24 @@ AddTest(
     GLOB flow_fully_saturated_pcs_0_ts_*.vtu NodalForces NodalForces 1e-15 0
 )
 
+AddTest(
+    NAME RichardsMechanics_square_1e2_flow_fully_saturated_anisotropic
+    PATH RichardsMechanics
+    EXECUTABLE ogs
+    EXECUTABLE_ARGS flow_fully_saturated_anisotropic.prj
+    WRAPPER time
+    TESTER vtkdiff
+    REQUIREMENTS NOT OGS_USE_MPI
+    DIFF_DATA
+    GLOB flow_fully_saturated_anisotropic_pcs_0_ts_*.vtu displacement displacement 2e-14 0
+    GLOB flow_fully_saturated_anisotropic_pcs_0_ts_*.vtu sigma sigma 1e-14 0
+    GLOB flow_fully_saturated_anisotropic_pcs_0_ts_*.vtu epsilon epsilon 1e-15 0
+    GLOB flow_fully_saturated_anisotropic_pcs_0_ts_*.vtu pressure pressure 1e-15 1e-15
+    GLOB flow_fully_saturated_anisotropic_pcs_0_ts_*.vtu velocity velocity 1e-15 1e-15
+    GLOB flow_fully_saturated_anisotropic_pcs_0_ts_*.vtu HydraulicFlow HydraulicFlow 1e-15 0
+    GLOB flow_fully_saturated_anisotropic_pcs_0_ts_*.vtu NodalForces NodalForces 1e-15 0
+)
+
 AddTest(
     NAME RichardsMechanics_RichardsFlow_2d_small
     PATH RichardsMechanics
diff --git a/Tests/Data/RichardsMechanics/RichardsFlow_2d_quasinewton.prj b/Tests/Data/RichardsMechanics/RichardsFlow_2d_quasinewton.prj
index 39d5eb4bf1dcfed4dbd4dbb220895ce9931b7526..e92ab00575c2406c49e3bda627ec3f20dba12e06 100644
--- a/Tests/Data/RichardsMechanics/RichardsFlow_2d_quasinewton.prj
+++ b/Tests/Data/RichardsMechanics/RichardsFlow_2d_quasinewton.prj
@@ -18,7 +18,6 @@
                 <youngs_modulus>E</youngs_modulus>
                 <poissons_ratio>nu</poissons_ratio>
             </constitutive_relation>
-            <intrinsic_permeability>k</intrinsic_permeability>
             <solid_bulk_modulus>K_SR</solid_bulk_modulus>
             <fluid_bulk_modulus>K_LR</fluid_bulk_modulus>
             <biot_coefficient>alpha</biot_coefficient>
diff --git a/Tests/Data/RichardsMechanics/RichardsFlow_2d_small.prj b/Tests/Data/RichardsMechanics/RichardsFlow_2d_small.prj
index 76f932e81a12203bf0d47fbab9502ff7bce69b97..3faa5c9cdad4d8b7015bffebeb6e2e68a3ea1d76 100644
--- a/Tests/Data/RichardsMechanics/RichardsFlow_2d_small.prj
+++ b/Tests/Data/RichardsMechanics/RichardsFlow_2d_small.prj
@@ -13,7 +13,6 @@
                 <youngs_modulus>E</youngs_modulus>
                 <poissons_ratio>nu</poissons_ratio>
             </constitutive_relation>
-            <intrinsic_permeability>k</intrinsic_permeability>
             <solid_bulk_modulus>K_SR</solid_bulk_modulus>
             <fluid_bulk_modulus>K_LR</fluid_bulk_modulus>
             <biot_coefficient>alpha</biot_coefficient>
diff --git a/Tests/Data/RichardsMechanics/confined_compression_fully_saturated.prj b/Tests/Data/RichardsMechanics/confined_compression_fully_saturated.prj
index 803d2011607593dc89e1cd7b91bcc8f5da92cb93..829d6feb2a7b83a9f8540fd57ba8a16c22ca2208 100644
--- a/Tests/Data/RichardsMechanics/confined_compression_fully_saturated.prj
+++ b/Tests/Data/RichardsMechanics/confined_compression_fully_saturated.prj
@@ -20,7 +20,6 @@
                 <youngs_modulus>E</youngs_modulus>
                 <poissons_ratio>nu</poissons_ratio>
             </constitutive_relation>
-            <intrinsic_permeability>k</intrinsic_permeability>
             <solid_bulk_modulus>K_SR</solid_bulk_modulus>
             <fluid_bulk_modulus>K_LR</fluid_bulk_modulus>
             <biot_coefficient>alpha</biot_coefficient>
diff --git a/Tests/Data/RichardsMechanics/flow_fully_saturated.prj b/Tests/Data/RichardsMechanics/flow_fully_saturated.prj
index 0de1a5157f567c7b9fa23cca84d7680cebae2d15..f919f42a3c8285910d92a6ce00bbaa94da99b651 100644
--- a/Tests/Data/RichardsMechanics/flow_fully_saturated.prj
+++ b/Tests/Data/RichardsMechanics/flow_fully_saturated.prj
@@ -13,7 +13,6 @@
                 <youngs_modulus>E</youngs_modulus>
                 <poissons_ratio>nu</poissons_ratio>
             </constitutive_relation>
-            <intrinsic_permeability>k</intrinsic_permeability>
             <solid_bulk_modulus>K_SR</solid_bulk_modulus>
             <fluid_bulk_modulus>K_LR</fluid_bulk_modulus>
             <biot_coefficient>alpha</biot_coefficient>
diff --git a/Tests/Data/RichardsMechanics/flow_fully_saturated_anisotropic.prj b/Tests/Data/RichardsMechanics/flow_fully_saturated_anisotropic.prj
new file mode 100644
index 0000000000000000000000000000000000000000..341afe8a88fb812c00b4f0cc8cfe7494311467f7
--- /dev/null
+++ b/Tests/Data/RichardsMechanics/flow_fully_saturated_anisotropic.prj
@@ -0,0 +1,276 @@
+<?xml version="1.0" encoding="ISO-8859-1"?>
+<OpenGeoSysProject>
+    <mesh>square_1x1_quad8_1e2.vtu</mesh>
+    <geometry>square_1x1.gml</geometry>
+    <processes>
+        <process>
+            <name>RM</name>
+            <type>RICHARDS_MECHANICS</type>
+            <integration_order>3</integration_order>
+            <dimension>2</dimension>
+            <constitutive_relation>
+                <type>LinearElasticIsotropic</type>
+                <youngs_modulus>E</youngs_modulus>
+                <poissons_ratio>nu</poissons_ratio>
+            </constitutive_relation>
+            <solid_bulk_modulus>K_SR</solid_bulk_modulus>
+            <fluid_bulk_modulus>K_LR</fluid_bulk_modulus>
+            <biot_coefficient>alpha</biot_coefficient>
+            <solid_density>rho_sr</solid_density>
+            <process_variables>
+                <displacement>displacement</displacement>
+                <pressure>pressure</pressure>
+            </process_variables>
+            <material_property>
+                <fluid>
+                    <density>
+                        <type>Constant</type>
+                        <value>1</value>
+                    </density>
+                    <viscosity>
+                        <type>Constant</type>
+                        <value>1e-9</value>
+                    </viscosity>
+                </fluid>
+                <porous_medium>
+                    <porous_medium id="0">
+                        <permeability>
+                            <permeability_tensor_entries>k</permeability_tensor_entries>
+                            <type>Constant</type>
+                        </permeability>
+                        <porosity>
+                            <type>Constant</type>
+                            <porosity_parameter>phi</porosity_parameter>
+                        </porosity>
+                        <storage>
+                            <type>Constant</type>
+                            <value>1e+300</value>
+                        </storage>
+                        <capillary_pressure>
+                            <type>Curve</type>
+                            <curve>
+                                <coords>1 0.98684211 0.97368421 0.96052632 0.94894737 0.92105263 0.89473684 0.86842105 0.84210526 0.81578947 0.78947368 0.65789474 0.63157895 0.57894737 0.52631579 0.39473684 0.26315789
+                                </coords>
+                                <values>1421.99634 1707.55369 2050.45507 2462.21599 2892.29884 3232.80969 3590.71602 3988.24638 4429.78757 4920.21206 5464.93174 9238.24221 10261.0136 12658.7874 15616.8683 26399.6732 44627.5611
+                                </values>
+                            </curve>
+                        </capillary_pressure>
+                        <relative_permeability>
+                            <type>Curve</type>
+                            <curve>
+                                <coords>0.26315789 0.39473684 0.52631579 0.57894737 0.63157895 0.65789474 0.78947368 0.81578947 0.84210526 0.86842105 0.89473684 0.92105263 0.94894737 0.96052632 0.97368421 0.98684211 1
+                                    </coords>
+                                <values>4.43E-05 0.00026547 0.00159003 0.00325358 0.00665757 0.00952343 0.05704014 0.08159396 0.11671736 0.16696017 0.23883078 0.34163922 0.49931406 0.58449912 0.69907308 0.8361059 1
+                                    </values>
+                            </curve>
+                        </relative_permeability>
+                    </porous_medium>
+                </porous_medium>
+            </material_property>
+            <secondary_variables>
+                <secondary_variable type="static" internal_name="sigma" output_name="sigma"/>
+                <secondary_variable type="static" internal_name="epsilon" output_name="epsilon"/>
+                <secondary_variable type="static" internal_name="velocity" output_name="velocity"/>
+                <secondary_variable type="static" internal_name="saturation" output_name="saturation"/>
+            </secondary_variables>
+            <specific_body_force>0 0</specific_body_force>
+            <temperature> temp </temperature>
+        </process>
+    </processes>
+    <time_loop>
+        <processes>
+            <process ref="RM">
+                <nonlinear_solver>basic_newton</nonlinear_solver>
+                <convergence_criterion>
+                    <type>PerComponentDeltaX</type>
+                    <norm_type>NORM2</norm_type>
+                    <abstols>1e-14 1e-16 1e-16</abstols>
+                </convergence_criterion>
+                <time_discretization>
+                    <type>BackwardEuler</type>
+                </time_discretization>
+                <time_stepping>
+                    <type>FixedTimeStepping</type>
+                    <t_initial>0</t_initial>
+                    <t_end>2</t_end>
+                    <timesteps>
+                        <pair>
+                            <repeat>2</repeat>
+                            <delta_t>1</delta_t>
+                        </pair>
+                    </timesteps>
+                </time_stepping>
+            </process>
+        </processes>
+        <output>
+            <type>VTK</type>
+            <prefix>flow_fully_saturated_anisotropic</prefix>
+            <timesteps>
+                <pair>
+                    <repeat>10000</repeat>
+                    <each_steps>1</each_steps>
+                </pair>
+            </timesteps>
+            <variables>
+                <variable>displacement</variable>
+                <variable>pressure</variable>
+                <variable>sigma</variable>
+                <variable>epsilon</variable>
+                <variable>velocity</variable>
+                <variable>saturation</variable>
+            </variables>
+        </output>
+    </time_loop>
+    <parameters>
+        <!-- Mechanics -->
+        <parameter>
+            <name>K_SR</name>
+            <type>Constant</type>
+            <value>1e100</value>
+        </parameter>
+        <parameter>
+            <name>E</name>
+            <type>Constant</type>
+            <value>1e10</value>
+        </parameter>
+        <parameter>
+            <name>nu</name>
+            <type>Constant</type>
+            <value>.3</value>
+        </parameter>
+        <!-- Model parameters -->
+        <parameter>
+            <name>k</name>
+            <type>Constant</type>
+            <values>7.75e-13 3.89711e-13 3.89711e-13 3.25e-13</values>
+        </parameter>
+        <parameter>
+            <name>mu</name>
+            <type>Constant</type>
+            <value>1e-9</value>
+        </parameter>
+        <parameter>
+            <name>alpha</name>
+            <type>Constant</type>
+            <value>1</value>
+        </parameter>
+        <parameter>
+            <name>phi</name>
+            <type>Constant</type>
+            <value>0.8</value>
+        </parameter>
+        <parameter>
+            <name>rho_sr</name>
+            <type>Constant</type>
+            <value>1</value>
+        </parameter>
+        <parameter>
+            <name>K_LR</name>
+            <type>Constant</type>
+            <value>1e100</value>
+        </parameter>
+        <parameter>
+            <name>temp</name>
+            <type>Constant</type>
+            <value>293.15</value>
+        </parameter>
+        <parameter>
+            <name>displacement0</name>
+            <type>Constant</type>
+            <values>0 0</values>
+        </parameter>
+        <parameter>
+            <name>pressure_ic</name>
+            <type>Constant</type>
+            <values>0</values>
+        </parameter>
+        <parameter>
+            <name>dirichlet0</name>
+            <type>Constant</type>
+            <value>0</value>
+        </parameter>
+        <parameter>
+            <name>dirichlet1</name>
+            <type>Constant</type>
+            <value>1</value>
+        </parameter>
+    </parameters>
+    <process_variables>
+        <process_variable>
+            <name>displacement</name>
+            <components>2</components>
+            <order>2</order>
+            <initial_condition>displacement0</initial_condition>
+            <boundary_conditions>
+                <boundary_condition>
+                    <geometrical_set>square_1x1_geometry</geometrical_set>
+                    <geometry>left</geometry>
+                    <type>Dirichlet</type>
+                    <component>0</component>
+                    <parameter>dirichlet0</parameter>
+                </boundary_condition>
+                <boundary_condition>
+                    <geometrical_set>square_1x1_geometry</geometrical_set>
+                    <geometry>right</geometry>
+                    <type>Dirichlet</type>
+                    <component>0</component>
+                    <parameter>dirichlet0</parameter>
+                </boundary_condition>
+                <boundary_condition>
+                    <geometrical_set>square_1x1_geometry</geometrical_set>
+                    <geometry>bottom</geometry>
+                    <type>Dirichlet</type>
+                    <component>1</component>
+                    <parameter>dirichlet0</parameter>
+                </boundary_condition>
+                <boundary_condition>
+                    <geometrical_set>square_1x1_geometry</geometrical_set>
+                    <geometry>top</geometry>
+                    <type>Dirichlet</type>
+                    <component>1</component>
+                    <parameter>dirichlet0</parameter>
+                </boundary_condition>
+            </boundary_conditions>
+        </process_variable>
+        <process_variable>
+            <name>pressure</name>
+            <components>1</components>
+            <order>1</order>
+            <initial_condition>pressure_ic</initial_condition>
+            <boundary_conditions>
+                <boundary_condition>
+                    <geometrical_set>square_1x1_geometry</geometrical_set>
+                    <geometry>left</geometry>
+                    <type>Dirichlet</type>
+                    <parameter>dirichlet0</parameter>
+                </boundary_condition>
+                <boundary_condition>
+                    <geometrical_set>square_1x1_geometry</geometrical_set>
+                    <geometry>right</geometry>
+                    <type>Dirichlet</type>
+                    <parameter>dirichlet1</parameter>
+                </boundary_condition>
+            </boundary_conditions>
+        </process_variable>
+    </process_variables>
+    <nonlinear_solvers>
+        <nonlinear_solver>
+            <name>basic_newton</name>
+            <type>Newton</type>
+            <max_iter>50</max_iter>
+            <linear_solver>general_linear_solver</linear_solver>
+        </nonlinear_solver>
+    </nonlinear_solvers>
+    <linear_solvers>
+        <linear_solver>
+            <name>general_linear_solver</name>
+            <eigen>
+                <solver_type>SparseLU</solver_type>
+                <scaling>true</scaling>
+                <!--
+                <solver_type>PardisoLU</solver_type>
+                -->
+            </eigen>
+        </linear_solver>
+    </linear_solvers>
+</OpenGeoSysProject>
diff --git a/Tests/Data/RichardsMechanics/flow_fully_saturated_anisotropic_pcs_0_ts_1_t_1.000000.vtu b/Tests/Data/RichardsMechanics/flow_fully_saturated_anisotropic_pcs_0_ts_1_t_1.000000.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..7fafcc63019aec2a1c065fe5bd12c9131f1fb393
--- /dev/null
+++ b/Tests/Data/RichardsMechanics/flow_fully_saturated_anisotropic_pcs_0_ts_1_t_1.000000.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:60c64026c77113bb842cfbd9d5a629a4a2fe60c46a74583c2f95d69c793536f5
+size 53706
diff --git a/Tests/Data/RichardsMechanics/flow_fully_saturated_anisotropic_pcs_0_ts_2_t_2.000000.vtu b/Tests/Data/RichardsMechanics/flow_fully_saturated_anisotropic_pcs_0_ts_2_t_2.000000.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..d22b5547c5c93aba17069c7a6e4cc46b98e4ec9e
--- /dev/null
+++ b/Tests/Data/RichardsMechanics/flow_fully_saturated_anisotropic_pcs_0_ts_2_t_2.000000.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:7be5b3dcbaf0c8b1ac025f43c1af92fc50beb5101e37c26f3207910fa34820e9
+size 53668
diff --git a/Tests/Data/RichardsMechanics/gravity.prj b/Tests/Data/RichardsMechanics/gravity.prj
index 1871179477170087e99626a8e0938386b3b7ece5..3e492acb986dc113fd3eac96603874453a9d56b2 100644
--- a/Tests/Data/RichardsMechanics/gravity.prj
+++ b/Tests/Data/RichardsMechanics/gravity.prj
@@ -20,7 +20,6 @@
                 <youngs_modulus>E</youngs_modulus>
                 <poissons_ratio>nu</poissons_ratio>
             </constitutive_relation>
-            <intrinsic_permeability>k</intrinsic_permeability>
             <solid_bulk_modulus>K_SR</solid_bulk_modulus>
             <fluid_bulk_modulus>K_LR</fluid_bulk_modulus>
             <biot_coefficient>alpha</biot_coefficient>
@@ -43,7 +42,7 @@
                 <porous_medium>
                     <porous_medium id="0">
                         <permeability>
-                            <permeability_tensor_entries>kappa1</permeability_tensor_entries>
+                            <permeability_tensor_entries>k</permeability_tensor_entries>
                             <type>Constant</type>
                         </permeability>
                         <porosity>
@@ -200,11 +199,6 @@
             <type>Constant</type>
             <value>293.15</value>
         </parameter>
-        <parameter>
-            <name>kappa1</name>
-            <type>Constant</type>
-            <values>4.46e-13</values>
-        </parameter>
         <parameter>
             <name>displacement0</name>
             <type>Constant</type>
diff --git a/Tests/Data/RichardsMechanics/mechanics_linear.prj b/Tests/Data/RichardsMechanics/mechanics_linear.prj
index e38e3d1fe87787aaaa76bd086690432b07866301..74dfa7758729a16dedec7a434c589f26f3bcb6f3 100644
--- a/Tests/Data/RichardsMechanics/mechanics_linear.prj
+++ b/Tests/Data/RichardsMechanics/mechanics_linear.prj
@@ -26,7 +26,6 @@
                 <youngs_modulus>E</youngs_modulus>
                 <poissons_ratio>nu</poissons_ratio>
             </constitutive_relation>
-            <intrinsic_permeability>k</intrinsic_permeability>
             <solid_bulk_modulus>K_SR</solid_bulk_modulus>
             <fluid_bulk_modulus>K_LR</fluid_bulk_modulus>
             <biot_coefficient>alpha</biot_coefficient>