diff --git a/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h b/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h
index fa7c57299a8c73d88c5fbb38878948c622184b23..38f09a600852914be665fefb871965609833f55c 100644
--- a/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h
+++ b/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h
@@ -216,6 +216,7 @@ struct DerivativesData
     FW4LWpCDerivativeData<DisplacementDim> dfW_4_LWpC;
     FT1DerivativeData dfT_1;
     FT2DerivativeData<DisplacementDim> dfT_2;
+    FU1KUTDerivativeData<DisplacementDim> dfu_1_KuT;
     FU2KUpCDerivativeData dfu_2_KupC;
 };
 
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveModels.h b/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveModels.h
index c79d9e81fa825c3ba5890b34f674438d6fbbe7c4..277093d023f50498109c4f8c4955b90b9fb4cf79 100644
--- a/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveModels.h
+++ b/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveModels.h
@@ -116,6 +116,7 @@ struct ConstitutiveModels
     FT2Model<DisplacementDim> fT_2_model;
     FT3Model<DisplacementDim> fT_3_model;
 
+    FU1KUTModel<DisplacementDim> fu_1_KuT_model;
     FU2KUpCModel fu_2_KupC_model;
 };
 }  // namespace ConstitutiveRelations
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/UEquation.cpp b/ProcessLib/TH2M/ConstitutiveRelations/UEquation.cpp
index 07f391f74d88d7c1534a714d91756e9df518658a..dcea844e0a42ee7603f51c5099bcde21c6170cab 100644
--- a/ProcessLib/TH2M/ConstitutiveRelations/UEquation.cpp
+++ b/ProcessLib/TH2M/ConstitutiveRelations/UEquation.cpp
@@ -13,6 +13,19 @@ namespace ProcessLib::TH2M
 {
 namespace ConstitutiveRelations
 {
+template <int DisplacementDim>
+void FU1KUTModel<DisplacementDim>::dEval(
+    SolidMechanicsDataStateless<DisplacementDim> const& s_mech_data,
+    SolidThermalExpansionData<DisplacementDim> const& s_therm_exp_data,
+    FU1KUTDerivativeData<DisplacementDim>& dfu_1_KuT) const
+{
+    dfu_1_KuT.dT = s_mech_data.stiffness_tensor *
+                   s_therm_exp_data.solid_linear_thermal_expansivity_vector;
+}
+
+template struct FU1KUTModel<2>;
+template struct FU1KUTModel<3>;
+
 void FU2KUpCModel::eval(BiotData const& biot_data,
                         BishopsData const& chi_S_L,
                         FU2KUpCData& fu_2_KupC) const
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/UEquation.h b/ProcessLib/TH2M/ConstitutiveRelations/UEquation.h
index bffa79f508b57944cf18655d9a9846c29e57dbbc..70388bf50d2228018e76958c5aa2898b45dbdaa2 100644
--- a/ProcessLib/TH2M/ConstitutiveRelations/UEquation.h
+++ b/ProcessLib/TH2M/ConstitutiveRelations/UEquation.h
@@ -13,11 +13,31 @@
 #include "Biot.h"
 #include "Bishops.h"
 #include "Saturation.h"
+#include "SolidMechanics.h"
+#include "SolidThermalExpansion.h"
 
 namespace ProcessLib::TH2M
 {
 namespace ConstitutiveRelations
 {
+template <int DisplacementDim>
+struct FU1KUTDerivativeData
+{
+    KelvinVector<DisplacementDim> dT;
+};
+
+template <int DisplacementDim>
+struct FU1KUTModel
+{
+    void dEval(
+        SolidMechanicsDataStateless<DisplacementDim> const& s_mech_data,
+        SolidThermalExpansionData<DisplacementDim> const& s_therm_exp_data,
+        FU1KUTDerivativeData<DisplacementDim>& dfu_1_KuT) const;
+};
+
+extern template struct FU1KUTModel<2>;
+extern template struct FU1KUTModel<3>;
+
 struct FU2KUpCData
 {
     double m = nan;
diff --git a/ProcessLib/TH2M/TH2MFEM-impl.h b/ProcessLib/TH2M/TH2MFEM-impl.h
index 3e0beb96ac31e7d6465a447dd3307337e4d5166c..a571f7ff8ea0eb73466938d49fe98040da9cf6b3 100644
--- a/ProcessLib/TH2M/TH2MFEM-impl.h
+++ b/ProcessLib/TH2M/TH2MFEM-impl.h
@@ -489,6 +489,9 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
     updateConstitutiveVariablesDerivatives(
         Eigen::VectorXd const& local_x, Eigen::VectorXd const& local_x_prev,
         double const t, double const dt,
+        std::vector<
+            ConstitutiveRelations::ConstitutiveData<DisplacementDim>> const&
+            ip_constitutive_data,
         std::vector<
             ConstitutiveRelations::ConstitutiveTempData<DisplacementDim>> const&
             ip_constitutive_variables,
@@ -524,6 +527,7 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
     {
         auto const& ip_data = _ip_data[ip];
         auto& ip_dd = ip_d_data[ip];
+        auto const& ip_cd = ip_constitutive_data[ip];
         auto const& ip_cv = ip_constitutive_variables[ip];
         auto const& ip_out = this->output_data_[ip];
         auto const& current_state = this->current_states_[ip];
@@ -714,6 +718,9 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
             ip_cv.viscosity_data,
             ip_dd.dfT_2);
 
+        models.fu_1_KuT_model.dEval(ip_cd.s_mech_data, ip_cv.s_therm_exp_data,
+                                    ip_dd.dfu_1_KuT);
+
         models.fu_2_KupC_model.dEval(ip_cv.biot_data,
                                      ip_cv.chi_S_L,
                                      pCap_data,
@@ -1339,7 +1346,7 @@ void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
         Eigen::Map<Eigen::VectorXd const>(local_x.data(), local_x.size()),
         Eigen::Map<Eigen::VectorXd const>(local_x_prev.data(),
                                           local_x_prev.size()),
-        t, dt, ip_constitutive_variables, models);
+        t, dt, ip_constitutive_data, ip_constitutive_variables, models);
 
     for (unsigned int_point = 0; int_point < n_integration_points; int_point++)
     {
@@ -1781,11 +1788,7 @@ void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
         local_Jac
             .template block<displacement_size, temperature_size>(
                 displacement_index, temperature_index)
-            .noalias() -=
-            BuT *
-            (ip_cd.s_mech_data.stiffness_tensor *
-             ip_cv.s_therm_exp_data.solid_linear_thermal_expansivity_vector) *
-            NT * w;
+            .noalias() -= BuT * ip_dd.dfu_1_KuT.dT * NT * w;
 
         /* TODO (naumov) Test with gravity needed to check this Jacobian part.
         local_Jac
diff --git a/ProcessLib/TH2M/TH2MFEM.h b/ProcessLib/TH2M/TH2MFEM.h
index 1ea5fcb35598f3bf1149c4270a8bd5bcb30d95d2..bad9ef32fd5fc973b2f19b53b8b47deff84e9d0c 100644
--- a/ProcessLib/TH2M/TH2MFEM.h
+++ b/ProcessLib/TH2M/TH2MFEM.h
@@ -202,6 +202,9 @@ private:
     updateConstitutiveVariablesDerivatives(
         Eigen::VectorXd const& local_x, Eigen::VectorXd const& local_x_prev,
         double const t, double const dt,
+        std::vector<
+            ConstitutiveRelations::ConstitutiveData<DisplacementDim>> const&
+            ip_constitutive_data,
         std::vector<
             ConstitutiveRelations::ConstitutiveTempData<DisplacementDim>> const&
             ip_constitutive_variables,