diff --git a/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.cpp b/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.cpp
index 42367c1b05eb6268d97259cd46d9cb15da9bc345..93589087db50fc1df5fbfb73451f19a4262b4186 100644
--- a/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.cpp
+++ b/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.cpp
@@ -76,7 +76,9 @@ void GroundwaterFlowProcess::assembleConcreteProcess(const double t,
                                                      GlobalVector const& x,
                                                      GlobalMatrix& M,
                                                      GlobalMatrix& K,
-                                                     GlobalVector& b)
+                                                     GlobalVector& b,
+                                                     StaggeredCouplingTerm
+                                                     const& /*coupled_term*/)
 {
     DBUG("Assemble GroundwaterFlowProcess.");
 
@@ -89,7 +91,8 @@ void GroundwaterFlowProcess::assembleConcreteProcess(const double t,
 void GroundwaterFlowProcess::assembleWithJacobianConcreteProcess(
     const double t, GlobalVector const& x, GlobalVector const& xdot,
     const double dxdot_dx, const double dx_dx, GlobalMatrix& M, GlobalMatrix& K,
-    GlobalVector& b, GlobalMatrix& Jac)
+    GlobalVector& b, GlobalMatrix& Jac,
+    StaggeredCouplingTerm const& coupled_term)
 {
     DBUG("AssembleWithJacobian GroundwaterFlowProcess.");
 
@@ -97,7 +100,7 @@ void GroundwaterFlowProcess::assembleWithJacobianConcreteProcess(
     GlobalExecutor::executeMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
         _local_assemblers, *_local_to_global_index_map, t, x, xdot, dxdot_dx,
-        dx_dx, M, K, b, Jac);
+        dx_dx, M, K, b, Jac, coupled_term);
 }
 
 
diff --git a/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.h b/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.h
index eaebcab80d6eb21a93d83ce274ed96255b8645a8..69bb1a685a9c569cacc5c0c5a7efccd326f8af09 100644
--- a/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.h
+++ b/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.h
@@ -101,12 +101,15 @@ private:
 
     void assembleConcreteProcess(const double t, GlobalVector const& x,
                                  GlobalMatrix& M, GlobalMatrix& K,
-                                 GlobalVector& b) override;
+                                 GlobalVector& b,
+                                 StaggeredCouplingTerm const& coupled_term
+                                ) override;
 
     void assembleWithJacobianConcreteProcess(
         const double t, GlobalVector const& x, GlobalVector const& xdot,
         const double dxdot_dx, const double dx_dx, GlobalMatrix& M,
-        GlobalMatrix& K, GlobalVector& b, GlobalMatrix& Jac) override;
+        GlobalMatrix& K, GlobalVector& b, GlobalMatrix& Jac,
+        StaggeredCouplingTerm const& coupled_term) override;
 
     GroundwaterFlowProcessData _process_data;
 
diff --git a/ProcessLib/HT/HTProcess.cpp b/ProcessLib/HT/HTProcess.cpp
index 2ce02f4b28c1ee5b73199b31c2afa86053e5ce08..d8f4eddb63aa472d42d75701784dc1d0ac6c6aee 100644
--- a/ProcessLib/HT/HTProcess.cpp
+++ b/ProcessLib/HT/HTProcess.cpp
@@ -71,7 +71,9 @@ void HTProcess::assembleConcreteProcess(const double t,
                                         GlobalVector const& x,
                                         GlobalMatrix& M,
                                         GlobalMatrix& K,
-                                        GlobalVector& b)
+                                        GlobalVector& b,
+                                        StaggeredCouplingTerm const&
+                                        /*coupled_term*/)
 {
     DBUG("Assemble HTProcess.");
 
@@ -84,7 +86,8 @@ void HTProcess::assembleConcreteProcess(const double t,
 void HTProcess::assembleWithJacobianConcreteProcess(
     const double t, GlobalVector const& x, GlobalVector const& xdot,
     const double dxdot_dx, const double dx_dx, GlobalMatrix& M, GlobalMatrix& K,
-    GlobalVector& b, GlobalMatrix& Jac)
+    GlobalVector& b, GlobalMatrix& Jac,
+    StaggeredCouplingTerm const& coupled_term)
 {
     DBUG("AssembleWithJacobian HTProcess.");
 
@@ -92,7 +95,7 @@ void HTProcess::assembleWithJacobianConcreteProcess(
     GlobalExecutor::executeMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
         _local_assemblers, *_local_to_global_index_map, t, x, xdot, dxdot_dx,
-        dx_dx, M, K, b, Jac);
+        dx_dx, M, K, b, Jac, coupled_term);
 }
 
 }  // namespace HT
diff --git a/ProcessLib/HT/HTProcess.h b/ProcessLib/HT/HTProcess.h
index 63a9f058dd002250c06f504c48bd850e456f0143..1da30247709dc883431382913ed455fbbb1a952a 100644
--- a/ProcessLib/HT/HTProcess.h
+++ b/ProcessLib/HT/HTProcess.h
@@ -68,12 +68,15 @@ private:
 
     void assembleConcreteProcess(const double t, GlobalVector const& x,
                                  GlobalMatrix& M, GlobalMatrix& K,
-                                 GlobalVector& b) override;
+                                 GlobalVector& b,
+                                 StaggeredCouplingTerm const& coupled_term
+                                ) override;
 
     void assembleWithJacobianConcreteProcess(
         const double t, GlobalVector const& x, GlobalVector const& xdot,
         const double dxdot_dx, const double dx_dx, GlobalMatrix& M,
-        GlobalMatrix& K, GlobalVector& b, GlobalMatrix& Jac) override;
+        GlobalMatrix& K, GlobalVector& b, GlobalMatrix& Jac,
+        StaggeredCouplingTerm const& coupled_term) override;
 
     HTProcessData _process_data;
 
diff --git a/ProcessLib/HeatConduction/HeatConductionProcess.cpp b/ProcessLib/HeatConduction/HeatConductionProcess.cpp
index b82273f8b6f06d6c31e6048d00bba724320e7119..3ceed6b6453b7e9534bab6a6f26877541c8511ed 100644
--- a/ProcessLib/HeatConduction/HeatConductionProcess.cpp
+++ b/ProcessLib/HeatConduction/HeatConductionProcess.cpp
@@ -72,7 +72,9 @@ void HeatConductionProcess::assembleConcreteProcess(const double t,
                                                     GlobalVector const& x,
                                                     GlobalMatrix& M,
                                                     GlobalMatrix& K,
-                                                    GlobalVector& b)
+                                                    GlobalVector& b,
+                                                    StaggeredCouplingTerm
+                                                    const& /*coupled_term*/)
 {
     DBUG("Assemble HeatConductionProcess.");
 
@@ -85,7 +87,8 @@ void HeatConductionProcess::assembleConcreteProcess(const double t,
 void HeatConductionProcess::assembleWithJacobianConcreteProcess(
     const double t, GlobalVector const& x, GlobalVector const& xdot,
     const double dxdot_dx, const double dx_dx, GlobalMatrix& M, GlobalMatrix& K,
-    GlobalVector& b, GlobalMatrix& Jac)
+    GlobalVector& b, GlobalMatrix& Jac,
+    StaggeredCouplingTerm const& coupled_term)
 {
     DBUG("AssembleWithJacobian HeatConductionProcess.");
 
@@ -93,7 +96,7 @@ void HeatConductionProcess::assembleWithJacobianConcreteProcess(
     GlobalExecutor::executeMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
         _local_assemblers, *_local_to_global_index_map, t, x, xdot, dxdot_dx,
-        dx_dx, M, K, b, Jac);
+        dx_dx, M, K, b, Jac, coupled_term);
 }
 
 void HeatConductionProcess::computeSecondaryVariableConcrete(const double t,
diff --git a/ProcessLib/HeatConduction/HeatConductionProcess.h b/ProcessLib/HeatConduction/HeatConductionProcess.h
index 4f59addd53de25ebc1ef84c6546e12c3418d773e..5c94f4f95b60c12a5884015b5ec8b99e37c1bb17 100644
--- a/ProcessLib/HeatConduction/HeatConductionProcess.h
+++ b/ProcessLib/HeatConduction/HeatConductionProcess.h
@@ -52,12 +52,15 @@ private:
 
     void assembleConcreteProcess(const double t, GlobalVector const& x,
                                  GlobalMatrix& M, GlobalMatrix& K,
-                                 GlobalVector& b) override;
+                                 GlobalVector& b,
+                                 StaggeredCouplingTerm const& coupled_term
+                                ) override;
 
     void assembleWithJacobianConcreteProcess(
         const double t, GlobalVector const& x, GlobalVector const& xdot,
         const double dxdot_dx, const double dx_dx, GlobalMatrix& M,
-        GlobalMatrix& K, GlobalVector& b, GlobalMatrix& Jac) override;
+        GlobalMatrix& K, GlobalVector& b, GlobalMatrix& Jac,
+        StaggeredCouplingTerm const& coupled_term) override;
 
     HeatConductionProcessData _process_data;
 
diff --git a/ProcessLib/HydroMechanics/HydroMechanicsFEM.h b/ProcessLib/HydroMechanics/HydroMechanicsFEM.h
index 97e3f0d0eabb8632963620caf89e07a81c1c4f30..72540162dfb8c12c71127679e1f74ed4493f8b1e 100644
--- a/ProcessLib/HydroMechanics/HydroMechanicsFEM.h
+++ b/ProcessLib/HydroMechanics/HydroMechanicsFEM.h
@@ -200,7 +200,9 @@ public:
                               std::vector<double>& /*local_M_data*/,
                               std::vector<double>& /*local_K_data*/,
                               std::vector<double>& local_rhs_data,
-                              std::vector<double>& local_Jac_data) override
+                              std::vector<double>& local_Jac_data,
+                              LocalCouplingTerm const& /*coupled_term*/
+                             ) override
     {
         assert(local_x.size() == pressure_size + displacement_size);
 
diff --git a/ProcessLib/HydroMechanics/HydroMechanicsProcess.h b/ProcessLib/HydroMechanics/HydroMechanicsProcess.h
index 8f0fd00cc455ed6b16013ae6e30ebadc3619b364..0a6dc8c707a9cdb46bda71ff4271b763205b8f2c 100644
--- a/ProcessLib/HydroMechanics/HydroMechanicsProcess.h
+++ b/ProcessLib/HydroMechanics/HydroMechanicsProcess.h
@@ -109,7 +109,9 @@ private:
 
     void assembleConcreteProcess(const double t, GlobalVector const& x,
                                  GlobalMatrix& M, GlobalMatrix& K,
-                                 GlobalVector& b) override
+                                 GlobalVector& b,
+                                 StaggeredCouplingTerm const&
+                                 /*coupled_term*/) override
     {
         DBUG("Assemble HydroMechanicsProcess.");
 
@@ -122,7 +124,8 @@ private:
     void assembleWithJacobianConcreteProcess(
         const double t, GlobalVector const& x, GlobalVector const& xdot,
         const double dxdot_dx, const double dx_dx, GlobalMatrix& M,
-        GlobalMatrix& K, GlobalVector& b, GlobalMatrix& Jac) override
+        GlobalMatrix& K, GlobalVector& b, GlobalMatrix& Jac,
+        StaggeredCouplingTerm const& coupled_term) override
     {
         DBUG("AssembleJacobian HydroMechanicsProcess.");
 
@@ -130,7 +133,7 @@ private:
         GlobalExecutor::executeMemberDereferenced(
             _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
             _local_assemblers, *_local_to_global_index_map, t, x, xdot,
-            dxdot_dx, dx_dx, M, K, b, Jac);
+            dxdot_dx, dx_dx, M, K, b, Jac, coupled_term);
     }
 
     void preTimestepConcreteProcess(GlobalVector const& x, double const t,
diff --git a/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcess.h b/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcess.h
index 93dd346344602c60e4971566274580f675698b04..18f260b49bb51aaef7bdc8906af34a8b155a5358 100644
--- a/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcess.h
+++ b/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcess.h
@@ -70,7 +70,9 @@ private:
 
     void assembleConcreteProcess(const double t, GlobalVector const& x,
                                  GlobalMatrix& M, GlobalMatrix& K,
-                                 GlobalVector& b) override
+                                 GlobalVector& b,
+                                 StaggeredCouplingTerm const&
+                                 /*coupled_term*/) override
     {
         DBUG("Assemble HydroMechanicsProcess.");
 
diff --git a/ProcessLib/LIE/SmallDeformation/SmallDeformationProcess.h b/ProcessLib/LIE/SmallDeformation/SmallDeformationProcess.h
index 03adda654a987e28590c7fc0d24f5c9060cf73ff..efd4da63603719d69ab8fc87db2c9345a1f34d51 100644
--- a/ProcessLib/LIE/SmallDeformation/SmallDeformationProcess.h
+++ b/ProcessLib/LIE/SmallDeformation/SmallDeformationProcess.h
@@ -68,7 +68,9 @@ private:
 
     void assembleConcreteProcess(const double t, GlobalVector const& x,
                                  GlobalMatrix& M, GlobalMatrix& K,
-                                 GlobalVector& b) override
+                                 GlobalVector& b,
+                                 StaggeredCouplingTerm const&
+                                 /*coupled_term*/) override
     {
         DBUG("Assemble SmallDeformationProcess.");
 
@@ -81,7 +83,8 @@ private:
     void assembleWithJacobianConcreteProcess(
         const double t, GlobalVector const& x, GlobalVector const& xdot,
         const double dxdot_dx, const double dx_dx, GlobalMatrix& M,
-        GlobalMatrix& K, GlobalVector& b, GlobalMatrix& Jac) override
+        GlobalMatrix& K, GlobalVector& b, GlobalMatrix& Jac,
+        StaggeredCouplingTerm const& coupled_term) override
     {
         DBUG("AssembleWithJacobian SmallDeformationProcess.");
 
@@ -89,7 +92,7 @@ private:
         GlobalExecutor::executeMemberDereferenced(
             _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
             _local_assemblers, *_local_to_global_index_map, t, x, xdot,
-            dxdot_dx, dx_dx, M, K, b, Jac);
+            dxdot_dx, dx_dx, M, K, b, Jac, coupled_term);
     }
 
     void preTimestepConcreteProcess(GlobalVector const& x, double const t,
diff --git a/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp b/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp
index 1d9d9389dbc718e5f69603708d74015db923ffef..b5f351e009a37199b06c3b3e64dd2a6ed44ab4b0 100644
--- a/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp
+++ b/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp
@@ -92,7 +92,9 @@ void LiquidFlowProcess::assembleConcreteProcess(const double t,
                                                 GlobalVector const& x,
                                                 GlobalMatrix& M,
                                                 GlobalMatrix& K,
-                                                GlobalVector& b)
+                                                GlobalVector& b,
+                                                StaggeredCouplingTerm const&
+                                                /*coupled_term*/)
 {
     DBUG("Assemble LiquidFlowProcess.");
     // Call global assembler for each local assembly item.
@@ -104,7 +106,8 @@ void LiquidFlowProcess::assembleConcreteProcess(const double t,
 void LiquidFlowProcess::assembleWithJacobianConcreteProcess(
     const double t, GlobalVector const& x, GlobalVector const& xdot,
     const double dxdot_dx, const double dx_dx, GlobalMatrix& M, GlobalMatrix& K,
-    GlobalVector& b, GlobalMatrix& Jac)
+    GlobalVector& b, GlobalMatrix& Jac,
+    StaggeredCouplingTerm const& coupled_term)
 {
     DBUG("AssembleWithJacobian LiquidFlowProcess.");
 
@@ -112,7 +115,7 @@ void LiquidFlowProcess::assembleWithJacobianConcreteProcess(
     GlobalExecutor::executeMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
         _local_assemblers, *_local_to_global_index_map, t, x, xdot, dxdot_dx,
-        dx_dx, M, K, b, Jac);
+        dx_dx, M, K, b, Jac, coupled_term);
 }
 
 void LiquidFlowProcess::computeSecondaryVariableConcrete(const double t,
diff --git a/ProcessLib/LiquidFlow/LiquidFlowProcess.h b/ProcessLib/LiquidFlow/LiquidFlowProcess.h
index ec69f532a3cc75d1d8c8030291c5e25ff14c7cfe..d7431f6258922f6b1b1e7b1da96e5c1386ef23df 100644
--- a/ProcessLib/LiquidFlow/LiquidFlowProcess.h
+++ b/ProcessLib/LiquidFlow/LiquidFlowProcess.h
@@ -80,11 +80,6 @@ public:
     ProcessType getProcessType() const override
                      {return ProcessLib::ProcessType::LiquidFlowProcess;}
 
-    MaterialLib::Fluid::FluidProperties* getFluidProperties() const
-    {
-        return _material_properties->getFluidProperties();
-    }
-
 private:
     void initializeConcreteProcess(
         NumLib::LocalToGlobalIndexMap const& dof_table,
@@ -92,12 +87,15 @@ private:
 
     void assembleConcreteProcess(const double t, GlobalVector const& x,
                                  GlobalMatrix& M, GlobalMatrix& K,
-                                 GlobalVector& b) override;
+                                 GlobalVector& b,
+                                 StaggeredCouplingTerm const& coupled_term
+                                ) override;
 
     void assembleWithJacobianConcreteProcess(
         const double t, GlobalVector const& x, GlobalVector const& xdot,
         const double dxdot_dx, const double dx_dx, GlobalMatrix& M,
-        GlobalMatrix& K, GlobalVector& b, GlobalMatrix& Jac) override;
+        GlobalMatrix& K, GlobalVector& b, GlobalMatrix& Jac,
+        StaggeredCouplingTerm const& coupled_term) override;
 
     const int _gravitational_axis_id;
     const double _gravitational_acceleration;
diff --git a/ProcessLib/Process.cpp b/ProcessLib/Process.cpp
index 4e69b98a00c44ef8d281d0935184ed3e8e43e3f7..f41e79a61af8918b4fb436dc1500557a375088d1 100644
--- a/ProcessLib/Process.cpp
+++ b/ProcessLib/Process.cpp
@@ -122,7 +122,8 @@ void Process::assemble(const double t, GlobalVector const& x, GlobalMatrix& M,
                        GlobalMatrix& K, GlobalVector& b,
                        StaggeredCouplingTerm const& coupled_term)
 {
-    assembleConcreteProcess(t, x, M, K, b);
+    assembleConcreteProcess(t, x, M, K, b, coupled_term);
+
     _boundary_conditions.applyNaturalBC(t, x, K, b);
 }
 
@@ -133,7 +134,8 @@ void Process::assembleWithJacobian(const double t, GlobalVector const& x,
                                    GlobalVector& b, GlobalMatrix& Jac,
                                    StaggeredCouplingTerm const& coupled_term)
 {
-    assembleWithJacobianConcreteProcess(t, x, xdot, dxdot_dx, dx_dx, M, K, b, Jac);
+    assembleWithJacobianConcreteProcess(t, x, xdot, dxdot_dx, dx_dx, M, K, b,
+                                        Jac, coupled_term);
 
     // TODO apply BCs to Jacobian.
     _boundary_conditions.applyNaturalBC(t, x, K, b);
diff --git a/ProcessLib/Process.h b/ProcessLib/Process.h
index 44df54f18ca60cc1a681558b8b6f13963aa97f9c..1f55b9e9b6411a4ded91aa7c88ba5d602eb446d2 100644
--- a/ProcessLib/Process.h
+++ b/ProcessLib/Process.h
@@ -139,13 +139,16 @@ private:
 
     virtual void assembleConcreteProcess(const double t, GlobalVector const& x,
                                          GlobalMatrix& M, GlobalMatrix& K,
-                                         GlobalVector& b) = 0;
+                                         GlobalVector& b,
+                                         StaggeredCouplingTerm const&
+                                         coupled_term) = 0;
 
     virtual void assembleWithJacobianConcreteProcess(
         const double t, GlobalVector const& x,
         GlobalVector const& xdot, const double dxdot_dx,
         const double dx_dx, GlobalMatrix& M, GlobalMatrix& K,
-        GlobalVector& b, GlobalMatrix& Jac) = 0;
+        GlobalVector& b, GlobalMatrix& Jac,
+        StaggeredCouplingTerm const& coupled_term) = 0;
 
     virtual void preTimestepConcreteProcess(GlobalVector const& /*x*/,
                                             const double /*t*/,
diff --git a/ProcessLib/RichardsFlow/RichardsFlowProcess.cpp b/ProcessLib/RichardsFlow/RichardsFlowProcess.cpp
index 0838ccf7d6e1d8e2123f465a87c416a9762b3fa5..f51c2d734e8ed4dabe749975742747fe09b93763 100644
--- a/ProcessLib/RichardsFlow/RichardsFlowProcess.cpp
+++ b/ProcessLib/RichardsFlow/RichardsFlowProcess.cpp
@@ -55,7 +55,9 @@ void RichardsFlowProcess::assembleConcreteProcess(const double t,
                                                   GlobalVector const& x,
                                                   GlobalMatrix& M,
                                                   GlobalMatrix& K,
-                                                  GlobalVector& b)
+                                                  GlobalVector& b,
+                                                  StaggeredCouplingTerm const&
+                                                  /*coupled_term*/)
 {
     DBUG("Assemble RichardsFlowProcess.");
 
@@ -68,7 +70,8 @@ void RichardsFlowProcess::assembleConcreteProcess(const double t,
 void RichardsFlowProcess::assembleWithJacobianConcreteProcess(
     const double t, GlobalVector const& x, GlobalVector const& xdot,
     const double dxdot_dx, const double dx_dx, GlobalMatrix& M, GlobalMatrix& K,
-    GlobalVector& b, GlobalMatrix& Jac)
+    GlobalVector& b, GlobalMatrix& Jac,
+    StaggeredCouplingTerm const& coupled_term)
 {
     DBUG("AssembleWithJacobian RichardsFlowProcess.");
 
@@ -76,7 +79,7 @@ void RichardsFlowProcess::assembleWithJacobianConcreteProcess(
     GlobalExecutor::executeMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
         _local_assemblers, *_local_to_global_index_map, t, x, xdot, dxdot_dx,
-        dx_dx, M, K, b, Jac);
+        dx_dx, M, K, b, Jac, coupled_term);
 }
 
 }  // namespace RichardsFlow
diff --git a/ProcessLib/RichardsFlow/RichardsFlowProcess.h b/ProcessLib/RichardsFlow/RichardsFlowProcess.h
index a0ec59047318d56bfad622f126d4149d19af0b8b..bcab9778a861dcd23b59d94c52d423e4a5551079 100644
--- a/ProcessLib/RichardsFlow/RichardsFlowProcess.h
+++ b/ProcessLib/RichardsFlow/RichardsFlowProcess.h
@@ -52,12 +52,15 @@ private:
 
     void assembleConcreteProcess(const double t, GlobalVector const& x,
                                  GlobalMatrix& M, GlobalMatrix& K,
-                                 GlobalVector& b) override;
+                                 GlobalVector& b,
+                                 StaggeredCouplingTerm const& coupled_term
+                                ) override;
 
     void assembleWithJacobianConcreteProcess(
         const double t, GlobalVector const& x, GlobalVector const& xdot,
         const double dxdot_dx, const double dx_dx, GlobalMatrix& M,
-        GlobalMatrix& K, GlobalVector& b, GlobalMatrix& Jac) override;
+        GlobalMatrix& K, GlobalVector& b, GlobalMatrix& Jac,
+        StaggeredCouplingTerm const& coupled_term) override;
 
     RichardsFlowProcessData _process_data;
 
diff --git a/ProcessLib/SmallDeformation/SmallDeformationProcess.h b/ProcessLib/SmallDeformation/SmallDeformationProcess.h
index 875f892c8fbaa57c8bfc9fbd03e950ce88faa60e..e0e4bf81528cce3acd2aba2ae28e0019b6a3df04 100644
--- a/ProcessLib/SmallDeformation/SmallDeformationProcess.h
+++ b/ProcessLib/SmallDeformation/SmallDeformationProcess.h
@@ -146,7 +146,9 @@ private:
 
     void assembleConcreteProcess(const double t, GlobalVector const& x,
                                  GlobalMatrix& M, GlobalMatrix& K,
-                                 GlobalVector& b) override
+                                 GlobalVector& b,
+                                 StaggeredCouplingTerm const&
+                                 /*coupled_term*/) override
     {
         DBUG("Assemble SmallDeformationProcess.");
 
@@ -159,7 +161,8 @@ private:
     void assembleWithJacobianConcreteProcess(
         const double t, GlobalVector const& x, GlobalVector const& xdot,
         const double dxdot_dx, const double dx_dx, GlobalMatrix& M,
-        GlobalMatrix& K, GlobalVector& b, GlobalMatrix& Jac) override
+        GlobalMatrix& K, GlobalVector& b, GlobalMatrix& Jac,
+        StaggeredCouplingTerm const& coupled_term) override
     {
         DBUG("AssembleWithJacobian SmallDeformationProcess.");
 
@@ -167,7 +170,7 @@ private:
         GlobalExecutor::executeMemberDereferenced(
             _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
             _local_assemblers, *_local_to_global_index_map, t, x, xdot,
-            dxdot_dx, dx_dx, M, K, b, Jac);
+            dxdot_dx, dx_dx, M, K, b, Jac, coupled_term);
     }
 
     void preTimestepConcreteProcess(GlobalVector const& x, double const t,
diff --git a/ProcessLib/StaggeredCouplingTerm.cpp b/ProcessLib/StaggeredCouplingTerm.cpp
index b588710308b800bc4e9a36b4ebee5923513dfe90..99e207857d53da8cea559e068b16a7dc022fef27 100644
--- a/ProcessLib/StaggeredCouplingTerm.cpp
+++ b/ProcessLib/StaggeredCouplingTerm.cpp
@@ -22,5 +22,12 @@ const StaggeredCouplingTerm createVoidStaggeredCouplingTerm()
     return StaggeredCouplingTerm(coupled_pcsss, coupled_xs);
 }
 
+const LocalCouplingTerm createVoidLocalCouplingTerm()
+{
+    std::map<ProcessType, Process const&> coupled_pcsss;
+    std::map<ProcessType, const std::vector<double>> local_coupled_xs;
+    return LocalCouplingTerm(coupled_pcsss, std::move(local_coupled_xs));
+}
+
 } // end of ProcessLib
 
diff --git a/ProcessLib/StaggeredCouplingTerm.h b/ProcessLib/StaggeredCouplingTerm.h
index c7a36166d8ee66fa3cfa0914e0076887b18a6186..76d056f14cba05808a5666617f6d97e11cf6105c 100644
--- a/ProcessLib/StaggeredCouplingTerm.h
+++ b/ProcessLib/StaggeredCouplingTerm.h
@@ -35,7 +35,22 @@ struct StaggeredCouplingTerm
     std::map<ProcessType, GlobalVector const*> const& coupled_xs;
 };
 
+struct LocalCouplingTerm
+{
+    LocalCouplingTerm(
+        std::map<ProcessType, Process const&> const& coupled_processes_,
+        std::map<ProcessType, const std::vector<double>>&& local_coupled_xs_)
+    : coupled_processes(coupled_processes_),
+      local_coupled_xs(std::move(local_coupled_xs_))
+    {
+    }
+
+    std::map<ProcessType, Process const&> const& coupled_processes;
+    std::map<ProcessType, const std::vector<double>> const local_coupled_xs;
+};
+
 const StaggeredCouplingTerm createVoidStaggeredCouplingTerm();
+const LocalCouplingTerm createVoidLocalCouplingTerm();
 
 } // end of ProcessLib
 #endif /* OGS_STAGGERED_COUPLING_TERM_H */
diff --git a/ProcessLib/TES/TESProcess.cpp b/ProcessLib/TES/TESProcess.cpp
index bad86ba65553db5d0867fb1798ce3a3397ce370b..c5f2ccc466ecfba70a9a8899bc44810a5b8bc723 100644
--- a/ProcessLib/TES/TESProcess.cpp
+++ b/ProcessLib/TES/TESProcess.cpp
@@ -226,7 +226,9 @@ void TESProcess::assembleConcreteProcess(const double t,
                                          GlobalVector const& x,
                                          GlobalMatrix& M,
                                          GlobalMatrix& K,
-                                         GlobalVector& b)
+                                         GlobalVector& b,
+                                         StaggeredCouplingTerm
+                                         const& /*coupled_term*/)
 {
     DBUG("Assemble TESProcess.");
 
@@ -239,12 +241,13 @@ void TESProcess::assembleConcreteProcess(const double t,
 void TESProcess::assembleWithJacobianConcreteProcess(
     const double t, GlobalVector const& x, GlobalVector const& xdot,
     const double dxdot_dx, const double dx_dx, GlobalMatrix& M, GlobalMatrix& K,
-    GlobalVector& b, GlobalMatrix& Jac)
+    GlobalVector& b, GlobalMatrix& Jac,
+    StaggeredCouplingTerm const& coupled_term)
 {
     GlobalExecutor::executeMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
         _local_assemblers, *_local_to_global_index_map, t, x, xdot, dxdot_dx,
-        dx_dx, M, K, b, Jac);
+        dx_dx, M, K, b, Jac, coupled_term);
 }
 
 void TESProcess::preTimestepConcreteProcess(GlobalVector const& x,
diff --git a/ProcessLib/TES/TESProcess.h b/ProcessLib/TES/TESProcess.h
index bed0ee6f6ac521a2d82d38afe0eb77f05b828d04..5a340fbb6c72f7c8bd2ba62e55a496bc78766138 100644
--- a/ProcessLib/TES/TESProcess.h
+++ b/ProcessLib/TES/TESProcess.h
@@ -59,14 +59,17 @@ private:
 
     void assembleConcreteProcess(const double t, GlobalVector const& x,
                                  GlobalMatrix& M, GlobalMatrix& K,
-                                 GlobalVector& b) override;
+                                 GlobalVector& b,
+                                 StaggeredCouplingTerm const& coupled_term
+                                ) override;
 
     void initializeSecondaryVariables();
 
     void assembleWithJacobianConcreteProcess(
         const double t, GlobalVector const& x, GlobalVector const& xdot,
         const double dxdot_dx, const double dx_dx, GlobalMatrix& M,
-        GlobalMatrix& K, GlobalVector& b, GlobalMatrix& Jac) override;
+        GlobalMatrix& K, GlobalVector& b, GlobalMatrix& Jac,
+        StaggeredCouplingTerm const& coupled_term) override;
 
     GlobalVector const& computeVapourPartialPressure(
         GlobalVector const& x,
diff --git a/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPProcess.cpp b/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPProcess.cpp
index 669fa71c3f6d74a4118f6aa66bab614ffc98a174..61e4c559d550733926bc515956ae186c4c9aad51 100644
--- a/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPProcess.cpp
+++ b/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPProcess.cpp
@@ -72,7 +72,9 @@ void TwoPhaseFlowWithPPProcess::assembleConcreteProcess(const double t,
                                                         GlobalVector const& x,
                                                         GlobalMatrix& M,
                                                         GlobalMatrix& K,
-                                                        GlobalVector& b)
+                                                        GlobalVector& b,
+                                                        StaggeredCouplingTerm
+                                                        const& /*coupled_term*/)
 {
     DBUG("Assemble TwoPhaseFlowWithPPProcess.");
     // Call global assembler for each local assembly item.
@@ -84,7 +86,8 @@ void TwoPhaseFlowWithPPProcess::assembleConcreteProcess(const double t,
 void TwoPhaseFlowWithPPProcess::assembleWithJacobianConcreteProcess(
     const double t, GlobalVector const& x, GlobalVector const& xdot,
     const double dxdot_dx, const double dx_dx, GlobalMatrix& M, GlobalMatrix& K,
-    GlobalVector& b, GlobalMatrix& Jac)
+    GlobalVector& b, GlobalMatrix& Jac,
+    StaggeredCouplingTerm const& coupled_term)
 {
     DBUG("AssembleWithJacobian TwoPhaseFlowWithPPProcess.");
 
@@ -92,7 +95,7 @@ void TwoPhaseFlowWithPPProcess::assembleWithJacobianConcreteProcess(
     GlobalExecutor::executeMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
         _local_assemblers, *_local_to_global_index_map, t, x, xdot, dxdot_dx,
-        dx_dx, M, K, b, Jac);
+        dx_dx, M, K, b, Jac, coupled_term);
 }
 
 }  // end of namespace
diff --git a/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPProcess.h b/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPProcess.h
index b3030be39876dfee803332f6d250658691325bc0..155deaef90169914a3d07e6b2383915297d58d5c 100644
--- a/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPProcess.h
+++ b/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPProcess.h
@@ -63,12 +63,15 @@ private:
 
     void assembleConcreteProcess(const double t, GlobalVector const& x,
                                  GlobalMatrix& M, GlobalMatrix& K,
-                                 GlobalVector& b) override;
+                                 GlobalVector& b,
+                                 StaggeredCouplingTerm const& coupled_term
+                                ) override;
 
     void assembleWithJacobianConcreteProcess(
         const double t, GlobalVector const& x, GlobalVector const& xdot,
         const double dxdot_dx, const double dx_dx, GlobalMatrix& M,
-        GlobalMatrix& K, GlobalVector& b, GlobalMatrix& Jac) override;
+        GlobalMatrix& K, GlobalVector& b, GlobalMatrix& Jac,
+        StaggeredCouplingTerm const& coupled_term) override;
 
     TwoPhaseFlowWithPPProcessData _process_data;
 
diff --git a/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoProcess.cpp b/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoProcess.cpp
index 2125e18045db29eab3a448e399572febb158a2f6..9af42f4c48f52d35a9a9ac71fc6e21bd05e69505 100644
--- a/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoProcess.cpp
+++ b/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoProcess.cpp
@@ -72,19 +72,22 @@ void TwoPhaseFlowWithPrhoProcess::assembleConcreteProcess(const double t,
                                                           GlobalVector const& x,
                                                           GlobalMatrix& M,
                                                           GlobalMatrix& K,
-                                                          GlobalVector& b)
+                                                          GlobalVector& b,
+                                                          StaggeredCouplingTerm
+                                                          const& coupled_term)
 {
     DBUG("Assemble TwoPhaseFlowWithPrhoProcess.");
     // Call global assembler for each local assembly item.
     GlobalExecutor::executeMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers,
-        *_local_to_global_index_map, t, x, M, K, b);
+        *_local_to_global_index_map, t, x, M, K, b, coupled_term);
 }
 
 void TwoPhaseFlowWithPrhoProcess::assembleWithJacobianConcreteProcess(
     const double t, GlobalVector const& x, GlobalVector const& xdot,
     const double dxdot_dx, const double dx_dx, GlobalMatrix& M, GlobalMatrix& K,
-    GlobalVector& b, GlobalMatrix& Jac)
+    GlobalVector& b, GlobalMatrix& Jac,
+    StaggeredCouplingTerm const& coupled_term)
 {
     DBUG("AssembleWithJacobian TwoPhaseFlowWithPrhoProcess.");
 
@@ -92,7 +95,7 @@ void TwoPhaseFlowWithPrhoProcess::assembleWithJacobianConcreteProcess(
     GlobalExecutor::executeMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
         _local_assemblers, *_local_to_global_index_map, t, x, xdot, dxdot_dx,
-        dx_dx, M, K, b, Jac);
+        dx_dx, M, K, b, Jac, coupled_term);
 }
 void TwoPhaseFlowWithPrhoProcess::preTimestepConcreteProcess(
     GlobalVector const& x, double const t, double const dt)
diff --git a/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoProcess.h b/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoProcess.h
index ccf2163602fa2908caa449cd5e9213e6bc1ee8dc..decfa1103e663efa3515c674d22944766574e7ed 100644
--- a/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoProcess.h
+++ b/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoProcess.h
@@ -56,12 +56,14 @@ private:
 
     void assembleConcreteProcess(const double t, GlobalVector const& x,
                                  GlobalMatrix& M, GlobalMatrix& K,
-                                 GlobalVector& b) override;
+                                 GlobalVector& b,
+                                 StaggeredCouplingTerm const& coupled_term) override;
 
     void assembleWithJacobianConcreteProcess(
         const double t, GlobalVector const& x, GlobalVector const& xdot,
         const double dxdot_dx, const double dx_dx, GlobalMatrix& M,
-        GlobalMatrix& K, GlobalVector& b, GlobalMatrix& Jac) override;
+        GlobalMatrix& K, GlobalVector& b, GlobalMatrix& Jac,
+        StaggeredCouplingTerm const& coupled_term) override;
 
     void preTimestepConcreteProcess(GlobalVector const& x, const double t,
                                     const double delta_t) override;
diff --git a/ProcessLib/VectorMatrixAssembler.cpp b/ProcessLib/VectorMatrixAssembler.cpp
index dda159cf6ab5061d69c2537c4db04d0b1899ef38..00628c642e7a208b8a40907dd779cf1435d45190 100644
--- a/ProcessLib/VectorMatrixAssembler.cpp
+++ b/ProcessLib/VectorMatrixAssembler.cpp
@@ -15,6 +15,8 @@
 #include "MathLib/LinAlg/Eigen/EigenMapTools.h"
 #include "LocalAssemblerInterface.h"
 
+#include "Process.h"
+
 namespace ProcessLib
 {
 VectorMatrixAssembler::VectorMatrixAssembler(
@@ -26,7 +28,8 @@ VectorMatrixAssembler::VectorMatrixAssembler(
 void VectorMatrixAssembler::assemble(
     const std::size_t mesh_item_id, LocalAssemblerInterface& local_assembler,
     const NumLib::LocalToGlobalIndexMap& dof_table, const double t,
-    const GlobalVector& x, GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b)
+    const GlobalVector& x, GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b,
+    const StaggeredCouplingTerm& coupled_term)
 {
     auto const indices = NumLib::getIndices(mesh_item_id, dof_table);
     auto const local_x = x.get(indices);
@@ -35,21 +38,50 @@ void VectorMatrixAssembler::assemble(
     _local_K_data.clear();
     _local_b_data.clear();
 
-    local_assembler.assemble(t, local_x, _local_M_data, _local_K_data, _local_b_data);
+    if (coupled_term.coupled_processes.empty())
+    {
+        local_assembler.assemble(t, local_x, _local_M_data, _local_K_data,
+                                 _local_b_data,
+                                 ProcessLib::createVoidLocalCouplingTerm());
+    }
+    else
+    {
+        std::map<ProcessLib::ProcessType, const std::vector<double>>
+            local_coupled_xs;
+        auto it = coupled_term.coupled_xs.begin();
+        while (it != coupled_term.coupled_xs.end())
+        {
+            GlobalVector const* coupled_x = it->second;
+            auto const local_coupled_x = coupled_x->get(indices);
+            BaseLib::insertMapIfKeyUniqueElseError(local_coupled_xs, it->first,
+                                                   local_coupled_x,
+                                                   "local_coupled_x");
+            it++;
+        }
+
+        ProcessLib::LocalCouplingTerm local_coupling_term(
+            coupled_term.coupled_processes, std::move(local_coupled_xs));
+
+        local_assembler.assemble(t, local_x, _local_M_data, _local_K_data,
+                                 _local_b_data, local_coupling_term);
+    }
 
     auto const num_r_c = indices.size();
     auto const r_c_indices =
         NumLib::LocalToGlobalIndexMap::RowColumnIndices(indices, indices);
 
-    if (!_local_M_data.empty()) {
+    if (!_local_M_data.empty())
+    {
         auto const local_M = MathLib::toMatrix(_local_M_data, num_r_c, num_r_c);
         M.add(r_c_indices, local_M);
     }
-    if (!_local_K_data.empty()) {
+    if (!_local_K_data.empty())
+    {
         auto const local_K = MathLib::toMatrix(_local_K_data, num_r_c, num_r_c);
         K.add(r_c_indices, local_K);
     }
-    if (!_local_b_data.empty()) {
+    if (!_local_b_data.empty())
+    {
         assert(_local_b_data.size() == num_r_c);
         b.add(indices, _local_b_data);
     }
@@ -60,7 +92,7 @@ void VectorMatrixAssembler::assembleWithJacobian(
     NumLib::LocalToGlobalIndexMap const& dof_table, const double t,
     GlobalVector const& x, GlobalVector const& xdot, const double dxdot_dx,
     const double dx_dx, GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b,
-    GlobalMatrix& Jac)
+    GlobalMatrix& Jac, const StaggeredCouplingTerm& coupled_term)
 {
     auto const indices = NumLib::getIndices(mesh_item_id, dof_table);
     auto const local_x = x.get(indices);
@@ -71,31 +103,64 @@ void VectorMatrixAssembler::assembleWithJacobian(
     _local_b_data.clear();
     _local_Jac_data.clear();
 
-    _jacobian_assembler->assembleWithJacobian(
-        local_assembler, t, local_x, local_xdot, dxdot_dx, dx_dx, _local_M_data,
-        _local_K_data, _local_b_data, _local_Jac_data);
+    if (coupled_term.coupled_processes.empty())
+    {
+        _jacobian_assembler->assembleWithJacobian(
+            local_assembler, t, local_x, local_xdot, dxdot_dx, dx_dx,
+            _local_M_data, _local_K_data, _local_b_data, _local_Jac_data,
+            ProcessLib::createVoidLocalCouplingTerm());
+    }
+    else
+    {
+        std::map<ProcessLib::ProcessType, const std::vector<double>>
+            local_coupled_xs;
+        auto it = coupled_term.coupled_xs.begin();
+        while (it != coupled_term.coupled_xs.end())
+        {
+            GlobalVector const* coupled_x = it->second;
+            auto const local_coupled_x = coupled_x->get(indices);
+            BaseLib::insertMapIfKeyUniqueElseError(local_coupled_xs, it->first,
+                                                   local_coupled_x,
+                                                   "local_coupled_x");
+            it++;
+        }
+
+        ProcessLib::LocalCouplingTerm local_coupling_term(
+            coupled_term.coupled_processes, std::move(local_coupled_xs));
+
+        _jacobian_assembler->assembleWithJacobian(
+            local_assembler, t, local_x, local_xdot, dxdot_dx, dx_dx,
+            _local_M_data, _local_K_data, _local_b_data, _local_Jac_data,
+            local_coupling_term);
+    }
 
     auto const num_r_c = indices.size();
     auto const r_c_indices =
         NumLib::LocalToGlobalIndexMap::RowColumnIndices(indices, indices);
 
-    if (!_local_M_data.empty()) {
+    if (!_local_M_data.empty())
+    {
         auto const local_M = MathLib::toMatrix(_local_M_data, num_r_c, num_r_c);
         M.add(r_c_indices, local_M);
     }
-    if (!_local_K_data.empty()) {
+    if (!_local_K_data.empty())
+    {
         auto const local_K = MathLib::toMatrix(_local_K_data, num_r_c, num_r_c);
         K.add(r_c_indices, local_K);
     }
-    if (!_local_b_data.empty()) {
+    if (!_local_b_data.empty())
+    {
         assert(_local_b_data.size() == num_r_c);
         b.add(indices, _local_b_data);
     }
-    if (!_local_Jac_data.empty()) {
+    if (!_local_Jac_data.empty())
+    {
         auto const local_Jac =
             MathLib::toMatrix(_local_Jac_data, num_r_c, num_r_c);
         Jac.add(r_c_indices, local_Jac);
-    } else {
+    }
+    else
+    {
         OGS_FATAL(
             "No Jacobian has been assembled! This might be due to programming "
             "errors in the local assembler of the current process.");
diff --git a/ProcessLib/VectorMatrixAssembler.h b/ProcessLib/VectorMatrixAssembler.h
index a017a00ab0050b4420d898588e0bf3a2106ee72c..5466d6d71c008aff6937f077fcf5f3038e7e3020 100644
--- a/ProcessLib/VectorMatrixAssembler.h
+++ b/ProcessLib/VectorMatrixAssembler.h
@@ -12,6 +12,7 @@
 #include <vector>
 #include "NumLib/NumericsConfig.h"
 #include "AbstractJacobianAssembler.h"
+#include "StaggeredCouplingTerm.h"
 
 namespace NumLib
 {
@@ -38,7 +39,8 @@ public:
                   LocalAssemblerInterface& local_assembler,
                   NumLib::LocalToGlobalIndexMap const& dof_table,
                   double const t, GlobalVector const& x, GlobalMatrix& M,
-                  GlobalMatrix& K, GlobalVector& b);
+                  GlobalMatrix& K, GlobalVector& b,
+                  const StaggeredCouplingTerm& coupled_term);
 
     //! Assembles \c M, \c K, \c b, and the Jacobian \c Jac of the residual.
     //! \note The Jacobian must be assembled.
@@ -49,7 +51,8 @@ public:
                               GlobalVector const& xdot, const double dxdot_dx,
                               const double dx_dx, GlobalMatrix& M,
                               GlobalMatrix& K, GlobalVector& b,
-                              GlobalMatrix& Jac);
+                              GlobalMatrix& Jac,
+                              const StaggeredCouplingTerm& coupled_term);
 
 private:
     // temporary data only stored here in order to avoid frequent memory