diff --git a/ProcessLib/TH2M/TH2MFEM-impl.h b/ProcessLib/TH2M/TH2MFEM-impl.h
index d7763bb9713851f88159e1d015e94619d4dbed3a..115b76ed8692b39d7799718cb80dda7de3448016 100644
--- a/ProcessLib/TH2M/TH2MFEM-impl.h
+++ b/ProcessLib/TH2M/TH2MFEM-impl.h
@@ -182,7 +182,41 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
             medium.property(MPL::PropertyType::permeability)
                 .value(vars, pos, t, dt));
 
+        auto const Bu =
+            LinearBMatrix::computeBMatrix<DisplacementDim,
+                                          ShapeFunctionDisplacement::NPOINTS,
+                                          typename BMatricesType::BMatrixType>(
+                gradNu, Nu, x_coord, _is_axially_symmetric);
+
+        auto& eps = ip_data.eps;
+        eps.noalias() = Bu * displacement;
+
         // relative permeability
+        // Set mechanical variables for the intrinsic permeability model
+        // For stress dependent permeability.
+        {
+            // Note: if Bishop model is available, ip_data.s_L in the following
+            // computation should be replaced with the Bishop value.
+            auto const sigma_total =
+                (_ip_data[ip].sigma_eff - ip_data.alpha_B *
+                                              (pGR - ip_data.s_L * pCap) *
+                                              Invariants::identity2)
+                    .eval();
+
+            vars[static_cast<int>(MPL::Variable::total_stress)]
+                .emplace<SymmetricTensor>(
+                    MathLib::KelvinVector::kelvinVectorToSymmetricTensor(
+                        sigma_total));
+        }
+        // For strain dependent permeability
+        vars[static_cast<int>(MPL::Variable::volumetric_strain)] =
+            Invariants::trace(eps);
+        vars[static_cast<int>(MPL::Variable::equivalent_plastic_strain)] =
+            _ip_data[ip].material_state_variables->getEquivalentPlasticStrain();
+        vars[static_cast<int>(MPL::Variable::mechanical_strain)]
+            .emplace<MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
+                eps);
+
         ip_data.k_rel_G =
             medium
                 .property(
@@ -211,15 +245,6 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
         MathLib::KelvinVector::KelvinVectorType<DisplacementDim> const
             dthermal_strain = ip_data.alpha_T_SR * T_dot * dt;
 
-        auto const Bu =
-            LinearBMatrix::computeBMatrix<DisplacementDim,
-                                          ShapeFunctionDisplacement::NPOINTS,
-                                          typename BMatricesType::BMatrixType>(
-                gradNu, Nu, x_coord, _is_axially_symmetric);
-
-        auto& eps = ip_data.eps;
-        eps.noalias() = Bu * displacement;
-
         auto& eps_prev = ip_data.eps_prev;
         auto& eps_m = ip_data.eps_m;
         auto& eps_m_prev = ip_data.eps_m_prev;
diff --git a/ProcessLib/TH2M/TH2MFEM.h b/ProcessLib/TH2M/TH2MFEM.h
index b770cb1b83824aeb76203aa24009b6b299871bbf..a62a3267368d48c7617a143dfd0149d640eb6386 100644
--- a/ProcessLib/TH2M/TH2MFEM.h
+++ b/ProcessLib/TH2M/TH2MFEM.h
@@ -70,6 +70,8 @@ public:
 
     static int const KelvinVectorSize =
         MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim);
+    using SymmetricTensor = Eigen::Matrix<double, KelvinVectorSize, 1>;
+
     using Invariants = MathLib::KelvinVector::Invariants<KelvinVectorSize>;
 
     TH2MLocalAssembler(TH2MLocalAssembler const&) = delete;