diff --git a/ProcessLib/TH2M/ConstitutiveRelations/Base.h b/ProcessLib/TH2M/ConstitutiveRelations/Base.h
index 8739cfbab3877861baa34edea4de1b6c28db9878..e5ec55f00b49947babc5bbb64c95ea1a5de9dd3b 100644
--- a/ProcessLib/TH2M/ConstitutiveRelations/Base.h
+++ b/ProcessLib/TH2M/ConstitutiveRelations/Base.h
@@ -59,6 +59,14 @@ using ReferenceTemperatureData =
 using GasPressureData = BaseLib::StrongType<double, struct GasPressureTag>;
 using CapillaryPressureData =
     BaseLib::StrongType<double, struct CapillaryPressureTag>;
+template <int DisplacementDim>
+using GasPressureGradientData =
+    BaseLib::StrongType<GlobalDimVector<DisplacementDim>,
+                        struct GasPressureGradientTag>;
+template <int DisplacementDim>
+using CapillaryPressureGradientData =
+    BaseLib::StrongType<GlobalDimVector<DisplacementDim>,
+                        struct CapillaryPressureGradientTag>;
 
 template <int DisplacementDim>
 using SpecificBodyForceData =
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveModels.h b/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveModels.h
index df7bee1cd8a9b7fbd249839e58bd4247a35870e8..3ac2a4840765c8fa7bf730eb159cc3c0a068ed78 100644
--- a/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveModels.h
+++ b/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveModels.h
@@ -85,6 +85,7 @@ struct ConstitutiveModels
     InternalEnergyModel internal_energy_model;
     EffectiveVolumetricEnthalpyModel effective_volumetric_enthalpy_model;
     GravityModel<DisplacementDim> gravity_model;
+    DarcyVelocityModel<DisplacementDim> darcy_velocity_model;
     FC1Model<DisplacementDim> fC_1_model;
     FC2aModel fC_2a_model;
     FC3aModel fC_3a_model;
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/DarcyVelocity.cpp b/ProcessLib/TH2M/ConstitutiveRelations/DarcyVelocity.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..5083d8c65c25e801bdcddc94ba51c3c2da8dcb59
--- /dev/null
+++ b/ProcessLib/TH2M/ConstitutiveRelations/DarcyVelocity.cpp
@@ -0,0 +1,43 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2024, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ */
+
+#include "DarcyVelocity.h"
+
+namespace ProcessLib::TH2M
+{
+namespace ConstitutiveRelations
+{
+template <int DisplacementDim>
+void DarcyVelocityModel<DisplacementDim>::eval(
+    CapillaryPressureGradientData<DisplacementDim> const& grad_p_cap,
+    FluidDensityData const& fluid_density_data,
+    GasPressureGradientData<DisplacementDim> const& grad_p_GR,
+    PermeabilityData<DisplacementDim> const& permeability_data,
+    SpecificBodyForceData<DisplacementDim> const& specific_body_force,
+    ViscosityData const& viscosity_data,
+    DarcyVelocityData<DisplacementDim>& darcy_velocity_data) const
+{
+    auto const k_over_mu_G =
+        permeability_data.Ki * permeability_data.k_rel_G / viscosity_data.mu_GR;
+    auto const k_over_mu_L =
+        permeability_data.Ki * permeability_data.k_rel_L / viscosity_data.mu_LR;
+
+    darcy_velocity_data.w_GS =
+        k_over_mu_G * fluid_density_data.rho_GR * specific_body_force() -
+        k_over_mu_G * grad_p_GR();
+    darcy_velocity_data.w_LS =
+        k_over_mu_L * grad_p_cap() +
+        k_over_mu_L * fluid_density_data.rho_LR * specific_body_force() -
+        k_over_mu_L * grad_p_GR();
+}
+
+template struct DarcyVelocityModel<2>;
+template struct DarcyVelocityModel<3>;
+}  // namespace ConstitutiveRelations
+}  // namespace ProcessLib::TH2M
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/DarcyVelocity.h b/ProcessLib/TH2M/ConstitutiveRelations/DarcyVelocity.h
index 84e66d31d02bcf744f6d5affa6a7183695d2392b..91751b830f485a309fea551dd463be2e9d7769d9 100644
--- a/ProcessLib/TH2M/ConstitutiveRelations/DarcyVelocity.h
+++ b/ProcessLib/TH2M/ConstitutiveRelations/DarcyVelocity.h
@@ -10,7 +10,10 @@
 #pragma once
 
 #include "Base.h"
+#include "FluidDensity.h"
+#include "PermeabilityData.h"
 #include "ProcessLib/Reflection/ReflectionData.h"
+#include "Viscosity.h"
 
 namespace ProcessLib::TH2M
 {
@@ -32,5 +35,20 @@ struct DarcyVelocityData
             R::makeReflectionData("velocity_liquid", &Self::w_LS)};
     }
 };
+
+template <int DisplacementDim>
+struct DarcyVelocityModel
+{
+    void eval(CapillaryPressureGradientData<DisplacementDim> const& grad_p_cap,
+              FluidDensityData const& fluid_density_data,
+              GasPressureGradientData<DisplacementDim> const& grad_p_GR,
+              PermeabilityData<DisplacementDim> const& permeability_data,
+              SpecificBodyForceData<DisplacementDim> const& specific_body_force,
+              ViscosityData const& viscosity_data,
+              DarcyVelocityData<DisplacementDim>& darcy_velocity_data) const;
+};
+
+extern template struct DarcyVelocityModel<2>;
+extern template struct DarcyVelocityModel<3>;
 }  // namespace ConstitutiveRelations
 }  // namespace ProcessLib::TH2M
diff --git a/ProcessLib/TH2M/TH2MFEM-impl.h b/ProcessLib/TH2M/TH2MFEM-impl.h
index cc118746435e2c1f2c8451d2828908357c2d6e84..c87165fd500875929b65544adb9737859ba26bd7 100644
--- a/ProcessLib/TH2M/TH2MFEM-impl.h
+++ b/ProcessLib/TH2M/TH2MFEM-impl.h
@@ -162,6 +162,10 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
         GlobalDimVectorType const gradpGR = gradNp * gas_pressure;
         GlobalDimVectorType const gradpCap = gradNp * capillary_pressure;
         GlobalDimVectorType const gradT = gradNp * temperature;
+        ConstitutiveRelations::GasPressureGradientData<DisplacementDim> const
+            grad_p_GR{gradpGR};
+        ConstitutiveRelations::CapillaryPressureGradientData<
+            DisplacementDim> const grad_p_cap{gradpCap};
 
         // medium properties
         models.elastic_tangent_stiffness_model.eval({pos, t, dt}, T_data,
@@ -487,21 +491,15 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
         // Derivatives for Jacobian
         // ---------------------------------------------------------------------
 
-        auto const& b = this->process_data_.specific_body_force;
-        GlobalDimMatrixType const k_over_mu_G =
-            ip_out.permeability_data.Ki * ip_out.permeability_data.k_rel_G /
-            ip_cv.viscosity_data.mu_GR;
-        GlobalDimMatrixType const k_over_mu_L =
-            ip_out.permeability_data.Ki * ip_out.permeability_data.k_rel_L /
-            ip_cv.viscosity_data.mu_LR;
-
-        ip_out.darcy_velocity_data.w_GS =
-            k_over_mu_G * ip_out.fluid_density_data.rho_GR * b -
-            k_over_mu_G * gradpGR;
-        ip_out.darcy_velocity_data.w_LS =
-            k_over_mu_L * gradpCap +
-            k_over_mu_L * ip_out.fluid_density_data.rho_LR * b -
-            k_over_mu_L * gradpGR;
+        models.darcy_velocity_model.eval(
+            grad_p_cap,
+            ip_out.fluid_density_data,
+            grad_p_GR,
+            ip_out.permeability_data,
+            ConstitutiveRelations::SpecificBodyForceData<DisplacementDim>{
+                this->process_data_.specific_body_force},
+            ip_cv.viscosity_data,
+            ip_out.darcy_velocity_data);
 
         models.fT_2_model.eval(ip_out.darcy_velocity_data,
                                ip_out.enthalpy_data,