diff --git a/ProcessLib/HeatTransportBHE/BHE/BHECommonCoaxial.cpp b/ProcessLib/HeatTransportBHE/BHE/BHECommonCoaxial.cpp
index 945ae14126ed2dd45357a24e0e9cc76b0ea02205..5f3b2cc6627a154d19f2bf7a35a6039d5cad1b56 100644
--- a/ProcessLib/HeatTransportBHE/BHE/BHECommonCoaxial.cpp
+++ b/ProcessLib/HeatTransportBHE/BHE/BHECommonCoaxial.cpp
@@ -20,6 +20,26 @@ namespace HeatTransportBHE
 {
 namespace BHE
 {
+BHECommonCoaxial::BHECommonCoaxial(
+    BoreholeGeometry const& borehole,
+    RefrigerantProperties const& refrigerant,
+    GroutParameters const& grout,
+    FlowAndTemperatureControl const& flowAndTemperatureControl,
+    PipeConfigurationCoaxial const& pipes,
+    bool const use_python_bcs)
+    : BHECommon{borehole, refrigerant, grout, flowAndTemperatureControl,
+                use_python_bcs},
+      _pipes(pipes)
+{
+    cross_section_area_inner_pipe = _pipes.inner_pipe.area();
+    cross_section_area_annulus =
+        _pipes.outer_pipe.area() - _pipes.inner_pipe.outsideArea();
+    cross_section_area_grout =
+        borehole_geometry.area() - _pipes.outer_pipe.outsideArea();
+
+    _thermal_resistances.fill(std::numeric_limits<double>::quiet_NaN());
+}
+
 std::array<double, BHECommonCoaxial::number_of_unknowns>
 BHECommonCoaxial::pipeHeatCapacities() const
 {
diff --git a/ProcessLib/HeatTransportBHE/BHE/BHECommonCoaxial.h b/ProcessLib/HeatTransportBHE/BHE/BHECommonCoaxial.h
index 6b42e434001103120188432dd7a7411b5932c1f9..b033891a34e223708f2c39ffeeab2b1b1f4aa31a 100644
--- a/ProcessLib/HeatTransportBHE/BHE/BHECommonCoaxial.h
+++ b/ProcessLib/HeatTransportBHE/BHE/BHECommonCoaxial.h
@@ -28,19 +28,7 @@ public:
                      GroutParameters const& grout,
                      FlowAndTemperatureControl const& flowAndTemperatureControl,
                      PipeConfigurationCoaxial const& pipes,
-                     bool const use_python_bcs)
-        : BHECommon{borehole, refrigerant, grout, flowAndTemperatureControl,
-                    use_python_bcs},
-          _pipes(pipes)
-    {
-        cross_section_area_inner_pipe = _pipes.inner_pipe.area();
-        cross_section_area_annulus =
-            _pipes.outer_pipe.area() - _pipes.inner_pipe.outsideArea();
-        cross_section_area_grout =
-            borehole_geometry.area() - _pipes.outer_pipe.outsideArea();
-
-        _thermal_resistances.fill(std::numeric_limits<double>::quiet_NaN());
-    }
+                     bool const use_python_bcs);
 
     static constexpr int number_of_unknowns = 3;
     static constexpr int number_of_grout_zones = 1;
diff --git a/ProcessLib/HeatTransportBHE/BHE/BHE_1U.cpp b/ProcessLib/HeatTransportBHE/BHE/BHE_1U.cpp
index a492ff5e01e238fcc3702b35d34692799f10042f..c5451e452d24fc05b09ff67d06eb6c7c5ea298af 100644
--- a/ProcessLib/HeatTransportBHE/BHE/BHE_1U.cpp
+++ b/ProcessLib/HeatTransportBHE/BHE/BHE_1U.cpp
@@ -21,6 +21,27 @@ namespace HeatTransportBHE
 {
 namespace BHE
 {
+BHE_1U::BHE_1U(BoreholeGeometry const& borehole,
+               RefrigerantProperties const& refrigerant,
+               GroutParameters const& grout,
+               FlowAndTemperatureControl const& flowAndTemperatureControl,
+               PipeConfigurationUType const& pipes,
+               bool const use_python_bcs)
+    : BHECommonUType{borehole, refrigerant,   grout, flowAndTemperatureControl,
+                     pipes,    use_python_bcs}
+{
+    _thermal_resistances.fill(std::numeric_limits<double>::quiet_NaN());
+
+    // Initialize thermal resistances.
+    auto values = visit(
+        [&](auto const& control) {
+            return control(refrigerant.reference_temperature,
+                           0. /* initial time */);
+        },
+        flowAndTemperatureControl);
+    updateHeatTransferCoefficients(values.flow_rate);
+}
+
 std::array<double, BHE_1U::number_of_unknowns> BHE_1U::pipeHeatCapacities()
     const
 {
@@ -207,6 +228,13 @@ std::array<double, BHE_1U::number_of_unknowns> BHE_1U::calcThermalResistances(
     // -------------------------------------------------------------------------
 }
 
+std::array<double, BHE_1U::number_of_unknowns> BHE_1U::crossSectionAreas() const
+{
+    return {{_pipes.inlet.area(), _pipes.outlet.area(),
+             borehole_geometry.area() / 2 - _pipes.inlet.area(),
+             borehole_geometry.area() / 2 - _pipes.outlet.area()}};
+}
+
 double BHE_1U::updateFlowRateAndTemperature(double const T_out,
                                             double const current_time)
 {
diff --git a/ProcessLib/HeatTransportBHE/BHE/BHE_1U.h b/ProcessLib/HeatTransportBHE/BHE/BHE_1U.h
index 1c310007b2059818867ddf7e48f82d666b90b8d0..6702d7322a67c87d4f5d750c758f8b5dd7fe28d9 100644
--- a/ProcessLib/HeatTransportBHE/BHE/BHE_1U.h
+++ b/ProcessLib/HeatTransportBHE/BHE/BHE_1U.h
@@ -46,22 +46,7 @@ public:
            GroutParameters const& grout,
            FlowAndTemperatureControl const& flowAndTemperatureControl,
            PipeConfigurationUType const& pipes,
-           bool const use_python_bcs)
-        : BHECommonUType{borehole, refrigerant,
-                         grout,    flowAndTemperatureControl,
-                         pipes,    use_python_bcs}
-    {
-        _thermal_resistances.fill(std::numeric_limits<double>::quiet_NaN());
-
-        // Initialize thermal resistances.
-        auto values = visit(
-            [&](auto const& control) {
-                return control(refrigerant.reference_temperature,
-                               0. /* initial time */);
-            },
-            flowAndTemperatureControl);
-        updateHeatTransferCoefficients(values.flow_rate);
-    }
+           bool const use_python_bcs);
 
     static constexpr int number_of_unknowns = 4;
     static constexpr int number_of_grout_zones = 2;
@@ -154,12 +139,7 @@ public:
         {0, 1}};
 
 public:
-    std::array<double, number_of_unknowns> crossSectionAreas() const
-    {
-        return {{_pipes.inlet.area(), _pipes.outlet.area(),
-                 borehole_geometry.area() / 2 - _pipes.inlet.area(),
-                 borehole_geometry.area() / 2 - _pipes.outlet.area()}};
-    }
+    std::array<double, number_of_unknowns> crossSectionAreas() const;
 
     void updateHeatTransferCoefficients(double const flow_rate);
 
diff --git a/ProcessLib/HeatTransportBHE/BHE/BHE_2U.cpp b/ProcessLib/HeatTransportBHE/BHE/BHE_2U.cpp
index ab30a1941d00e54d44f7b5e2b23a8b88a602af2b..734058abeabf52fd5e06cd5e07f5884f0ce69172 100644
--- a/ProcessLib/HeatTransportBHE/BHE/BHE_2U.cpp
+++ b/ProcessLib/HeatTransportBHE/BHE/BHE_2U.cpp
@@ -21,6 +21,27 @@ namespace HeatTransportBHE
 {
 namespace BHE
 {
+BHE_2U::BHE_2U(BoreholeGeometry const& borehole,
+               RefrigerantProperties const& refrigerant,
+               GroutParameters const& grout,
+               FlowAndTemperatureControl const& flowAndTemperatureControl,
+               PipeConfigurationUType const& pipes,
+               bool const use_python_bcs)
+    : BHECommonUType{borehole, refrigerant,   grout, flowAndTemperatureControl,
+                     pipes,    use_python_bcs}
+{
+    _thermal_resistances.fill(std::numeric_limits<double>::quiet_NaN());
+
+    // Initialize thermal resistances.
+    auto values = visit(
+        [&](auto const& control) {
+            return control(refrigerant.reference_temperature,
+                           0. /* initial time */);
+        },
+        flowAndTemperatureControl);
+    updateHeatTransferCoefficients(values.flow_rate);
+}
+
 std::array<double, BHE_2U::number_of_unknowns> BHE_2U::pipeHeatCapacities()
     const
 {
@@ -237,6 +258,20 @@ std::array<double, BHE_2U::number_of_unknowns> BHE_2U::calcThermalResistances(
     // -------------------------------------------------------------------------
 }
 
+std::array<double, BHE_2U::number_of_unknowns> BHE_2U::crossSectionAreas() const
+{
+    return {{
+        _pipes.inlet.area(),
+        _pipes.inlet.area(),
+        _pipes.outlet.area(),
+        _pipes.outlet.area(),
+        borehole_geometry.area() / 4 - _pipes.inlet.outsideArea(),
+        borehole_geometry.area() / 4 - _pipes.inlet.outsideArea(),
+        borehole_geometry.area() / 4 - _pipes.outlet.outsideArea(),
+        borehole_geometry.area() / 4 - _pipes.outlet.outsideArea(),
+    }};
+}
+
 double BHE_2U::updateFlowRateAndTemperature(double const T_out,
                                             double const current_time)
 {
diff --git a/ProcessLib/HeatTransportBHE/BHE/BHE_2U.h b/ProcessLib/HeatTransportBHE/BHE/BHE_2U.h
index 0eca5582707761334014dabd5a79ef9aafd59924..c63ff205d3833d960d6bdde40d045486b69604c4 100644
--- a/ProcessLib/HeatTransportBHE/BHE/BHE_2U.h
+++ b/ProcessLib/HeatTransportBHE/BHE/BHE_2U.h
@@ -46,22 +46,7 @@ public:
            GroutParameters const& grout,
            FlowAndTemperatureControl const& flowAndTemperatureControl,
            PipeConfigurationUType const& pipes,
-           bool const use_python_bcs)
-        : BHECommonUType{borehole, refrigerant,
-                         grout,    flowAndTemperatureControl,
-                         pipes,    use_python_bcs}
-    {
-        _thermal_resistances.fill(std::numeric_limits<double>::quiet_NaN());
-
-        // Initialize thermal resistances.
-        auto values = visit(
-            [&](auto const& control) {
-                return control(refrigerant.reference_temperature,
-                               0. /* initial time */);
-            },
-            flowAndTemperatureControl);
-        updateHeatTransferCoefficients(values.flow_rate);
-    }
+           bool const use_python_bcs);
 
     static constexpr int number_of_unknowns = 8;
     static constexpr int number_of_grout_zones = 4;
@@ -216,19 +201,7 @@ public:
         {0, 2}, {1, 3}};
 
 public:
-    std::array<double, number_of_unknowns> crossSectionAreas() const
-    {
-        return {{
-            _pipes.inlet.area(),
-            _pipes.inlet.area(),
-            _pipes.outlet.area(),
-            _pipes.outlet.area(),
-            borehole_geometry.area() / 4 - _pipes.inlet.outsideArea(),
-            borehole_geometry.area() / 4 - _pipes.inlet.outsideArea(),
-            borehole_geometry.area() / 4 - _pipes.outlet.outsideArea(),
-            borehole_geometry.area() / 4 - _pipes.outlet.outsideArea(),
-        }};
-    }
+    std::array<double, number_of_unknowns> crossSectionAreas() const;
 
     void updateHeatTransferCoefficients(double const flow_rate);
 
diff --git a/ProcessLib/HeatTransportBHE/HeatTransportBHEProcess.cpp b/ProcessLib/HeatTransportBHE/HeatTransportBHEProcess.cpp
index ac343e2f9ca8529018fdad5e880650ef8bd7006b..20d39b268bbc18693ab3d5403e586c6f25ea4a3c 100644
--- a/ProcessLib/HeatTransportBHE/HeatTransportBHEProcess.cpp
+++ b/ProcessLib/HeatTransportBHE/HeatTransportBHEProcess.cpp
@@ -283,28 +283,51 @@ void HeatTransportBHEProcess::createBHEBoundaryConditionTopBottom(
         // the BHE temperature is therefore bhe_i + 1
         const int variable_id = bhe_i + 1;
 
-        // Bottom and top nodes w.r.t. the z coordinate.
-        auto const bottom_top_nodes = std::minmax_element(
-            begin(bhe_nodes), end(bhe_nodes),
-            [&](auto const& a, auto const& b) {
-                return a->getCoords()[2] < b->getCoords()[2];
-            });
-        auto const bc_bottom_node_id = (*bottom_top_nodes.first)->getID();
-        auto const bc_top_node_id = (*bottom_top_nodes.second)->getID();
+        std::vector<MeshLib::Node*> bhe_boundary_nodes;
+        // cherry-pick the boundary nodes according to
+        // the number of connected line elements.
+        for (auto const& bhe_node : bhe_nodes)
+        {
+            // Count number of 1d elements connected with every BHE node.
+            const std::size_t n_line_elements = std::count_if(
+                bhe_node->getElements().begin(), bhe_node->getElements().end(),
+                [](MeshLib::Element const* elem) {
+                    return (elem->getDimension() == 1);
+                });
+
+            if (n_line_elements == 1)
+            {
+                bhe_boundary_nodes.push_back(bhe_node);
+            }
+        }
+
+        if (bhe_boundary_nodes.size() != 2)
+        {
+            OGS_FATAL(
+                "Error!!! The BHE boundary nodes are not correctly found, "
+                "for every single BHE, there should be 2 boundary nodes.");
+        }
+        auto get_global_index =
+            [&](std::size_t const node_id, int const component) {
+                return _local_to_global_index_map->getGlobalIndex(
+                    {_mesh.getID(), MeshLib::MeshItemType::Node, node_id},
+                    variable_id, component);
+            };
 
         auto get_global_bhe_bc_indices =
-            [&](std::size_t const node_id,
-                std::pair<int, int> const& in_out_component_id) {
+            [&](std::array<
+                std::pair<std::size_t /*node_id*/, int /*component*/>, 2>
+                    nodes_and_components) {
                 return std::make_pair(
-                    _local_to_global_index_map->getGlobalIndex(
-                        {_mesh.getID(), MeshLib::MeshItemType::Node, node_id},
-                        variable_id, in_out_component_id.first),
-                    _local_to_global_index_map->getGlobalIndex(
-                        {_mesh.getID(), MeshLib::MeshItemType::Node, node_id},
-                        variable_id, in_out_component_id.second));
+                    get_global_index(nodes_and_components[0].first,
+                                     nodes_and_components[0].second),
+                    get_global_index(nodes_and_components[1].first,
+                                     nodes_and_components[1].second));
             };
 
-        auto createBCs = [&](auto& bhe) {
+        auto createBCs = [&, bc_top_node_id = bhe_boundary_nodes[0]->getID(),
+                          bc_bottom_node_id =
+                              bhe_boundary_nodes[1]->getID()](auto& bhe) {
             for (auto const& in_out_component_id :
                  bhe.inflow_outflow_bc_component_ids)
             {
@@ -317,8 +340,11 @@ void HeatTransportBHEProcess::createBHEBoundaryConditionTopBottom(
                         // apply the customized top, inflow BC.
                         bcs.addBoundaryCondition(
                             ProcessLib::createBHEInflowPythonBoundaryCondition(
-                                get_global_bhe_bc_indices(bc_top_node_id,
-                                                          in_out_component_id),
+                                get_global_bhe_bc_indices(
+                                    {{{bc_top_node_id,
+                                       in_out_component_id.first},
+                                      {bc_top_node_id,
+                                       in_out_component_id.second}}}),
                                 bhe,
                                 *(_process_data.py_bc_object)));
                     }
@@ -339,8 +365,10 @@ void HeatTransportBHEProcess::createBHEBoundaryConditionTopBottom(
                     // Top, inflow, normal case
                     bcs.addBoundaryCondition(
                         createBHEInflowDirichletBoundaryCondition(
-                            get_global_bhe_bc_indices(bc_top_node_id,
-                                                      in_out_component_id),
+                            get_global_bhe_bc_indices(
+                                {{{bc_top_node_id, in_out_component_id.first},
+                                  {bc_top_node_id,
+                                   in_out_component_id.second}}}),
                             [&bhe](double const T, double const t) {
                                 return bhe.updateFlowRateAndTemperature(T, t);
                             }));
@@ -348,8 +376,10 @@ void HeatTransportBHEProcess::createBHEBoundaryConditionTopBottom(
                 // Bottom, outflow, all cases
                 bcs.addBoundaryCondition(
                     createBHEBottomDirichletBoundaryCondition(
-                        get_global_bhe_bc_indices(bc_bottom_node_id,
-                                                  in_out_component_id)));
+                        get_global_bhe_bc_indices(
+                            {{{bc_bottom_node_id, in_out_component_id.first},
+                              {bc_bottom_node_id,
+                               in_out_component_id.second}}})));
             }
         };
         visit(createBCs, _process_data._vec_BHE_property[bhe_i]);