diff --git a/ProcessLib/HeatConduction/HeatConductionProcess.h b/ProcessLib/HeatConduction/HeatConductionProcess.h
index 7b754a2d36f94812a91398f66a076838f4ba992c..5d180db20656d9b4e3796033dc1ef44002f142d6 100644
--- a/ProcessLib/HeatConduction/HeatConductionProcess.h
+++ b/ProcessLib/HeatConduction/HeatConductionProcess.h
@@ -38,7 +38,9 @@ public:
     //! @{
 
     bool isLinear() const override { return true; }
-    //! @}
+
+    void computeSecondaryVariableConcrete(double const t,
+                                          GlobalVector const& x) override;
 
 private:
     void initializeConcreteProcess(
diff --git a/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h
index 6886303a067ec6116c31a67d5e6647e7f5c5d712..de4116f0bab853e0bec6c81347f2e6d1ba9f5985 100644
--- a/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h
+++ b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h
@@ -154,7 +154,8 @@ void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
 template <typename ShapeFunction, typename IntegrationMethod,
           unsigned GlobalDim>
 void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
-    computeSecondaryVariableConcrete(std::vector<double> const& local_x)
+    computeSecondaryVariableConcrete(double const t,
+                                     std::vector<double> const& local_x)
 {
     auto const local_matrix_size = local_x.size();
     assert(local_matrix_size == ShapeFunction::NPOINTS);
@@ -194,7 +195,7 @@ void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
             GlobalDimVectorType darcy_velocity = -K * sm.dNdx * local_p_vec;
             // gravity term
             if (_gravitational_axis_id >= 0)
-                darcy_velocity[GlobalDim - 1] -= K * rho_g;
+                darcy_velocity[_gravitational_axis_id] -= K * rho_g;
             for (unsigned d = 0; d < GlobalDim; ++d)
             {
                 _darcy_velocities[d][ip] = darcy_velocity[d];
@@ -203,10 +204,12 @@ void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
         else
         {
             // Compute the velocity
-            GlobalDimVectorType darcy_velocity = -perm * sm.dNdx * local_p_vec / mu;
+            GlobalDimVectorType darcy_velocity
+                    = -perm * sm.dNdx * local_p_vec / mu;
             if (_gravitational_axis_id >= 0)
             {
-                darcy_velocity.noalias() -=  rho_g * perm.col(GlobalDim - 1) / mu;
+                darcy_velocity.noalias() -=
+                        rho_g * perm.col(_gravitational_axis_id) / mu;
             }
             for (unsigned d = 0; d < GlobalDim; ++d)
             {
diff --git a/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h
index b75a21fecd3f7a23d42cc2c9199650ecc6923bb0..ac6264efc881414a8b17ff609a2bf32cdafd6421 100644
--- a/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h
+++ b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h
@@ -85,7 +85,8 @@ public:
                   std::vector<double>& local_K_data,
                   std::vector<double>& local_b_data) override;
 
-    void computeSecondaryVariableConcrete(std::vector<double> const& local_x) override;
+    void computeSecondaryVariableConcrete(double const /*t*/,
+                                std::vector<double> const& local_x) 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 ef4ca9047bc305a057f8c88197730f8fa4dbba07..07dbb24c6c80df4251461144e62bfe7e78b6e1cc 100644
--- a/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp
+++ b/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp
@@ -113,12 +113,13 @@ void LiquidFlowProcess::assembleWithJacobianConcreteProcess(
         dx_dx, M, K, b, Jac);
 }
 
-void LiquidFlowProcess::computeSecondaryVariableConcrete(GlobalVector const& x)
+void LiquidFlowProcess::computeSecondaryVariableConcrete(const double t,
+                                                         GlobalVector const& x)
 {
     DBUG("Compute the velocity for LiquidFlowProcess.");
     GlobalExecutor::executeMemberOnDereferenced(
             &LiquidFlowLocalAssemblerInterface::computeSecondaryVariable,
-            _local_assemblers, *_local_to_global_index_map, x);
+            _local_assemblers, *_local_to_global_index_map, t, x);
 }
         
 }  // end of namespace
diff --git a/ProcessLib/LiquidFlow/LiquidFlowProcess.h b/ProcessLib/LiquidFlow/LiquidFlowProcess.h
index 550d6038f16c1020541b6653236d6be799b4c993..60f0151e3f925b9a32ddacc55da4b2ac23d60421 100644
--- a/ProcessLib/LiquidFlow/LiquidFlowProcess.h
+++ b/ProcessLib/LiquidFlow/LiquidFlowProcess.h
@@ -69,7 +69,8 @@ public:
         double const gravitational_acceleration,
         BaseLib::ConfigTree const& config);
 
-    void computeSecondaryVariableConcrete(GlobalVector const& x) override;
+    void computeSecondaryVariableConcrete(double const t,
+                                          GlobalVector const& x) override;
 
     bool isLinear() const override { return true; }
 private:
diff --git a/ProcessLib/LocalAssemblerInterface.cpp b/ProcessLib/LocalAssemblerInterface.cpp
index e756ae1d92e10f44670deefa148e865bb6d9924b..c35897c1f4a981ebe1348ede53bdf09486247c69 100644
--- a/ProcessLib/LocalAssemblerInterface.cpp
+++ b/ProcessLib/LocalAssemblerInterface.cpp
@@ -26,13 +26,14 @@ void LocalAssemblerInterface::assembleWithJacobian(
         "assembler.");
 }
 
-void LocalAssemblerInterface::computeSecondaryVariable(std::size_t const mesh_item_id,
+void LocalAssemblerInterface::computeSecondaryVariable(
+                              std::size_t const mesh_item_id,
                               NumLib::LocalToGlobalIndexMap const& dof_table,
-                              GlobalVector const& x)
+                              double const t, GlobalVector const& x)
 {
     auto const indices = NumLib::getIndices(mesh_item_id, dof_table);
     auto const local_x = x.get(indices);
-    computeSecondaryVariableConcrete(local_x);
+    computeSecondaryVariableConcrete(t, local_x);
 }
 
 void LocalAssemblerInterface::preTimestep(
diff --git a/ProcessLib/LocalAssemblerInterface.h b/ProcessLib/LocalAssemblerInterface.h
index 372b6101fc273b288098e7c533a52589daf657a0..959ab56d9f5ea91fefb62d7660d6797032c2f324 100644
--- a/ProcessLib/LocalAssemblerInterface.h
+++ b/ProcessLib/LocalAssemblerInterface.h
@@ -47,7 +47,7 @@ public:
 
     virtual void computeSecondaryVariable(std::size_t const mesh_item_id,
                               NumLib::LocalToGlobalIndexMap const& dof_table,
-                              GlobalVector const& x);
+                              const double t, GlobalVector const& x);
 
     virtual void preTimestep(std::size_t const mesh_item_id,
                              NumLib::LocalToGlobalIndexMap const& dof_table,
@@ -76,7 +76,7 @@ private:
     virtual void postTimestepConcrete(std::vector<double> const& /*local_x*/) {}
 
     virtual void computeSecondaryVariableConcrete
-                         (std::vector<double> const& /*local_x*/) {}
+                (double const /*t*/, std::vector<double> const& /*local_x*/) {}
 };
 
 } // namespace ProcessLib
diff --git a/ProcessLib/Process.cpp b/ProcessLib/Process.cpp
index 620ae625af48a2b2bd3838c459aa890cd8689b00..0e906ca47769d21f685cbf24121f2bf4a167544b 100644
--- a/ProcessLib/Process.cpp
+++ b/ProcessLib/Process.cpp
@@ -242,9 +242,9 @@ void Process::postTimestep(GlobalVector const& x)
     postTimestepConcreteProcess(x);
 }
 
-void Process::computeSecondaryVariable(GlobalVector const& x)
+void Process::computeSecondaryVariable( const double t, GlobalVector const& x)
 {
-    computeSecondaryVariableConcrete(x);
+    computeSecondaryVariableConcrete(t, x);
 }
 
 void Process::preIteration(const unsigned iter, const GlobalVector &x)
diff --git a/ProcessLib/Process.h b/ProcessLib/Process.h
index 42d6b50a5098df4a3d31f6bf4cb36eade3180500..dda1c820decd84fcf28722c874fb8e80ae121ad0 100644
--- a/ProcessLib/Process.h
+++ b/ProcessLib/Process.h
@@ -59,8 +59,8 @@ public:
     void preIteration(const unsigned iter,
                       GlobalVector const& x) override final;
 
-    /// computeSecondaryVariable for the coupled equations or for output.
-    void computeSecondaryVariable(GlobalVector const& x);
+    /// compute secondary variables for the coupled equations or for output.
+    void computeSecondaryVariable(const double t, GlobalVector const& x);
 
     NumLib::IterationResult postIteration(GlobalVector const& x) override final;
 
@@ -151,7 +151,8 @@ private:
     virtual void preIterationConcreteProcess(const unsigned /*iter*/,
                                              GlobalVector const& /*x*/){}
 
-    virtual void computeSecondaryVariableConcrete(GlobalVector const& /*x*/) {};
+    virtual void computeSecondaryVariableConcrete(const double /*t*/,
+                                                  GlobalVector const& /*x*/) {};
 
     virtual NumLib::IterationResult postIterationConcreteProcess(
         GlobalVector const& /*x*/)
diff --git a/ProcessLib/UncoupledProcessesTimeLoop.cpp b/ProcessLib/UncoupledProcessesTimeLoop.cpp
index fce35e90c7b4959298206308ff3627dfe813359c..8ba27a8827246727f5c1441983b9d7b74ff4de43 100644
--- a/ProcessLib/UncoupledProcessesTimeLoop.cpp
+++ b/ProcessLib/UncoupledProcessesTimeLoop.cpp
@@ -505,7 +505,7 @@ bool UncoupledProcessesTimeLoop::loop()
             nonlinear_solver_succeeded = solveOneTimeStepOneProcess(
                 x, timestep, t, delta_t, *spd, *_output);
 
-            pcs.computeSecondaryVariable(x);
+            pcs.computeSecondaryVariable(t, x);
 
             INFO("[time] Solving process #%u took %g s in timestep #%u.",
                  timestep, time_timestep.elapsed(), pcs_idx);