diff --git a/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h b/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h
index 2cd3860282d84a67feef5e0595dbba2076823ffa..bbb0a0426790c6a3cb42ce3f360519cc8737164d 100644
--- a/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h
+++ b/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h
@@ -149,7 +149,6 @@ struct ConstitutiveTempData
     ElasticTangentStiffnessData<DisplacementDim> C_el_data;
     BiotData biot_data;
     SolidCompressibilityData beta_p_SR;
-    SaturationDataDeriv dS_L_dp_cap;
     BishopsData chi_S_L;
     SolidThermalExpansionData<DisplacementDim> s_therm_exp_data;
     TotalStressData<DisplacementDim> total_stress_data;
@@ -196,6 +195,7 @@ struct ConstitutiveTempData
 template <int DisplacementDim>
 struct DerivativesData
 {
+    SaturationDataDeriv dS_L_dp_cap;
     AdvectionDerivativeData<DisplacementDim> advection_d_data;
     ThermalConductivityDerivativeData<DisplacementDim>
         thermal_conductivity_d_data;
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/Saturation.cpp b/ProcessLib/TH2M/ConstitutiveRelations/Saturation.cpp
index 5f4693174a6e9b7f4cd0c88f2a5c9b3efbfd67ce..8ca73adaa1b7820caccfd10db4d4fb64f0fee18e 100644
--- a/ProcessLib/TH2M/ConstitutiveRelations/Saturation.cpp
+++ b/ProcessLib/TH2M/ConstitutiveRelations/Saturation.cpp
@@ -16,8 +16,7 @@ namespace ConstitutiveRelations
 void SaturationModel::eval(SpaceTimeData const& x_t,
                            MediaData const& media_data,
                            CapillaryPressureData const& p_cap,
-                           SaturationData& S_L_data,
-                           SaturationDataDeriv& dS_L_data) const
+                           SaturationData& S_L_data) const
 {
     namespace MPL = MaterialPropertyLib;
     MPL::VariableArray variables;
@@ -27,6 +26,18 @@ void SaturationModel::eval(SpaceTimeData const& x_t,
 
     S_L_data.S_L = medium.property(MPL::PropertyType::saturation)
                        .template value<double>(variables, x_t.x, x_t.t, x_t.dt);
+}
+
+void SaturationModel::dEval(SpaceTimeData const& x_t,
+                            MediaData const& media_data,
+                            CapillaryPressureData const& p_cap,
+                            SaturationDataDeriv& dS_L_data) const
+{
+    namespace MPL = MaterialPropertyLib;
+    MPL::VariableArray variables;
+    variables.capillary_pressure = p_cap();
+
+    auto const& medium = media_data.medium;
 
     dS_L_data() = medium.property(MPL::PropertyType::saturation)
                       .template dValue<double>(
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/Saturation.h b/ProcessLib/TH2M/ConstitutiveRelations/Saturation.h
index f8ab3071673e39bce24b1a67305bdc3551f8db17..66dbbea5dcc486d838923da8cc2130a3c8a13773 100644
--- a/ProcessLib/TH2M/ConstitutiveRelations/Saturation.h
+++ b/ProcessLib/TH2M/ConstitutiveRelations/Saturation.h
@@ -34,8 +34,12 @@ using SaturationDataDeriv =
 struct SaturationModel
 {
     void eval(SpaceTimeData const& x_t, MediaData const& media_data,
-              CapillaryPressureData const& p_cap, SaturationData& S_L_data,
-              SaturationDataDeriv& dS_L_data) const;
+              CapillaryPressureData const& p_cap,
+              SaturationData& S_L_data) const;
+
+    void dEval(SpaceTimeData const& x_t, MediaData const& media_data,
+               CapillaryPressureData const& p_cap,
+               SaturationDataDeriv& dS_L_data) const;
 };
 }  // namespace ConstitutiveRelations
 }  // namespace ProcessLib::TH2M
diff --git a/ProcessLib/TH2M/TH2MFEM-impl.h b/ProcessLib/TH2M/TH2MFEM-impl.h
index 4eff6a4203ecc3ee0f3b141c0e602d4b2ab8fc6a..bfadf134ba24c5f7dd06b89974fec33010c4f3a1 100644
--- a/ProcessLib/TH2M/TH2MFEM-impl.h
+++ b/ProcessLib/TH2M/TH2MFEM-impl.h
@@ -180,7 +180,7 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
 
         ip_out.eps_data.eps.noalias() = Bu * displacement;
         models.S_L_model.eval({pos, t, dt}, media_data, pCap_data,
-                              current_state.S_L_data, ip_cv.dS_L_dp_cap);
+                              current_state.S_L_data);
 
         models.chi_S_L_model.eval({pos, t, dt}, media_data,
                                   current_state.S_L_data, ip_cv.chi_S_L);
@@ -546,16 +546,19 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
         ConstitutiveRelations::CapillaryPressureData const pCap_data{
             Np.dot(capillary_pressure)};
 
+        models.S_L_model.dEval({pos, t, dt}, media_data, pCap_data,
+                               ip_dd.dS_L_dp_cap);
+
         models.advection_model.dEval(current_state.constituent_density_data,
                                      ip_out.permeability_data,
                                      ip_cv.viscosity_data,
-                                     ip_cv.dS_L_dp_cap,
+                                     ip_dd.dS_L_dp_cap,
                                      ip_cv.phase_transition_data,
                                      ip_dd.advection_d_data);
 
         models.thermal_conductivity_model.dEval(
             {pos, t, dt}, media_data, T_data, ip_out.porosity_data,
-            ip_cv.porosity_d_data, current_state.S_L_data, ip_cv.dS_L_dp_cap,
+            ip_cv.porosity_d_data, current_state.S_L_data, ip_dd.dS_L_dp_cap,
             ip_dd.thermal_conductivity_d_data);
 
         models.solid_density_model.dEval({pos, t, dt}, media_data, T_data,
@@ -571,7 +574,7 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
             ip_out.porosity_data,
             ip_cv.porosity_d_data,
             current_state.S_L_data,
-            ip_cv.dS_L_dp_cap,
+            ip_dd.dS_L_dp_cap,
             ip_out.solid_density_data,
             ip_dd.solid_density_d_data,
             ip_out.solid_enthalpy_data,
@@ -585,7 +588,7 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
             ip_out.porosity_data,
             ip_cv.porosity_d_data,
             current_state.S_L_data,
-            ip_cv.dS_L_dp_cap,
+            ip_dd.dS_L_dp_cap,
             ip_out.solid_density_data,
             ip_dd.solid_density_d_data,
             ip_out.solid_enthalpy_data,
@@ -600,7 +603,7 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
                                      ip_out.porosity_data,
                                      ip_cv.porosity_d_data,
                                      current_state.S_L_data,
-                                     ip_cv.dS_L_dp_cap,
+                                     ip_dd.dS_L_dp_cap,
                                      ip_cv.beta_p_SR,
                                      ip_dd.dfC_2a);
         }
@@ -609,7 +612,7 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
                                  prev_state.constituent_density_data,
                                  ip_cv.phase_transition_data,
                                  current_state.S_L_data,
-                                 ip_cv.dS_L_dp_cap,
+                                 ip_dd.dS_L_dp_cap,
                                  ip_dd.dfC_3a);
 
         models.fC_4_LCpG_model.dEval(ip_out.permeability_data,
@@ -621,7 +624,7 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
         models.fC_4_LCpC_model.dEval(current_state.constituent_density_data,
                                      ip_out.permeability_data,
                                      ip_cv.phase_transition_data,
-                                     ip_cv.dS_L_dp_cap,
+                                     ip_dd.dS_L_dp_cap,
                                      ip_cv.viscosity_data,
                                      ip_dd.dfC_4_LCpC);
 
@@ -658,7 +661,7 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
                                     ip_cv.porosity_d_data,
                                     current_state.rho_W_LR,
                                     current_state.S_L_data,
-                                    ip_cv.dS_L_dp_cap,
+                                    ip_dd.dS_L_dp_cap,
                                     ip_cv.beta_p_SR,
                                     ip_dd.dfW_2);
         }
@@ -670,14 +673,14 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
                                  prev_state.rho_W_LR,
                                  current_state.rho_W_LR,
                                  current_state.S_L_data,
-                                 ip_cv.dS_L_dp_cap,
+                                 ip_dd.dS_L_dp_cap,
                                  ip_dd.dfW_3a);
 
         models.fW_4_LWpG_model.dEval(current_state.constituent_density_data,
                                      ip_out.permeability_data,
                                      ip_cv.phase_transition_data,
                                      current_state.rho_W_LR,
-                                     ip_cv.dS_L_dp_cap,
+                                     ip_dd.dS_L_dp_cap,
                                      ip_cv.viscosity_data,
                                      ip_dd.dfW_4_LWpG);
 
@@ -688,7 +691,7 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
                                      ip_out.porosity_data,
                                      current_state.rho_W_LR,
                                      current_state.S_L_data,
-                                     ip_cv.dS_L_dp_cap,
+                                     ip_dd.dS_L_dp_cap,
                                      ip_cv.viscosity_data,
                                      ip_dd.dfW_4_LWpC);
 
@@ -709,7 +712,7 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
         models.fu_2_KupC_model.dEval(ip_cv.biot_data,
                                      ip_cv.chi_S_L,
                                      pCap_data,
-                                     ip_cv.dS_L_dp_cap,
+                                     ip_dd.dS_L_dp_cap,
                                      ip_dd.dfu_2_KupC);
     }
 
@@ -1491,7 +1494,7 @@ void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
             local_Jac.template block<C_size, W_size>(C_index, W_index)
                 .noalias() += NpT *
                               (ip_dd.dfC_2a.dp_cap * s_L_dot +
-                               ip_cv.fC_2a.a * ip_cv.dS_L_dp_cap() / dt) *
+                               ip_cv.fC_2a.a * ip_dd.dS_L_dp_cap() / dt) *
                               Np * w;
 
             local_Jac
@@ -1587,7 +1590,7 @@ void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
             local_Jac.template block<W_size, W_size>(W_index, W_index)
                 .noalias() += NpT *
                               (ip_dd.dfW_2.dp_cap * s_L_dot +
-                               ip_cv.fW_2.a * ip_cv.dS_L_dp_cap() / dt) *
+                               ip_cv.fW_2.a * ip_dd.dS_L_dp_cap() / dt) *
                               Np * w;
 
             local_Jac