diff --git a/ProcessLib/RichardsMechanics/ConstitutiveRelations/ConstitutiveData.h b/ProcessLib/RichardsMechanics/ConstitutiveRelations/ConstitutiveData.h
index 242707fb1efebf9ba3a71fd91256facf5d5f2007..417faf5bc92ff737397136be278733d0a3334686 100644
--- a/ProcessLib/RichardsMechanics/ConstitutiveRelations/ConstitutiveData.h
+++ b/ProcessLib/RichardsMechanics/ConstitutiveRelations/ConstitutiveData.h
@@ -16,6 +16,7 @@
 #include "MicroSaturation.h"
 #include "ProcessLib/ConstitutiveRelations/Base.h"
 #include "ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/Biot.h"
+#include "ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/DarcyLaw.h"
 #include "ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/LiquidViscosity.h"
 #include "ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/PermeabilityData.h"
 #include "ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/Porosity.h"
@@ -49,7 +50,8 @@ using StatefulDataPrev = ProcessLib::ConstitutiveRelations::PrevStateOf<
 
 /// Data that is needed for output purposes, but not directly for the assembly.
 template <int DisplacementDim>
-using OutputData = std::tuple<>;
+using OutputData = std::tuple<
+    ProcessLib::ThermoRichardsMechanics::DarcyLawData<DisplacementDim>>;
 
 /// Data that is needed for the equation system assembly.
 template <int DisplacementDim>
diff --git a/ProcessLib/RichardsMechanics/IntegrationPointData.h b/ProcessLib/RichardsMechanics/IntegrationPointData.h
index 5c3908bdd707f09f1053ef8f8fa87e86f89c5e1c..ee813f5b7ba1697fa0a89b56dd15ec707cc632d3 100644
--- a/ProcessLib/RichardsMechanics/IntegrationPointData.h
+++ b/ProcessLib/RichardsMechanics/IntegrationPointData.h
@@ -40,8 +40,6 @@ struct IntegrationPointData final
     typename ShapeMatricesTypePressure::NodalRowVectorType N_p;
     typename ShapeMatricesTypePressure::GlobalDimNodalMatrixType dNdx_p;
 
-    typename ShapeMatricesTypePressure::GlobalDimVectorType v_darcy;
-
     double dry_density_solid = std::numeric_limits<double>::quiet_NaN();
     double dry_density_pellet_saturated =
         std::numeric_limits<double>::quiet_NaN();
diff --git a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h
index c76007b6d9b951b51e78a87b74586c20f8a6c727..99c9e940d94a9c4fe5752c830a7d107267ed1a4e 100644
--- a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h
+++ b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h
@@ -1718,7 +1718,9 @@ std::vector<double> const& RichardsMechanicsLocalAssembler<
 
     for (unsigned ip = 0; ip < n_integration_points; ip++)
     {
-        cache_matrix.col(ip).noalias() = _ip_data[ip].v_darcy;
+        cache_matrix.col(ip).noalias() = *std::get<
+            ProcessLib::ThermoRichardsMechanics::DarcyLawData<DisplacementDim>>(
+            this->output_data_[ip]);
     }
 
     return cache;
@@ -2247,8 +2249,10 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
 
         // Compute the velocity
         auto const& dNdx_p = _ip_data[ip].dNdx_p;
-        _ip_data[ip].v_darcy.noalias() =
-            -K_over_mu * dNdx_p * p_L + rho_LR * K_over_mu * b;
+        std::get<
+            ProcessLib::ThermoRichardsMechanics::DarcyLawData<DisplacementDim>>(
+            this->output_data_[ip])
+            ->noalias() = -K_over_mu * dNdx_p * p_L + rho_LR * K_over_mu * b;
 
         saturation_avg += S_L;
         porosity_avg += phi;