diff --git a/ProcessLib/ThermoRichardsMechanics/IntegrationPointData.h b/ProcessLib/ThermoRichardsMechanics/IntegrationPointData.h
index 89df44eb03fb6c87e428b00b03398a291a26db76..f6f388504dc59d345f0d58e5144a3d377af619a7 100644
--- a/ProcessLib/ThermoRichardsMechanics/IntegrationPointData.h
+++ b/ProcessLib/ThermoRichardsMechanics/IntegrationPointData.h
@@ -96,9 +96,8 @@ struct IntegrationPointData final
 
     MathLib::KelvinVector::KelvinMatrixType<DisplacementDim>
     computeElasticTangentStiffness(
-        double const t,
-        ParameterLib::SpatialPosition const& x_position,
-        double const dt,
+        double const t, ParameterLib::SpatialPosition const& x_position,
+        double const dt, double const temperature_prev,
         double const temperature)
     {
         namespace MPL = MaterialPropertyLib;
@@ -122,7 +121,7 @@ struct IntegrationPointData final
         variable_array_prev[static_cast<int>(MPL::Variable::mechanical_strain)]
             .emplace<KV>(KV::Zero());
         variable_array_prev[static_cast<int>(MPL::Variable::temperature)]
-            .emplace<double>(temperature);
+            .emplace<double>(temperature_prev);
 
         auto&& solution =
             solid_material.integrateStress(variable_array_prev, variable_array,
@@ -144,7 +143,7 @@ struct IntegrationPointData final
         double const t,
         ParameterLib::SpatialPosition const& x_position,
         double const dt,
-        double const temperature)
+        double const temperature_prev)
     {
         MaterialPropertyLib::VariableArray variable_array_prev;
         variable_array_prev[static_cast<int>(
@@ -157,7 +156,7 @@ struct IntegrationPointData final
                 eps_m_prev);
         variable_array_prev[static_cast<int>(
                                 MaterialPropertyLib::Variable::temperature)]
-            .emplace<double>(temperature);
+            .emplace<double>(temperature_prev);
 
         auto&& solution = solid_material.integrateStress(
             variable_array_prev, variable_array, t, x_position, dt,
diff --git a/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsFEM-impl.h b/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsFEM-impl.h
index f762e31ecf7c1ffaac7e59293639433c9a7b3792..2d4dd6350ee6f53ccc713b6cdd1de02b51d669d4 100644
--- a/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsFEM-impl.h
+++ b/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsFEM-impl.h
@@ -236,7 +236,7 @@ void ThermoRichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
         // Set eps_m_prev from potentially non-zero eps and sigma_sw from
         // restart.
         auto const C_el = ip_data_[ip].computeElasticTangentStiffness(
-            t, x_position, dt, T_ip);
+            t, x_position, dt, T_ip, T_ip);
         auto& eps = ip_data_[ip].eps;
         auto& sigma_sw = ip_data_[ip].sigma_sw;
         ip_data_[ip].eps_m_prev.noalias() =
@@ -371,6 +371,7 @@ void ThermoRichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
         NumLib::shapeFunctionInterpolate(T, N, T_ip);
         double T_dot_ip;
         NumLib::shapeFunctionInterpolate(T_dot, N, T_dot_ip);
+        double const dT = T_dot_ip * dt;
 
         double p_cap_ip;
         NumLib::shapeFunctionInterpolate(-p_L, N, p_cap_ip);
@@ -393,8 +394,9 @@ void ThermoRichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
             medium->property(MPL::PropertyType::biot_coefficient)
                 .template value<double>(variables, x_position, t, dt);
 
+        double const T_ip_prev = T_ip - dT;
         auto const C_el = ip_data_[ip].computeElasticTangentStiffness(
-            t, x_position, dt, T_ip);
+            t, x_position, dt, T_ip_prev, T_ip);
 
         auto const beta_SR =
             (1 - alpha) /
@@ -529,8 +531,7 @@ void ThermoRichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
                     .value(variables, x_position, t, dt));
 
         MathLib::KelvinVector::KelvinVectorType<DisplacementDim> const
-            dthermal_strain =
-                solid_linear_thermal_expansivity_vector * T_dot_ip * dt;
+            dthermal_strain = solid_linear_thermal_expansivity_vector * dT;
 
         eps_m.noalias() =
             solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate)
@@ -540,8 +541,8 @@ void ThermoRichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
             .emplace<MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
                 eps_m - dthermal_strain);
 
-        auto C = ip_data_[ip].updateConstitutiveRelation(variables, t,
-                                                         x_position, dt, T_ip);
+        auto C = ip_data_[ip].updateConstitutiveRelation(
+            variables, t, x_position, dt, T_ip_prev);
 
         local_Jac
             .template block<displacement_size, displacement_size>(
@@ -1144,6 +1145,7 @@ void ThermoRichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
         NumLib::shapeFunctionInterpolate(T, N, T_ip);
         double T_dot_ip;
         NumLib::shapeFunctionInterpolate(T_dot, N, T_dot_ip);
+        double const dT = T_dot_ip * dt;
 
         double p_cap_ip;
         NumLib::shapeFunctionInterpolate(-p_L, N, p_cap_ip);
@@ -1182,8 +1184,9 @@ void ThermoRichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
             medium->property(MPL::PropertyType::biot_coefficient)
                 .template value<double>(variables, x_position, t, dt);
 
+        double const T_ip_prev = T_ip - dT;
         auto const C_el = ip_data_[ip].computeElasticTangentStiffness(
-            t, x_position, dt, T_ip);
+            t, x_position, dt, T_ip_prev, T_ip);
 
         auto const beta_SR =
             (1 - alpha) /
@@ -1314,8 +1317,7 @@ void ThermoRichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
                     .value(variables, x_position, t, dt));
 
         MathLib::KelvinVector::KelvinVectorType<DisplacementDim> const
-            dthermal_strain =
-                solid_linear_thermal_expansivity_vector * T_dot_ip * dt;
+            dthermal_strain = solid_linear_thermal_expansivity_vector * dT;
 
         eps_m.noalias() =
             solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate)
@@ -1327,7 +1329,7 @@ void ThermoRichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
                 eps_m - dthermal_strain);
 
         ip_data_[ip].updateConstitutiveRelation(variables, t, x_position, dt,
-                                                T_ip);
+                                                T_ip_prev);
 
         auto const& b = process_data_.specific_body_force;