diff --git a/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h
index d35c32b83640ae384004c3e5c4ff667e7fe36253..8e09f59a2460bc396ae0ec7e34cce233f84eac64 100644
--- a/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h
+++ b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h
@@ -54,6 +54,47 @@ void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
     }
 }
 
+template <typename ShapeFunction, typename IntegrationMethod,
+          unsigned GlobalDim>
+Eigen::Vector3d
+LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::getFlux(
+    MathLib::Point3d const& p_local_coords, double const t,
+    std::vector<double> const& local_x) const
+{
+    // eval dNdx and invJ at p
+    auto const fe =
+        NumLib::createIsoparametricFiniteElement<ShapeFunction,
+                                                 ShapeMatricesType>(_element);
+
+    typename ShapeMatricesType::ShapeMatrices shape_matrices(
+        ShapeFunction::DIM, GlobalDim, ShapeFunction::NPOINTS);
+
+    // Note: Axial symmetry is set to false here, because we only need dNdx
+    // here, which is not affected by axial symmetry.
+    fe.computeShapeFunctions(p_local_coords.getCoords(), shape_matrices,
+                             GlobalDim, false);
+
+    // create pos object to access the correct media property
+    ParameterLib::SpatialPosition pos;
+    pos.setElementID(_element.getID());
+    const int material_id = _material_properties.getMaterialID(pos);
+
+    double pressure = 0.0;
+    NumLib::shapeFunctionInterpolate(local_x, shape_matrices.N, pressure);
+    const Eigen::MatrixXd& permeability = _material_properties.getPermeability(
+        material_id, t, pos, _element.getDimension(), pressure,
+        _reference_temperature);
+    const double mu =
+        _material_properties.getViscosity(pressure, _reference_temperature);
+
+    Eigen::Vector3d flux;
+    flux.head<GlobalDim>() =
+        -permeability / mu * shape_matrices.dNdx *
+        Eigen::Map<const NodalVectorType>(local_x.data(), local_x.size());
+
+    return flux;
+}
+
 template <typename ShapeFunction, typename IntegrationMethod,
           unsigned GlobalDim>
 template <typename LaplacianGravityVelocityCalculator>
diff --git a/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h
index 4cbd7559c7509eafa6779b3df2875bf3cebbaff1..7295ebc69905e4733d3594bc8673f39c226a8134 100644
--- a/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h
+++ b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h
@@ -128,6 +128,12 @@ public:
                   std::vector<double>& local_K_data,
                   std::vector<double>& local_b_data) override;
 
+    /// Computes the flux in the point \c p_local_coords that is given in local
+    /// coordinates using the values from \c local_x.
+    Eigen::Vector3d getFlux(MathLib::Point3d const& p_local_coords,
+                            double const t,
+                            std::vector<double> const& local_x) const override;
+
     Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
         const unsigned integration_point) const override
     {
diff --git a/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp b/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp
index 2f83c129e76744946d6d0d1ed4c152d42593d4ae..cf1ef6dc527dce85f5c9478142b0c6dbbfdf22bb 100644
--- a/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp
+++ b/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp
@@ -125,5 +125,21 @@ void LiquidFlowProcess::computeSecondaryVariableConcrete(const double t,
         x, _coupled_solutions);
 }
 
+Eigen::Vector3d LiquidFlowProcess::getFlux(
+    std::size_t const element_id,
+    MathLib::Point3d const& p,
+    double const t,
+    std::vector<GlobalVector*> const& x) const
+{
+    // fetch local_x from primary variable
+    std::vector<GlobalIndexType> indices_cache;
+    auto const r_c_indices = NumLib::getRowColumnIndices(
+        element_id, *_local_to_global_index_map, indices_cache);
+    constexpr int process_id = 0;  // monolithic scheme.
+    std::vector<double> local_x(x[process_id]->get(r_c_indices.rows));
+
+    return _local_assemblers[element_id]->getFlux(p, t, local_x);
+}
+
 }  // namespace LiquidFlow
 }  // namespace ProcessLib
diff --git a/ProcessLib/LiquidFlow/LiquidFlowProcess.h b/ProcessLib/LiquidFlow/LiquidFlowProcess.h
index 5501699f56c924200efc08c182e970c7881d4765..eea5c460f90a107eb28abe6c5d95a1812900f479 100644
--- a/ProcessLib/LiquidFlow/LiquidFlowProcess.h
+++ b/ProcessLib/LiquidFlow/LiquidFlowProcess.h
@@ -78,6 +78,12 @@ public:
                                           int const process_id) override;
 
     bool isLinear() const override { return true; }
+
+    Eigen::Vector3d getFlux(std::size_t const element_id,
+                            MathLib::Point3d const& p,
+                            double const t,
+                            std::vector<GlobalVector*> const& x) const override;
+
     int getGravitationalAxisID() const { return _gravitational_axis_id; }
     double getGravitationalAcceleration() const
     {