diff --git a/ProcessLib/TH2M/ConstitutiveRelations/Base.h b/ProcessLib/TH2M/ConstitutiveRelations/Base.h
index 7ee7b82d247822c6e4579e2bd4988fc797332889..d0f26028eca002d68e9daa55417801abd26c1742 100644
--- a/ProcessLib/TH2M/ConstitutiveRelations/Base.h
+++ b/ProcessLib/TH2M/ConstitutiveRelations/Base.h
@@ -51,6 +51,8 @@ struct TemperatureData
     double T_prev;
 };
 
+using ReferenceTemperatureData =
+    BaseLib::StrongType<double, struct ReferenceTemperatureTag>;
 using GasPressureData = BaseLib::StrongType<double, struct GasPressureTag>;
 using CapillaryPressureData =
     BaseLib::StrongType<double, struct CapillaryPressureTag>;
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/SolidThermalExpansion.cpp b/ProcessLib/TH2M/ConstitutiveRelations/SolidThermalExpansion.cpp
index 39ed38dd266a09329305829b90301a04d6759880..196d01b7904748f5c235d4f271d0399569d44004 100644
--- a/ProcessLib/TH2M/ConstitutiveRelations/SolidThermalExpansion.cpp
+++ b/ProcessLib/TH2M/ConstitutiveRelations/SolidThermalExpansion.cpp
@@ -18,6 +18,7 @@ namespace ConstitutiveRelations
 template <int DisplacementDim>
 void SolidThermalExpansionModel<DisplacementDim>::eval(
     SpaceTimeData const& x_t, MediaData const& media_data,
+    TemperatureData const& T_data, ReferenceTemperatureData T0,
     SolidThermalExpansionData<DisplacementDim>& out) const
 {
     namespace MPL = MaterialPropertyLib;
@@ -34,6 +35,9 @@ void SolidThermalExpansionModel<DisplacementDim>::eval(
     using Invariants = MathLib::KelvinVector::Invariants<KelvinVectorSize>;
     out.beta_T_SR =
         Invariants::trace(out.solid_linear_thermal_expansivity_vector);
+
+    double const delta_T = T_data.T - T0();
+    out.thermal_volume_strain = out.beta_T_SR * delta_T;
 }
 
 template struct SolidThermalExpansionModel<2>;
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/SolidThermalExpansion.h b/ProcessLib/TH2M/ConstitutiveRelations/SolidThermalExpansion.h
index 6b7d32057121bd6149b6a83e991ec6aa14cb414d..9f225f30f7d2691173ea2bdde69fd1f136db19af 100644
--- a/ProcessLib/TH2M/ConstitutiveRelations/SolidThermalExpansion.h
+++ b/ProcessLib/TH2M/ConstitutiveRelations/SolidThermalExpansion.h
@@ -21,12 +21,14 @@ struct SolidThermalExpansionData
     KelvinVector<DisplacementDim> solid_linear_thermal_expansivity_vector;
     double beta_T_SR;  /// Isotropic solid phase volumetric thermal expansion
                        /// coefficient.
+    double thermal_volume_strain = nan;
 };
 
 template <int DisplacementDim>
 struct SolidThermalExpansionModel
 {
     void eval(SpaceTimeData const& x_t, MediaData const& media_data,
+              TemperatureData const& T_data, ReferenceTemperatureData T0,
               SolidThermalExpansionData<DisplacementDim>& out) const;
 };
 
diff --git a/ProcessLib/TH2M/IntegrationPointData.h b/ProcessLib/TH2M/IntegrationPointData.h
index 288d0bbc109bb0eddf2816173c6f94eec410f710..cd79de8021a87d6767dceacde4ed059aff771145 100644
--- a/ProcessLib/TH2M/IntegrationPointData.h
+++ b/ProcessLib/TH2M/IntegrationPointData.h
@@ -64,8 +64,6 @@ struct IntegrationPointData final
     GlobalDimVectorType w_GS;
     GlobalDimVectorType w_LS;
 
-    double thermal_volume_strain = std::numeric_limits<double>::quiet_NaN();
-
     double integration_weight = std::numeric_limits<double>::quiet_NaN();
 
     void pushBackState() { rho_u_eff_prev = rho_u_eff; }
diff --git a/ProcessLib/TH2M/TH2MFEM-impl.h b/ProcessLib/TH2M/TH2MFEM-impl.h
index 860c34b15f0560bbd2dfe841047a14da1a302e0f..3550b5f609fd6de896115040e81ed3d262cfc401 100644
--- a/ProcessLib/TH2M/TH2MFEM-impl.h
+++ b/ProcessLib/TH2M/TH2MFEM-impl.h
@@ -160,6 +160,8 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
         ConstitutiveRelations::TemperatureData const T_data{T, T_prev};
         ConstitutiveRelations::GasPressureData const pGR_data{pGR};
         ConstitutiveRelations::CapillaryPressureData const pCap_data{pCap};
+        ConstitutiveRelations::ReferenceTemperatureData const T0{
+            this->process_data_.reference_temperature(t, pos)[0]};
         double const pLR = pGR - pCap;
         GlobalDimVectorType const gradpGR = gradNp * gas_pressure;
         GlobalDimVectorType const gradpCap = gradNp * capillary_pressure;
@@ -196,7 +198,7 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
             current_state.swelling_data, ip_cv.swelling_data);
 
         // solid phase linear thermal expansion coefficient
-        models.s_therm_exp_model.eval({pos, t, dt}, media_data,
+        models.s_therm_exp_model.eval({pos, t, dt}, media_data, T_data, T0,
                                       ip_cv.s_therm_exp_data);
 
         models.mechanical_strain_model.eval(
@@ -248,11 +250,6 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
                 .template value<double>(
                     vars, pos, t, std::numeric_limits<double>::quiet_NaN());
 
-        double const T0 = this->process_data_.reference_temperature(t, pos)[0];
-        double const delta_T(T - T0);
-        ip_data.thermal_volume_strain =
-            ip_cv.s_therm_exp_data.beta_T_SR * delta_T;
-
         // initial porosity
         auto const phi_0 = medium.property(MPL::PropertyType::porosity)
                                .template value<double>(vars, pos, t, dt);
@@ -263,8 +260,9 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
         auto const& m = Invariants::identity2;
         double const div_u = m.transpose() * eps;
 
-        const double phi_S = phi_S_0 * (1. + ip_data.thermal_volume_strain -
-                                        ip_cv.biot_data() * div_u);
+        const double phi_S =
+            phi_S_0 * (1. + ip_cv.s_therm_exp_data.thermal_volume_strain -
+                       ip_cv.biot_data() * div_u);
 #else   // NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
         const double phi_S = phi_S_0;
 #endif  // NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
@@ -275,8 +273,9 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
 
         // solid phase density
 #ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
-        auto const rhoSR = rho_ref_SR * (1. - ip_data.thermal_volume_strain +
-                                         (ip_cv.biot_data() - 1.) * div_u);
+        auto const rhoSR =
+            rho_ref_SR * (1. - ip_cv.s_therm_exp_data.thermal_volume_strain +
+                          (ip_cv.biot_data() - 1.) * div_u);
 #else   // NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
         auto const rhoSR = rho_ref_SR;
 #endif  // NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
@@ -364,7 +363,7 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
                     .template dValue<double>(vars, MPL::Variable::temperature,
                                              pos, t, dt)
 #ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
-                * (1. - ip_data.thermal_volume_strain +
+                * (1. - ip_cv.s_therm_exp_data.thermal_volume_strain +
                    (ip_cv.biot_data() - 1.) * div_u) -
             rho_ref_SR * ip_cv.s_therm_exp_data.beta_T_SR
 #endif
@@ -376,11 +375,12 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
                 vars, MPL::Variable::temperature, pos, t, dt);
 
         auto const dphi_S_0_dT = -dphi_0_dT;
-        const double dphi_S_dT = dphi_S_0_dT
+        const double dphi_S_dT =
+            dphi_S_0_dT
 #ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
-                                     * (1. + ip_data.thermal_volume_strain -
-                                        ip_cv.biot_data() * div_u) +
-                                 phi_S_0 * ip_cv.s_therm_exp_data.beta_T_SR
+                * (1. + ip_cv.s_therm_exp_data.thermal_volume_strain -
+                   ip_cv.biot_data() * div_u) +
+            phi_S_0 * ip_cv.s_therm_exp_data.beta_T_SR
 #endif
             ;