diff --git a/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h b/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h
index 78a648ac6c2e2723b83992c4148fa05b76b673ff..6c7a0cab1afdc22105e988d1345d7a6b69584c18 100644
--- a/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h
+++ b/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h
@@ -173,6 +173,7 @@ struct ConstitutiveTempData
     FW4LWpGData<DisplacementDim> fW_4_LWpG;
     FW4LWpCData<DisplacementDim> fW_4_LWpC;
     FW4LWTData<DisplacementDim> fW_4_LWT;
+    FW4MWpGData fW_4_MWpG;
 
     using DisplacementDimVector = Eigen::Matrix<double, DisplacementDim, 1>;
     using DisplacementDimMatrix =
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveModels.h b/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveModels.h
index 4acaa507183a921a5ead0d9f34521c292f5f1463..c5e975961a25d041cfb68592fee00271a76f4776 100644
--- a/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveModels.h
+++ b/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveModels.h
@@ -99,6 +99,7 @@ struct ConstitutiveModels
     FW4LWpGModel<DisplacementDim> fW_4_LWpG_model;
     FW4LWpCModel<DisplacementDim> fW_4_LWpC_model;
     FW4LWTModel<DisplacementDim> fW_4_LWT_model;
+    FW4MWpGModel fW_4_MWpG_model;
 };
 }  // namespace ConstitutiveRelations
 }  // namespace ProcessLib::TH2M
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/WEquation.cpp b/ProcessLib/TH2M/ConstitutiveRelations/WEquation.cpp
index 44b2b7acdb0748864379613958e3533062ed90ba..65056b82bf96f944b8344bef40d9f83b9e86d15d 100644
--- a/ProcessLib/TH2M/ConstitutiveRelations/WEquation.cpp
+++ b/ProcessLib/TH2M/ConstitutiveRelations/WEquation.cpp
@@ -375,5 +375,22 @@ void FW4LWTModel<DisplacementDim>::eval(
 template struct FW4LWTModel<2>;
 template struct FW4LWTModel<3>;
 
+void FW4MWpGModel::eval(BiotData const& biot_data,
+                        ConstituentDensityData const& constituent_density_data,
+                        PorosityData const& porosity_data,
+                        PureLiquidDensityData const& rho_W_LR,
+                        SaturationData const& S_L_data,
+                        SolidCompressibilityData const& beta_p_SR,
+                        FW4MWpGData& fW_4_MWpG) const
+{
+    double const S_L = S_L_data.S_L;
+    double 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_MWpG.m = rho_W_FR * (biot_data() - porosity_data.phi) * beta_p_SR();
+}
+
 }  // namespace ConstitutiveRelations
 }  // namespace ProcessLib::TH2M
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/WEquation.h b/ProcessLib/TH2M/ConstitutiveRelations/WEquation.h
index e4bccfd1705c6e4772de99f1fbe9d81ab86bb942..a574fda04c7554068a88c69fcd44c7b82a758408 100644
--- a/ProcessLib/TH2M/ConstitutiveRelations/WEquation.h
+++ b/ProcessLib/TH2M/ConstitutiveRelations/WEquation.h
@@ -204,5 +204,21 @@ struct FW4LWTModel
               FW4LWTData<DisplacementDim>& fW_4_LWT) const;
 };
 
+struct FW4MWpGData
+{
+    double m = nan;
+};
+
+struct FW4MWpGModel
+{
+    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,
+              SolidCompressibilityData const& beta_p_SR,
+              FW4MWpGData& fW_4_MWpG) const;
+};
+
 }  // namespace ConstitutiveRelations
 }  // namespace ProcessLib::TH2M
diff --git a/ProcessLib/TH2M/TH2MFEM-impl.h b/ProcessLib/TH2M/TH2MFEM-impl.h
index 944598803b64a9f901af72f75e1d25e5722d0d26..fe608037f394e90186ee05bd54ebadbfbdf5649d 100644
--- a/ProcessLib/TH2M/TH2MFEM-impl.h
+++ b/ProcessLib/TH2M/TH2MFEM-impl.h
@@ -393,6 +393,14 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
                                    current_state.S_L_data,
                                    ip_cv.fW_4_LWT);
 
+        models.fW_4_MWpG_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.beta_p_SR,
+                                    ip_cv.fW_4_MWpG);
+
         // for variable output
         auto const xmCL = 1. - ip_out.mass_mole_fractions_data.xmWL;
 
@@ -1137,9 +1145,7 @@ void TH2MLocalAssembler<
         // W-component equation
         // ---------------------------------------------------------------------
 
-        MWpG.noalias() += NpT * rho_W_FR *
-                          (alpha_B - ip_out.porosity_data.phi) * beta_p_SR *
-                          Np * w;
+        MWpG.noalias() += NpT * ip_cv.fW_4_MWpG.m * Np * w;
         MWpC.noalias() -= NpT * rho_W_FR *
                           (alpha_B - ip_out.porosity_data.phi) * beta_p_SR *
                           s_L * Np * w;
@@ -1596,9 +1602,7 @@ void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
         // W-component equation
         // ---------------------------------------------------------------------
 
-        MWpG.noalias() += NpT * rho_W_FR *
-                          (alpha_B - ip_out.porosity_data.phi) * beta_p_SR *
-                          Np * w;
+        MWpG.noalias() += NpT * ip_cv.fW_4_MWpG.m * Np * w;
         MWpC.noalias() -= NpT * rho_W_FR *
                           (alpha_B - ip_out.porosity_data.phi) * beta_p_SR *
                           s_L * Np * w;