diff --git a/MaterialLib/MPL/Properties/PermeabilityMohrCoulombFailureIndexModel.cpp b/MaterialLib/MPL/Properties/PermeabilityMohrCoulombFailureIndexModel.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..f04b9217993b42aabbc70477c35c3e5e5be4760e
--- /dev/null
+++ b/MaterialLib/MPL/Properties/PermeabilityMohrCoulombFailureIndexModel.cpp
@@ -0,0 +1,143 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2020, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ * Created on June 4, 2020, 10:13 AM
+ */
+
+#include "PermeabilityMohrCoulombFailureIndexModel.h"
+
+#include <boost/math/constants/constants.hpp>
+#include <cmath>
+#include <limits>
+
+#include "BaseLib/Error.h"
+#include "MaterialLib/MPL/Medium.h"
+#include "MaterialLib/MPL/Utils/FormEigenTensor.h"
+#include "MaterialLib/MPL/Utils/GetSymmetricTensor.h"
+#include "MathLib/KelvinVector.h"
+#include "MathLib/MathTools.h"
+#include "ParameterLib/CoordinateSystem.h"
+#include "ParameterLib/Parameter.h"
+
+namespace MaterialPropertyLib
+{
+template <int DisplacementDim>
+PermeabilityMohrCoulombFailureIndexModel<DisplacementDim>::
+ PermeabilityMohrCoulombFailureIndexModel(
+ std::string name, ParameterLib::Parameter<double> const& k0,
+ double const kr, double const b, double const c, double const phi,
+ double const k_max,
+ ParameterLib::CoordinateSystem const* const local_coordinate_system)
+ : k0_(k0),
+ kr_(kr),
+ b_(b),
+ c_(c),
+ phi_(boost::math::constants::degree<double>() * phi),
+ k_max_(k_max),
+ local_coordinate_system_(local_coordinate_system)
+{
+ name_ = std::move(name);
+}
+
+template <int DisplacementDim>
+void PermeabilityMohrCoulombFailureIndexModel<DisplacementDim>::checkScale()
+ const
+{
+ if (!std::holds_alternative<Medium*>(scale_))
+ {
+ OGS_FATAL(
+ "The property 'PermeabilityMohrCoulombFailureIndexModel' is "
+ "implemented on the 'medium' scale only.");
+ }
+}
+
+template <int DisplacementDim>
+PropertyDataType
+PermeabilityMohrCoulombFailureIndexModel<DisplacementDim>::value(
+ VariableArray const& variable_array,
+ ParameterLib::SpatialPosition const& pos, double const t,
+ double const /*dt*/) const
+{
+ auto const& stress_vector = std::get<SymmetricTensor<DisplacementDim>>(
+ variable_array[static_cast<int>(Variable::stress)]);
+
+ auto const& stress_tensor =
+ formEigenTensor<3>(static_cast<PropertyDataType>(stress_vector));
+
+ Eigen::SelfAdjointEigenSolver<Eigen::Matrix<double, 3, 3>>
+ eigenvalue_solver(stress_tensor);
+
+ // Principle stress
+ auto const sigma = eigenvalue_solver.eigenvalues();
+
+ auto k_data = k0_(t, pos);
+
+ double const max_sigma = std::max(std::fabs(sigma[0]), std::fabs(sigma[2]));
+
+ if (max_sigma < std::numeric_limits<double>::epsilon())
+ {
+ return fromVector(k_data);
+ }
+
+ double const sigma_m = 0.5 * (sigma[2] + sigma[0]);
+ double tau_tt = c_ * std::cos(phi_) - sigma_m * std::sin(phi_);
+ const double apex_cut_offset = 0.001;
+ if (std::fabs(tau_tt) < apex_cut_offset)
+ {
+ tau_tt = apex_cut_offset * c_ * std::cos(phi_);
+ }
+
+ // +- tau_t = tau_tt
+ double const f = 0.5 * std::fabs(sigma[2] - sigma[0]) / tau_tt;
+
+ if (f >= 1.0)
+ {
+ const double exp_value = std::exp(b_ * f);
+ for (auto& k_i : k_data)
+ {
+ k_i = std::min(k_i + kr_ * exp_value, k_max_);
+ }
+ }
+
+ // Local coordinate transformation is only applied for the case that the
+ // initial intrinsic permeability is given with orthotropic assumption.
+ if (local_coordinate_system_ && (k_data.size() == DisplacementDim))
+ {
+ Eigen::Matrix<double, DisplacementDim, DisplacementDim> const e =
+ local_coordinate_system_->transformation<DisplacementDim>(pos);
+ Eigen::Matrix<double, DisplacementDim, DisplacementDim> k =
+ Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Zero();
+
+ for (int i = 0; i < DisplacementDim; ++i)
+ {
+ Eigen::Matrix<double, DisplacementDim, DisplacementDim> const
+ ei_otimes_ei = e.col(i) * e.col(i).transpose();
+
+ k += k_data[i] * ei_otimes_ei;
+ }
+ return k;
+ }
+
+ return fromVector(k_data);
+}
+
+template <int DisplacementDim>
+PropertyDataType
+PermeabilityMohrCoulombFailureIndexModel<DisplacementDim>::dValue(
+ VariableArray const& /*variable_array*/, Variable const /*variable*/,
+ ParameterLib::SpatialPosition const& /*pos*/, double const /*t*/,
+ double const /*dt*/) const
+{
+ OGS_FATAL(
+ "The derivative of the intrinsic permeability k(sigma, ...) with "
+ "respect to stress tensor (sigma) is not implemented because that "
+ "dk/du is normally omitted.");
+}
+template class PermeabilityMohrCoulombFailureIndexModel<2>;
+template class PermeabilityMohrCoulombFailureIndexModel<3>;
+} // namespace MaterialPropertyLib
diff --git a/MaterialLib/MPL/Properties/PermeabilityMohrCoulombFailureIndexModel.h b/MaterialLib/MPL/Properties/PermeabilityMohrCoulombFailureIndexModel.h
new file mode 100644
index 0000000000000000000000000000000000000000..2a03b63660416d5ee2b0cc003a0206515e919f89
--- /dev/null
+++ b/MaterialLib/MPL/Properties/PermeabilityMohrCoulombFailureIndexModel.h
@@ -0,0 +1,98 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2020, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ * Created on June 4, 2020, 10:13 AM
+ */
+
+#pragma once
+
+#include "MaterialLib/MPL/Property.h"
+#include "MaterialLib/MPL/VariableType.h"
+
+namespace ParameterLib
+{
+struct CoordinateSystem;
+template <typename T>
+struct Parameter;
+} // namespace ParameterLib
+
+namespace MaterialPropertyLib
+{
+/**
+ * \brief A failure index dependent permeability model \cite wangetalPerm2020
+ *
+ * \f[ \mathbf{k} =
+ * \mathbf{k}_0+ H(f-1) k_\text{r} \mathrm{e}^{b f}\mathbf{I}\f]
+ *
+ * where
+ * \f$\mathbf{k}_0\f$ is the intrinsic permeability
+ * of the undamaged material,
+ * \f$H\f$ is the Heaviside step function, \f$f\f$ is the failure index,
+ * \f$k_\text{r}\f$ is a reference permeability,
+ * \f$b\f$ is a fitting parameter.
+ * \f$k_\text{r}\f$ and \f$b\f$ can be calibrated by experimental data.
+ * The failure index \f$f\f$ is calculated from the Mohr Coulomb failure
+ * criterion comparing an acting shear stress. The Mohr Coulomb failure
+ * criterion \cite labuz2012mohr takes
+ * the form
+ * \f[\tau(\sigma)=c-\sigma \mathrm{tan} \phi\f]
+ * with \f$\tau\f$ the shear stress, \f$c\f$ the cohesion, \f$\sigma\f$ the
+ * tensile stress, and \f$\phi\f$ the internal friction angle.
+ *
+ * The failure index is calculated by
+ * \f[
+ * f=\frac{\tau_m }{\cos(\phi)\tau(\sigma_m)}
+ * \f]
+ * with
+ * \f$\tau_m=(\sigma_3-\sigma_1)/2\f$
+ * and \f$\sigma_m=(\sigma_1+\sigma_3)/2\f$,
+ * where \f$\sigma_1\f$ and \f$\sigma_3\f$ are the minimum and maximum shear
+ * stress, respectively.
+ *
+ * Note: the conventional mechanics notations are used, which mean that tensile
+ * stress is positive.
+ *
+ */
+template <int DisplacementDim>
+class PermeabilityMohrCoulombFailureIndexModel final : public Property
+{
+public:
+ PermeabilityMohrCoulombFailureIndexModel(
+ std::string name, ParameterLib::Parameter<double> const& k0,
+ double const kr, double const b, double const c, double const phi,
+ double const k_max,
+ ParameterLib::CoordinateSystem const* const local_coordinate_system);
+
+ void checkScale() const override;
+
+ PropertyDataType value(VariableArray const& variable_array,
+ ParameterLib::SpatialPosition const& pos,
+ double const t, double const dt) const override;
+ PropertyDataType dValue(VariableArray const& variable_array,
+ Variable const variable,
+ ParameterLib::SpatialPosition const& pos,
+ double const t, double const dt) const override;
+
+private:
+ /// Intrinsic permeability for undamaged material.
+ ParameterLib::Parameter<double> const& k0_;
+ /// Reference permeability.
+ double const kr_;
+ /// Fitting parameter.
+ double const b_;
+ /// Cohesion.
+ double const c_;
+ /// Angle of internal friction.
+ double const phi_;
+
+ /// Maximum permeability.
+ double const k_max_;
+ ParameterLib::CoordinateSystem const* const local_coordinate_system_;
+};
+
+} // namespace MaterialPropertyLib