diff --git a/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp b/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp
index beacb9f51b7e3f5f22b9df9af282af69ab1682ab..d398f7fc5cbdf18762a377b9c935605a7f41febf 100644
--- a/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp
+++ b/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp
@@ -66,18 +66,20 @@ void LiquidFlowProcess::initializeConcreteProcess(
         pv.getShapeFunctionOrder(), _local_assemblers,
         mesh.isAxiallySymmetric(), integration_order, _gravitational_axis_id,
         _gravitational_acceleration, _reference_temperature,
-        *_material_properties);
+        *_material_properties, this->_coupling_term);
 
     _secondary_variables.addSecondaryVariable(
         "darcy_velocity",
-        makeExtrapolator(mesh.getDimension(),
-            getExtrapolator(), _local_assemblers,
+        makeExtrapolator(
+            mesh.getDimension(), getExtrapolator(), _local_assemblers,
             &LiquidFlowLocalAssemblerInterface::getIntPtDarcyVelocity));
 }
 
-void LiquidFlowProcess::assembleConcreteProcess(
-    const double t, GlobalVector const& x, GlobalMatrix& M, GlobalMatrix& K,
-    GlobalVector& b)
+void LiquidFlowProcess::assembleConcreteProcess(const double t,
+                                                GlobalVector const& x,
+                                                GlobalMatrix& M,
+                                                GlobalMatrix& K,
+                                                GlobalVector& b)
 {
     DBUG("Assemble LiquidFlowProcess.");
     // Call global assembler for each local assembly item.
@@ -100,15 +102,21 @@ void LiquidFlowProcess::assembleWithJacobianConcreteProcess(
         dx_dx, M, K, b, Jac, _coupling_term);
 }
 
-void LiquidFlowProcess::computeSecondaryVariableConcrete(
-    const double t,
-    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, t, x,
-        _coupling_term);
+        _local_assemblers, *_local_to_global_index_map, t, x, _coupling_term);
+}
+
+void LiquidFlowProcess::setStaggeredCouplingTermToLocalAssemblers()
+{
+    DBUG("Compute the velocity for LiquidFlowProcess.");
+    GlobalExecutor::executeMemberOnDereferenced(
+        &LiquidFlowLocalAssemblerInterface::setStaggeredCouplingTerm,
+        _local_assemblers, _coupling_term);
 }
 
 }  // end of namespace
diff --git a/ProcessLib/LiquidFlow/LiquidFlowProcess.h b/ProcessLib/LiquidFlow/LiquidFlowProcess.h
index 9bd7152f130953eca6dd3ad15d182ee9ae81e407..0ab5f4ac69be1c2b474d958753dd313863ac5574 100644
--- a/ProcessLib/LiquidFlow/LiquidFlowProcess.h
+++ b/ProcessLib/LiquidFlow/LiquidFlowProcess.h
@@ -73,8 +73,8 @@ public:
         double const reference_temperature,
         BaseLib::ConfigTree const& config);
 
-    void computeSecondaryVariableConcrete(
-        double const t, GlobalVector const& x) override;
+    void computeSecondaryVariableConcrete(double const t,
+                                          GlobalVector const& x) override;
 
     bool isLinear() const override { return true; }
     int getGravitationalAxisID() const { return _gravitational_axis_id; }
@@ -88,14 +88,16 @@ public:
         return _material_properties.get();
     }
 
+    void setStaggeredCouplingTermToLocalAssemblers() override;
+
 private:
     void initializeConcreteProcess(
         NumLib::LocalToGlobalIndexMap const& dof_table,
         MeshLib::Mesh const& mesh, unsigned const integration_order) override;
 
-    void assembleConcreteProcess(
-        const double t, GlobalVector const& x, GlobalMatrix& M, GlobalMatrix& K,
-        GlobalVector& b) override;
+    void assembleConcreteProcess(const double t, GlobalVector const& x,
+                                 GlobalMatrix& M, GlobalMatrix& K,
+                                 GlobalVector& b) override;
 
     void assembleWithJacobianConcreteProcess(
         const double t, GlobalVector const& x, GlobalVector const& xdot,
diff --git a/ProcessLib/Process.h b/ProcessLib/Process.h
index 1e0a15beae4987439d286fc89f13a95ac63d61ea..88ef33bf5ee60e203fba6c30ad7822211c3ec475 100644
--- a/ProcessLib/Process.h
+++ b/ProcessLib/Process.h
@@ -75,6 +75,7 @@ public:
         _coupling_term = coupling_term;
     }
 
+    virtual void setStaggeredCouplingTermToLocalAssemblers() {}
     void assemble(const double t, GlobalVector const& x, GlobalMatrix& M,
                   GlobalMatrix& K, GlobalVector& b) final;