diff --git a/ProcessLib/AbstractJacobianAssembler.h b/ProcessLib/AbstractJacobianAssembler.h
index 7989d383b16673deb59c37180702ddc72422efd2..b4267fadb7ab0d13e8a280dbf30e554d9bb7bffa 100644
--- a/ProcessLib/AbstractJacobianAssembler.h
+++ b/ProcessLib/AbstractJacobianAssembler.h
@@ -14,6 +14,7 @@
 namespace ProcessLib
 {
 class LocalAssemblerInterface;
+struct LocalCouplingTerm;
 
 //! Base class for Jacobian assemblers.
 class AbstractJacobianAssembler
@@ -29,6 +30,18 @@ public:
         std::vector<double>& local_K_data, std::vector<double>& local_b_data,
         std::vector<double>& local_Jac_data) = 0;
 
+    //! Assembles the Jacobian, the matrices \f$M\f$ and \f$K\f$, and the vector
+    //! \f$b\f$ with coupling.
+    virtual void coupling_assembleWithJacobian(
+        LocalAssemblerInterface& /*local_assembler*/, double const /*t*/,
+        std::vector<double> const& /*local_x*/,
+        std::vector<double> const& /*local_xdot*/, const double /*dxdot_dx*/,
+        const double /*dx_dx*/, std::vector<double>& /*local_M_data*/,
+        std::vector<double>& /*local_K_data*/,
+        std::vector<double>& /*local_b_data*/,
+        std::vector<double>& /*local_Jac_data*/,
+        LocalCouplingTerm const& /*coupled_term*/) {}
+
     virtual ~AbstractJacobianAssembler() = default;
 };
 
diff --git a/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.cpp b/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.cpp
index 93589087db50fc1df5fbfb73451f19a4262b4186..11370da33e06a9ba8cd316f7156f6b4111b37f8a 100644
--- a/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.cpp
+++ b/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.cpp
@@ -78,14 +78,14 @@ void GroundwaterFlowProcess::assembleConcreteProcess(const double t,
                                                      GlobalMatrix& K,
                                                      GlobalVector& b,
                                                      StaggeredCouplingTerm
-                                                     const& /*coupled_term*/)
+                                                     const& coupled_term)
 {
     DBUG("Assemble GroundwaterFlowProcess.");
 
     // 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 GroundwaterFlowProcess::assembleWithJacobianConcreteProcess(
diff --git a/ProcessLib/HT/HTProcess.cpp b/ProcessLib/HT/HTProcess.cpp
index d8f4eddb63aa472d42d75701784dc1d0ac6c6aee..a9d0b6d19e13b5a30a0239929c67a1b769fa2c53 100644
--- a/ProcessLib/HT/HTProcess.cpp
+++ b/ProcessLib/HT/HTProcess.cpp
@@ -73,14 +73,14 @@ void HTProcess::assembleConcreteProcess(const double t,
                                         GlobalMatrix& K,
                                         GlobalVector& b,
                                         StaggeredCouplingTerm const&
-                                        /*coupled_term*/)
+                                        coupled_term)
 {
     DBUG("Assemble HTProcess.");
 
     // 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 HTProcess::assembleWithJacobianConcreteProcess(
diff --git a/ProcessLib/HeatConduction/HeatConductionProcess.cpp b/ProcessLib/HeatConduction/HeatConductionProcess.cpp
index 3ceed6b6453b7e9534bab6a6f26877541c8511ed..d181c782f1071adcc7937a5a306c9a79cb73ff55 100644
--- a/ProcessLib/HeatConduction/HeatConductionProcess.cpp
+++ b/ProcessLib/HeatConduction/HeatConductionProcess.cpp
@@ -74,14 +74,14 @@ void HeatConductionProcess::assembleConcreteProcess(const double t,
                                                     GlobalMatrix& K,
                                                     GlobalVector& b,
                                                     StaggeredCouplingTerm
-                                                    const& /*coupled_term*/)
+                                                    const& coupled_term)
 {
     DBUG("Assemble HeatConductionProcess.");
 
     // 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 HeatConductionProcess::assembleWithJacobianConcreteProcess(
diff --git a/ProcessLib/HydroMechanics/HydroMechanicsFEM.h b/ProcessLib/HydroMechanics/HydroMechanicsFEM.h
index 72540162dfb8c12c71127679e1f74ed4493f8b1e..97e3f0d0eabb8632963620caf89e07a81c1c4f30 100644
--- a/ProcessLib/HydroMechanics/HydroMechanicsFEM.h
+++ b/ProcessLib/HydroMechanics/HydroMechanicsFEM.h
@@ -200,9 +200,7 @@ 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,
-                              LocalCouplingTerm const& /*coupled_term*/
-                             ) override
+                              std::vector<double>& local_Jac_data) override
     {
         assert(local_x.size() == pressure_size + displacement_size);
 
diff --git a/ProcessLib/HydroMechanics/HydroMechanicsProcess.h b/ProcessLib/HydroMechanics/HydroMechanicsProcess.h
index 0a6dc8c707a9cdb46bda71ff4271b763205b8f2c..88d59391b9a2d9efc3917113f091aedf78fddc85 100644
--- a/ProcessLib/HydroMechanics/HydroMechanicsProcess.h
+++ b/ProcessLib/HydroMechanics/HydroMechanicsProcess.h
@@ -111,14 +111,15 @@ private:
                                  GlobalMatrix& M, GlobalMatrix& K,
                                  GlobalVector& b,
                                  StaggeredCouplingTerm const&
-                                 /*coupled_term*/) override
+                                 coupled_term) override
     {
         DBUG("Assemble HydroMechanicsProcess.");
 
         // 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_assemblers, *_local_to_global_index_map, t, x, M, K, b,
+            coupled_term);
     }
 
     void assembleWithJacobianConcreteProcess(
diff --git a/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcess.h b/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcess.h
index 18f260b49bb51aaef7bdc8906af34a8b155a5358..8638c0fbf6e5c3ccbc1ad401390bef854743fc09 100644
--- a/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcess.h
+++ b/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcess.h
@@ -72,20 +72,22 @@ private:
                                  GlobalMatrix& M, GlobalMatrix& K,
                                  GlobalVector& b,
                                  StaggeredCouplingTerm const&
-                                 /*coupled_term*/) override
+                                 coupled_term) override
     {
         DBUG("Assemble HydroMechanicsProcess.");
 
         // 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_assemblers, *_local_to_global_index_map, t, x, M, K, b,
+            coupled_term);
     }
 
     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 HydroMechanicsProcess.");
 
@@ -93,7 +95,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/SmallDeformation/SmallDeformationProcess.h b/ProcessLib/LIE/SmallDeformation/SmallDeformationProcess.h
index efd4da63603719d69ab8fc87db2c9345a1f34d51..52fe3451935661e1602da6417be2c308f264c6bd 100644
--- a/ProcessLib/LIE/SmallDeformation/SmallDeformationProcess.h
+++ b/ProcessLib/LIE/SmallDeformation/SmallDeformationProcess.h
@@ -70,14 +70,15 @@ private:
                                  GlobalMatrix& M, GlobalMatrix& K,
                                  GlobalVector& b,
                                  StaggeredCouplingTerm const&
-                                 /*coupled_term*/) override
+                                 coupled_term) override
     {
         DBUG("Assemble SmallDeformationProcess.");
 
         // 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_assemblers, *_local_to_global_index_map, t, x, M, K, b,
+            coupled_term);
     }
 
     void assembleWithJacobianConcreteProcess(
diff --git a/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp b/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp
index b5f351e009a37199b06c3b3e64dd2a6ed44ab4b0..e9598f66f35bdc0256c16b8908b2c6899f85aefb 100644
--- a/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp
+++ b/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp
@@ -94,13 +94,13 @@ void LiquidFlowProcess::assembleConcreteProcess(const double t,
                                                 GlobalMatrix& K,
                                                 GlobalVector& b,
                                                 StaggeredCouplingTerm const&
-                                                /*coupled_term*/)
+                                                coupled_term)
 {
     DBUG("Assemble LiquidFlowProcess.");
     // 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 LiquidFlowProcess::assembleWithJacobianConcreteProcess(
diff --git a/ProcessLib/LocalAssemblerInterface.cpp b/ProcessLib/LocalAssemblerInterface.cpp
index 508f73b6bfad4ed17d8c9981cf90738e2862e7bc..ab83fdaef16c54192e507055b4b2a9c5d4a63f43 100644
--- a/ProcessLib/LocalAssemblerInterface.cpp
+++ b/ProcessLib/LocalAssemblerInterface.cpp
@@ -13,6 +13,19 @@
 
 namespace ProcessLib
 {
+
+void LocalAssemblerInterface::coupling_assemble(
+    double const /*t*/, std::vector<double> const& /*local_x*/,
+    std::vector<double>& /*local_M_data*/,
+    std::vector<double>& /*local_K_data*/,
+    std::vector<double>& /*local_b_data*/,
+    LocalCouplingTerm const& /*coupled_term*/)
+{
+    OGS_FATAL(
+        "The coupling_assemble() function is not implemented in the local "
+        "assembler.");
+}
+
 void LocalAssemblerInterface::assembleWithJacobian(
     double const /*t*/, std::vector<double> const& /*local_x*/,
     std::vector<double> const& /*local_xdot*/, const double /*dxdot_dx*/,
@@ -26,6 +39,20 @@ void LocalAssemblerInterface::assembleWithJacobian(
         "assembler.");
 }
 
+void LocalAssemblerInterface::coupling_assembleWithJacobian(
+    double const /*t*/, std::vector<double> const& /*local_x*/,
+    std::vector<double> const& /*local_xdot*/, const double /*dxdot_dx*/,
+    const double /*dx_dx*/, std::vector<double>& /*local_M_data*/,
+    std::vector<double>& /*local_K_data*/,
+    std::vector<double>& /*local_b_data*/,
+    std::vector<double>& /*local_Jac_data*/,
+    LocalCouplingTerm const& /*coupled_term*/)
+{
+    OGS_FATAL(
+        "The coupling_assembleWithJacobian() function is not implemented in"
+        " the local.");
+}
+
 void LocalAssemblerInterface::computeSecondaryVariable(
                               std::size_t const mesh_item_id,
                               NumLib::LocalToGlobalIndexMap const& dof_table,
diff --git a/ProcessLib/LocalAssemblerInterface.h b/ProcessLib/LocalAssemblerInterface.h
index 8176485f679d31959f251bc8d869ac0691917ed9..aeb7329343ddb73af011d18f51c260780f2b5ea9 100644
--- a/ProcessLib/LocalAssemblerInterface.h
+++ b/ProcessLib/LocalAssemblerInterface.h
@@ -11,6 +11,7 @@
 
 #include "NumLib/NumericsConfig.h"
 #include "MathLib/Point3d.h"
+#include "StaggeredCouplingTerm.h"
 
 namespace NumLib
 {
@@ -35,6 +36,12 @@ public:
         std::vector<double>& local_M_data, std::vector<double>& local_K_data,
         std::vector<double>& local_b_data) = 0;
 
+    virtual void coupling_assemble(
+        double const t, std::vector<double> const& local_x,
+        std::vector<double>& local_M_data, std::vector<double>& local_K_data,
+        std::vector<double>& local_b_data,
+        LocalCouplingTerm const& coupled_term);
+
     virtual void assembleWithJacobian(double const t,
                                       std::vector<double> const& local_x,
                                       std::vector<double> const& local_xdot,
@@ -44,6 +51,16 @@ public:
                                       std::vector<double>& local_b_data,
                                       std::vector<double>& local_Jac_data);
 
+    virtual void coupling_assembleWithJacobian(double const t,
+                                      std::vector<double> const& local_x,
+                                      std::vector<double> const& local_xdot,
+                                      const double dxdot_dx, const double dx_dx,
+                                      std::vector<double>& local_M_data,
+                                      std::vector<double>& local_K_data,
+                                      std::vector<double>& local_b_data,
+                                      std::vector<double>& local_Jac_data,
+                                      LocalCouplingTerm const& coupled_term);
+
     virtual void computeSecondaryVariable(std::size_t const mesh_item_id,
                               NumLib::LocalToGlobalIndexMap const& dof_table,
                               const double t, GlobalVector const& x);
diff --git a/ProcessLib/RichardsFlow/RichardsFlowProcess.cpp b/ProcessLib/RichardsFlow/RichardsFlowProcess.cpp
index f51c2d734e8ed4dabe749975742747fe09b93763..7e5858a6dc4f32d2177788efab9c5ca10bbc1fb5 100644
--- a/ProcessLib/RichardsFlow/RichardsFlowProcess.cpp
+++ b/ProcessLib/RichardsFlow/RichardsFlowProcess.cpp
@@ -57,14 +57,14 @@ void RichardsFlowProcess::assembleConcreteProcess(const double t,
                                                   GlobalMatrix& K,
                                                   GlobalVector& b,
                                                   StaggeredCouplingTerm const&
-                                                  /*coupled_term*/)
+                                                  coupled_term)
 {
     DBUG("Assemble RichardsFlowProcess.");
 
     // 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 RichardsFlowProcess::assembleWithJacobianConcreteProcess(
diff --git a/ProcessLib/SmallDeformation/SmallDeformationProcess.h b/ProcessLib/SmallDeformation/SmallDeformationProcess.h
index e0e4bf81528cce3acd2aba2ae28e0019b6a3df04..c94a7b001a4e90d902e0f18803771be087974f12 100644
--- a/ProcessLib/SmallDeformation/SmallDeformationProcess.h
+++ b/ProcessLib/SmallDeformation/SmallDeformationProcess.h
@@ -148,14 +148,15 @@ private:
                                  GlobalMatrix& M, GlobalMatrix& K,
                                  GlobalVector& b,
                                  StaggeredCouplingTerm const&
-                                 /*coupled_term*/) override
+                                 coupled_term) override
     {
         DBUG("Assemble SmallDeformationProcess.");
 
         // 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_assemblers, *_local_to_global_index_map, t, x, M, K, b,
+            coupled_term);
     }
 
     void assembleWithJacobianConcreteProcess(
diff --git a/ProcessLib/TES/TESProcess.cpp b/ProcessLib/TES/TESProcess.cpp
index c5f2ccc466ecfba70a9a8899bc44810a5b8bc723..1b9c77a9b9b27074f61ac66c0990808933f4d692 100644
--- a/ProcessLib/TES/TESProcess.cpp
+++ b/ProcessLib/TES/TESProcess.cpp
@@ -228,14 +228,14 @@ void TESProcess::assembleConcreteProcess(const double t,
                                          GlobalMatrix& K,
                                          GlobalVector& b,
                                          StaggeredCouplingTerm
-                                         const& /*coupled_term*/)
+                                         const& coupled_term)
 {
     DBUG("Assemble TESProcess.");
 
     // 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 TESProcess::assembleWithJacobianConcreteProcess(
diff --git a/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPProcess.cpp b/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPProcess.cpp
index 61e4c559d550733926bc515956ae186c4c9aad51..79b9619105fab9bbb69b32c9b0acb06151f8d76b 100644
--- a/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPProcess.cpp
+++ b/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPProcess.cpp
@@ -74,13 +74,13 @@ void TwoPhaseFlowWithPPProcess::assembleConcreteProcess(const double t,
                                                         GlobalMatrix& K,
                                                         GlobalVector& b,
                                                         StaggeredCouplingTerm
-                                                        const& /*coupled_term*/)
+                                                        const& coupled_term)
 {
     DBUG("Assemble TwoPhaseFlowWithPPProcess.");
     // 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 TwoPhaseFlowWithPPProcess::assembleWithJacobianConcreteProcess(
diff --git a/ProcessLib/VectorMatrixAssembler.cpp b/ProcessLib/VectorMatrixAssembler.cpp
index 00628c642e7a208b8a40907dd779cf1435d45190..313ac0e03016601e23140bc62b8e4f14045b8a64 100644
--- a/ProcessLib/VectorMatrixAssembler.cpp
+++ b/ProcessLib/VectorMatrixAssembler.cpp
@@ -41,8 +41,7 @@ void VectorMatrixAssembler::assemble(
     if (coupled_term.coupled_processes.empty())
     {
         local_assembler.assemble(t, local_x, _local_M_data, _local_K_data,
-                                 _local_b_data,
-                                 ProcessLib::createVoidLocalCouplingTerm());
+                                 _local_b_data);
     }
     else
     {
@@ -62,8 +61,9 @@ void VectorMatrixAssembler::assemble(
         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);
+        local_assembler.coupling_assemble(t, local_x, _local_M_data,
+                                          _local_K_data, _local_b_data,
+                                          local_coupling_term);
     }
 
     auto const num_r_c = indices.size();
@@ -107,8 +107,7 @@ void VectorMatrixAssembler::assembleWithJacobian(
     {
         _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());
+            _local_M_data, _local_K_data, _local_b_data, _local_Jac_data);
     }
     else
     {
@@ -128,7 +127,7 @@ void VectorMatrixAssembler::assembleWithJacobian(
         ProcessLib::LocalCouplingTerm local_coupling_term(
             coupled_term.coupled_processes, std::move(local_coupled_xs));
 
-        _jacobian_assembler->assembleWithJacobian(
+        _jacobian_assembler->coupling_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);