diff --git a/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h b/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h
index 7312f86bc8c424c0c16a23d99ec271727fe79efc..582d8f8215b90d2d235b09778fbde1fddf1286e8 100644
--- a/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h
+++ b/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h
@@ -175,6 +175,7 @@ struct ConstitutiveTempData
     FW4LWTData<DisplacementDim> fW_4_LWT;
     FW4MWpGData fW_4_MWpG;
     FW4MWpCData fW_4_MWpC;
+    FW4MWTData fW_4_MWT;
 
     using DisplacementDimVector = Eigen::Matrix<double, DisplacementDim, 1>;
     using DisplacementDimMatrix =
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveModels.h b/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveModels.h
index 4c4eec3ad793ebdb8917a052884ad87db7a532d4..01b100cbe3217e71c3fcafb22ef16fae9bdd4ab3 100644
--- a/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveModels.h
+++ b/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveModels.h
@@ -101,6 +101,7 @@ struct ConstitutiveModels
     FW4LWTModel<DisplacementDim> fW_4_LWT_model;
     FW4MWpGModel fW_4_MWpG_model;
     FW4MWpCModel fW_4_MWpC_model;
+    FW4MWTModel<DisplacementDim> fW_4_MWT_model;
 };
 }  // namespace ConstitutiveRelations
 }  // namespace ProcessLib::TH2M
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/WEquation.cpp b/ProcessLib/TH2M/ConstitutiveRelations/WEquation.cpp
index c269bad6fd81788c8dcf9b12d68c9966ea12a1b6..721c36ad78daa7128f5cef98f4d6642498eabb33 100644
--- a/ProcessLib/TH2M/ConstitutiveRelations/WEquation.cpp
+++ b/ProcessLib/TH2M/ConstitutiveRelations/WEquation.cpp
@@ -416,5 +416,27 @@ void FW4MWpCModel::eval(BiotData const& biot_data,
         (S_L - S_L_data_prev->S_L);
 }
 
+template <int DisplacementDim>
+void FW4MWTModel<DisplacementDim>::eval(
+    BiotData const& biot_data,
+    ConstituentDensityData const& constituent_density_data,
+    PorosityData const& porosity_data,
+    PureLiquidDensityData const& rho_W_LR,
+    SaturationData const& S_L_data,
+    SolidThermalExpansionData<DisplacementDim> const& s_therm_exp_data,
+    FW4MWTData& fW_4_MWT) const
+{
+    auto const S_L = S_L_data.S_L;
+    auto const S_G = 1. - S_L;
+    double const rho_W_FR =
+        S_G * constituent_density_data.rho_W_GR + S_L * rho_W_LR();
+
+    fW_4_MWT.m = -rho_W_FR * (biot_data() - porosity_data.phi) *
+                 s_therm_exp_data.beta_T_SR;
+}
+
+template struct FW4MWTModel<2>;
+template struct FW4MWTModel<3>;
+
 }  // namespace ConstitutiveRelations
 }  // namespace ProcessLib::TH2M
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/WEquation.h b/ProcessLib/TH2M/ConstitutiveRelations/WEquation.h
index ac2644843b0518ac28a3bac65bd8728c4ff48722..e2f99341be5d3af46a09d713e640cf72c236d64c 100644
--- a/ProcessLib/TH2M/ConstitutiveRelations/WEquation.h
+++ b/ProcessLib/TH2M/ConstitutiveRelations/WEquation.h
@@ -239,5 +239,26 @@ struct FW4MWpCModel
               FW4MWpCData& fW_4_MWpC) const;
 };
 
+struct FW4MWTData
+{
+    double m = nan;
+};
+
+template <int DisplacementDim>
+struct FW4MWTModel
+{
+    void eval(
+        BiotData const& biot_data,
+        ConstituentDensityData const& constituent_density_data,
+        PorosityData const& porosity_data,
+        PureLiquidDensityData const& rho_W_LR,
+        SaturationData const& S_L_data,
+        SolidThermalExpansionData<DisplacementDim> const& s_therm_exp_data,
+        FW4MWTData& fW_4_MWT) const;
+};
+
+extern template struct FW4MWTModel<2>;
+extern template struct FW4MWTModel<3>;
+
 }  // namespace ConstitutiveRelations
 }  // namespace ProcessLib::TH2M
diff --git a/ProcessLib/TH2M/TH2MFEM-impl.h b/ProcessLib/TH2M/TH2MFEM-impl.h
index c4dc4dc316cb8d2f31e641b92c8659b6fd538639..df7288e7eebf3dbdeb8eede279bdf64f16e07f54 100644
--- a/ProcessLib/TH2M/TH2MFEM-impl.h
+++ b/ProcessLib/TH2M/TH2MFEM-impl.h
@@ -411,6 +411,14 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
                                     ip_cv.beta_p_SR,
                                     ip_cv.fW_4_MWpC);
 
+        models.fW_4_MWT_model.eval(ip_cv.biot_data,
+                                   current_state.constituent_density_data,
+                                   ip_out.porosity_data,
+                                   current_state.rho_W_LR,
+                                   current_state.S_L_data,
+                                   ip_cv.s_therm_exp_data,
+                                   ip_cv.fW_4_MWT);
+
         // for variable output
         auto const xmCL = 1. - ip_out.mass_mole_fractions_data.xmWL;
 
@@ -1089,8 +1097,6 @@ void TH2MLocalAssembler<
         double const pCap = Np.dot(capillary_pressure);
         double const pCap_prev = Np.dot(capillary_pressure_prev);
 
-        double const beta_T_SR = ip_cv.s_therm_exp_data.beta_T_SR;
-
         auto const s_L = current_state.S_L_data.S_L;
         auto const s_G = 1. - s_L;
         auto const s_L_dot = (s_L - prev_state.S_L_data->S_L) / dt;
@@ -1166,8 +1172,7 @@ void TH2MLocalAssembler<
             }
         }
 
-        MWT.noalias() -= NpT * rho_W_FR * (alpha_B - ip_out.porosity_data.phi) *
-                         beta_T_SR * Np * w;
+        MWT.noalias() += NpT * ip_cv.fW_4_MWT.m * Np * w;
 
         MWu.noalias() += NpT * rho_W_FR * alpha_B * mT * Bu * w;
 
@@ -1451,7 +1456,6 @@ void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
         double const pGR_prev = Np.dot(gas_pressure_prev);
         double const pCap_prev = Np.dot(capillary_pressure_prev);
         double const T_prev = NT.dot(temperature_prev);
-        double const beta_T_SR = ip_cv.s_therm_exp_data.beta_T_SR;
 
         auto& s_L = current_state.S_L_data.S_L;
         auto const s_G = 1. - s_L;
@@ -1614,8 +1618,7 @@ void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
             }
         }
 
-        MWT.noalias() -= NpT * rho_W_FR * (alpha_B - ip_out.porosity_data.phi) *
-                         beta_T_SR * Np * w;
+        MWT.noalias() += NpT * ip_cv.fW_4_MWT.m * Np * w;
 
         MWu.noalias() += NpT * rho_W_FR * alpha_B * mT * Bu * w;