diff --git a/ProcessLib/HT/HTFEM.h b/ProcessLib/HT/HTFEM.h
index 8d9a61ae3b5d2182e2efe4ef5ee3da14ed8bf891..72f0d5a7580c246de592077a70ff46d9b70086d3 100644
--- a/ProcessLib/HT/HTFEM.h
+++ b/ProcessLib/HT/HTFEM.h
@@ -261,8 +261,54 @@ public:
             Bp += w * density_water_T * dNdx.transpose() * perm_over_visc * b;
             /* with Oberbeck-Boussing assumption density difference only exists
              * in buoyancy effects */
+        }
+    }
+
+    void computeSecondaryVariableConcrete(
+        double const t, std::vector<double> const& local_x) override
+    {
+        SpatialPosition pos;
+        pos.setElementID(_element.getID());
+
+        auto const K =
+            _process_data.porous_media_properties.getIntrinsicPermeability(t,
+                                                                           pos);
+        MaterialLib::Fluid::FluidProperty::ArrayType vars;
+
+        unsigned const n_integration_points =
+            _integration_method.getNumberOfPoints();
+
+        auto const p_nodal_values = Eigen::Map<const NodalVectorType>(
+            &local_x[ShapeFunction::NPOINTS], ShapeFunction::NPOINTS);
+
+        for (unsigned ip = 0; ip < n_integration_points; ++ip)
+        {
+            auto const& ip_data = _ip_data[ip];
+            auto const& N = ip_data.N;
+            auto const& dNdx = ip_data.dNdx;
+
+            double T_int_pt = 0.0;
+            double p_int_pt = 0.0;
+            NumLib::shapeFunctionInterpolate(local_x, N, T_int_pt, p_int_pt);
+            vars[static_cast<int>(
+                MaterialLib::Fluid::PropertyVariableType::T)] = T_int_pt;
+            vars[static_cast<int>(
+                MaterialLib::Fluid::PropertyVariableType::p)] = p_int_pt;
+
+            auto const mu = _process_data.fluid_properties->getValue(
+                MaterialLib::Fluid::FluidPropertyType::Viscosity, vars);
+            GlobalDimMatrixType const K_over_mu = K / mu;
+
+            GlobalDimVectorType velocity = -K_over_mu * dNdx * p_nodal_values;
+            if (_process_data.has_gravity)
+            {
+                auto const rho_w = _process_data.fluid_properties->getValue(
+                    MaterialLib::Fluid::FluidPropertyType::Density, vars);
+                auto const b = _process_data.specific_body_force;
+                // here it is assumed that the vector b is directed 'downwards'
+                velocity += K_over_mu * rho_w * b;
+            }
 
-            // velocity computed for output.
             for (unsigned d = 0; d < GlobalDim; ++d)
             {
                 _darcy_velocities[d][ip] = velocity[d];
diff --git a/ProcessLib/HT/HTProcess.cpp b/ProcessLib/HT/HTProcess.cpp
index 3cebe46de0224d92e778695a0e246a68b3e030fc..7a70c5ad7671133cc2b60e663432a0bba279e2f8 100644
--- a/ProcessLib/HT/HTProcess.cpp
+++ b/ProcessLib/HT/HTProcess.cpp
@@ -98,6 +98,16 @@ void HTProcess::assembleWithJacobianConcreteProcess(
         dx_dx, M, K, b, Jac, coupling_term);
 }
 
+void HTProcess::computeSecondaryVariableConcrete(
+    double const t, GlobalVector const& x,
+    StaggeredCouplingTerm const& coupling_term)
+{
+    DBUG("Compute the Darcy velocity for HTProcess.");
+    GlobalExecutor::executeMemberOnDereferenced(
+        &HTLocalAssemblerInterface::computeSecondaryVariable, _local_assemblers,
+        *_local_to_global_index_map, t, x, coupling_term);
+}
+
 }  // namespace HT
 }  // namespace ProcessLib
 
diff --git a/ProcessLib/HT/HTProcess.h b/ProcessLib/HT/HTProcess.h
index a5048e295cd1daaa2f4507d3d955a06bb0a8d8d6..6dcdaaeaf2669de46918240013aba189f7687cdb 100644
--- a/ProcessLib/HT/HTProcess.h
+++ b/ProcessLib/HT/HTProcess.h
@@ -56,6 +56,9 @@ public:
 
     bool isLinear() const override { return false; }
     //! @}
+    void computeSecondaryVariableConcrete(
+        double const t, GlobalVector const& x,
+        StaggeredCouplingTerm const& coupling_term) override;
 
 private:
     void initializeConcreteProcess(