diff --git a/ProcessLib/HT/HTProcess.cpp b/ProcessLib/HT/HTProcess.cpp
index d00b04b188d81d5df122d7cf78d75e3bdf627423..026847fb8cda2fc150ff4457434c8c31929d4e2b 100644
--- a/ProcessLib/HT/HTProcess.cpp
+++ b/ProcessLib/HT/HTProcess.cpp
@@ -75,11 +75,14 @@ void HTProcess::initializeConcreteProcess(
     }
     else
     {
+        const int heat_transport_process_id = 0;
+        const int hydraulic_process_id = 1;
+
         ProcessLib::createLocalAssemblers<StaggeredHTFEM>(
             mesh.getDimension(), mesh.getElements(), dof_table,
             pv.getShapeFunctionOrder(), _local_assemblers,
-            mesh.isAxiallySymmetric(), integration_order,
-            *_material_properties);
+            mesh.isAxiallySymmetric(), integration_order, *_material_properties,
+            heat_transport_process_id, hydraulic_process_id);
     }
 
     _secondary_variables.addSecondaryVariable(
diff --git a/ProcessLib/HT/StaggeredHTFEM-impl.h b/ProcessLib/HT/StaggeredHTFEM-impl.h
index 7c12558ec2bf51da46ad34b52f92f7c62832457a..e457b733714dd0d9b0679ae275331d22f39f08bb 100644
--- a/ProcessLib/HT/StaggeredHTFEM-impl.h
+++ b/ProcessLib/HT/StaggeredHTFEM-impl.h
@@ -28,7 +28,7 @@ void StaggeredHTFEM<ShapeFunction, IntegrationMethod, GlobalDim>::
                                std::vector<double>& local_b_data,
                                LocalCoupledSolutions const& coupled_xs)
 {
-    if (coupled_xs.process_id == 0)
+    if (coupled_xs.process_id == _heat_transport_process_id)
     {
         assembleHeatTransportEquation(t, local_M_data, local_K_data,
                                       local_b_data, coupled_xs);
@@ -47,14 +47,16 @@ void StaggeredHTFEM<ShapeFunction, IntegrationMethod, GlobalDim>::
                               std::vector<double>& local_b_data,
                               LocalCoupledSolutions const& coupled_xs)
 {
-    auto const& local_p = coupled_xs.local_coupled_xs[1];
+    auto const& local_p = coupled_xs.local_coupled_xs[_hydraulic_process_id];
     auto const local_matrix_size = local_p.size();
     // This assertion is valid only if all nodal d.o.f. use the same shape
     // matrices.
     assert(local_matrix_size == ShapeFunction::NPOINTS);
 
-    auto const& local_T1 = coupled_xs.local_coupled_xs[0];
-    auto const& local_T0 = coupled_xs.local_coupled_xs0[0];
+    auto const& local_T1 =
+        coupled_xs.local_coupled_xs[_heat_transport_process_id];
+    auto const& local_T0 =
+        coupled_xs.local_coupled_xs0[_heat_transport_process_id];
     const double dt = coupled_xs.dt;
 
     auto local_M = MathLib::createZeroedMatrix<LocalMatrixType>(
@@ -126,9 +128,10 @@ void StaggeredHTFEM<ShapeFunction, IntegrationMethod, GlobalDim>::
         GlobalDimMatrixType K_over_mu = intrinsic_permeability / viscosity;
 
         // matrix assembly
-        local_M.noalias() += w * (porosity * dfluid_density_dp / fluid_density +
-                                  specific_storage) *
-                             N.transpose() * N;
+        local_M.noalias() +=
+            w *
+            (porosity * dfluid_density_dp / fluid_density + specific_storage) *
+            N.transpose() * N;
 
         local_K.noalias() += w * dNdx.transpose() * K_over_mu * dNdx;
 
@@ -171,7 +174,7 @@ void StaggeredHTFEM<ShapeFunction, IntegrationMethod, GlobalDim>::
                                   std::vector<double>& /*local_b_data*/,
                                   LocalCoupledSolutions const& coupled_xs)
 {
-    auto const& local_p = coupled_xs.local_coupled_xs[1];
+    auto const& local_p = coupled_xs.local_coupled_xs[_hydraulic_process_id];
     auto const local_matrix_size = local_p.size();
     // This assertion is valid only if all nodal d.o.f. use the same shape
     // matrices.
@@ -180,7 +183,8 @@ void StaggeredHTFEM<ShapeFunction, IntegrationMethod, GlobalDim>::
     auto local_p_Eigen_type =
         Eigen::Map<const NodalVectorType>(&local_p[0], local_matrix_size);
 
-    auto const& local_T1 = coupled_xs.local_coupled_xs[0];
+    auto const& local_T1 =
+        coupled_xs.local_coupled_xs[_heat_transport_process_id];
 
     auto local_M = MathLib::createZeroedMatrix<LocalMatrixType>(
         local_M_data, local_matrix_size, local_matrix_size);
@@ -282,7 +286,9 @@ StaggeredHTFEM<ShapeFunction, IntegrationMethod, GlobalDim>::
     auto const local_xs = getCurrentLocalSolutions(
         *(this->_coupled_solutions), indices_of_all_coupled_processes);
 
-    return this->getIntPtDarcyVelocityLocal(t, local_xs[0], local_xs[1], cache);
+    return this->getIntPtDarcyVelocityLocal(
+        t, local_xs[_hydraulic_process_id],
+        local_xs[_heat_transport_process_id], cache);
 }
 }  // namespace HT
 }  // namespace ProcessLib
diff --git a/ProcessLib/HT/StaggeredHTFEM.h b/ProcessLib/HT/StaggeredHTFEM.h
index d9296be17add4161be26819fc6ae357818240d33..d54700bfeb9231ab07de5256eb1f0a23ba6de30d 100644
--- a/ProcessLib/HT/StaggeredHTFEM.h
+++ b/ProcessLib/HT/StaggeredHTFEM.h
@@ -55,10 +55,14 @@ public:
                    std::size_t const local_matrix_size,
                    bool is_axially_symmetric,
                    unsigned const integration_order,
-                   HTMaterialProperties const& material_properties)
+                   HTMaterialProperties const& material_properties,
+                   const int heat_transport_process_id,
+                   const int hydraulic_process_id)
         : HTFEM<ShapeFunction, IntegrationMethod, GlobalDim>(
               element, local_matrix_size, is_axially_symmetric,
-              integration_order, material_properties, 1)
+              integration_order, material_properties, 1),
+        _heat_transport_process_id(heat_transport_process_id),
+        _hydraulic_process_id(hydraulic_process_id)
     {
     }
 
@@ -84,6 +88,8 @@ private:
         double const t, std::vector<double>& local_M_data,
         std::vector<double>& local_K_data, std::vector<double>& local_b_data,
         LocalCoupledSolutions const& coupled_xs);
+    const int _heat_transport_process_id;
+    const int _hydraulic_process_id;
 };
 
 }  // namespace HT