diff --git a/ProcessLib/HeatTransportBHE/BHE/BHECommonCoaxial.cpp b/ProcessLib/HeatTransportBHE/BHE/BHECommonCoaxial.cpp
index 5f3b2cc6627a154d19f2bf7a35a6039d5cad1b56..f8a6d32eb6f46634ce3227c762ac496436562e08 100644
--- a/ProcessLib/HeatTransportBHE/BHE/BHECommonCoaxial.cpp
+++ b/ProcessLib/HeatTransportBHE/BHE/BHECommonCoaxial.cpp
@@ -87,7 +87,8 @@ BHECommonCoaxial::pipeHeatConductions() const
 }
 
 std::array<Eigen::Vector3d, BHECommonCoaxial::number_of_unknowns>
-BHECommonCoaxial::pipeAdvectionVectors() const
+    BHECommonCoaxial::pipeAdvectionVectors(
+        Eigen::Vector3d /*elem_direction_vec*/) const
 {
     double const rho_r = refrigerant.density;
     double const Cp_r = refrigerant.specific_heat_capacity;
diff --git a/ProcessLib/HeatTransportBHE/BHE/BHECommonCoaxial.h b/ProcessLib/HeatTransportBHE/BHE/BHECommonCoaxial.h
index b033891a34e223708f2c39ffeeab2b1b1f4aa31a..03264584ad6f8bf731e69f2dde65f73e51034cff 100644
--- a/ProcessLib/HeatTransportBHE/BHE/BHECommonCoaxial.h
+++ b/ProcessLib/HeatTransportBHE/BHE/BHECommonCoaxial.h
@@ -50,7 +50,8 @@ public:
 
     std::array<double, number_of_unknowns> pipeHeatConductions() const;
 
-    std::array<Eigen::Vector3d, number_of_unknowns> pipeAdvectionVectors()
+    std::array<Eigen::Vector3d, number_of_unknowns> pipeAdvectionVectors(
+        Eigen::Vector3d /*elem_direction_vec*/)
         const;
 
     double cross_section_area_inner_pipe, cross_section_area_annulus,
diff --git a/ProcessLib/HeatTransportBHE/BHE/BHE_1P.cpp b/ProcessLib/HeatTransportBHE/BHE/BHE_1P.cpp
index 9be0555074b296f4694a5e839c9ab988fe716071..59a8eb1abcde75c0a8008902293a392f7d9aec44 100644
--- a/ProcessLib/HeatTransportBHE/BHE/BHE_1P.cpp
+++ b/ProcessLib/HeatTransportBHE/BHE/BHE_1P.cpp
@@ -57,15 +57,17 @@ std::array<double, BHE_1P::number_of_unknowns> BHE_1P::pipeHeatConductions()
 }
 
 std::array<Eigen::Vector3d, BHE_1P::number_of_unknowns>
-BHE_1P::pipeAdvectionVectors() const
+BHE_1P::pipeAdvectionVectors(Eigen::Vector3d elem_direction_vec) const
 {
     double const& rho_r = refrigerant.density;
     double const& Cp_r = refrigerant.specific_heat_capacity;
+    Eigen::Vector3d adv_vector =
+        rho_r * Cp_r * _flow_velocity * elem_direction_vec;
 
-    return {{// pipe, Eq. 19
-             {0, 0, -rho_r * Cp_r * _flow_velocity},
-             // grout, Eq. 21
-             {0, 0, 0}}};
+    return {// pipe, Eq. 19
+            adv_vector,
+            // grout, Eq. 21
+            {0, 0, 0}};
 }
 
 double BHE_1P::compute_R_gs(double const chi, double const R_g)
diff --git a/ProcessLib/HeatTransportBHE/BHE/BHE_1P.h b/ProcessLib/HeatTransportBHE/BHE/BHE_1P.h
index 82f9b9e2b41676637e284b452edd1a9805698b49..d98916309b67d84233a41a349f250767f0547140 100644
--- a/ProcessLib/HeatTransportBHE/BHE/BHE_1P.h
+++ b/ProcessLib/HeatTransportBHE/BHE/BHE_1P.h
@@ -66,8 +66,8 @@ public:
 
     std::array<double, number_of_unknowns> pipeHeatConductions() const;
 
-    std::array<Eigen::Vector3d, number_of_unknowns> pipeAdvectionVectors()
-        const;
+    std::array<Eigen::Vector3d, number_of_unknowns> pipeAdvectionVectors(
+        Eigen::Vector3d elem_direction_vec) const;
 
     template <int NPoints,
               typename SingleUnknownMatrixType,
diff --git a/ProcessLib/HeatTransportBHE/BHE/BHE_1U.cpp b/ProcessLib/HeatTransportBHE/BHE/BHE_1U.cpp
index c5451e452d24fc05b09ff67d06eb6c7c5ea298af..915c9c20fb52f9dc0aa2f73a33f1313ee3c484fe 100644
--- a/ProcessLib/HeatTransportBHE/BHE/BHE_1U.cpp
+++ b/ProcessLib/HeatTransportBHE/BHE/BHE_1U.cpp
@@ -84,7 +84,7 @@ std::array<double, BHE_1U::number_of_unknowns> BHE_1U::pipeHeatConductions()
 }
 
 std::array<Eigen::Vector3d, BHE_1U::number_of_unknowns>
-BHE_1U::pipeAdvectionVectors() const
+    BHE_1U::pipeAdvectionVectors(Eigen::Vector3d /*elem_direction_vec*/) const
 {
     double const& rho_r = refrigerant.density;
     double const& Cp_r = refrigerant.specific_heat_capacity;
diff --git a/ProcessLib/HeatTransportBHE/BHE/BHE_1U.h b/ProcessLib/HeatTransportBHE/BHE/BHE_1U.h
index 6702d7322a67c87d4f5d750c758f8b5dd7fe28d9..5b3b1d1b315f7d5ebcf4d0b0b2138e7e1363ff32 100644
--- a/ProcessLib/HeatTransportBHE/BHE/BHE_1U.h
+++ b/ProcessLib/HeatTransportBHE/BHE/BHE_1U.h
@@ -55,7 +55,8 @@ public:
 
     std::array<double, number_of_unknowns> pipeHeatConductions() const;
 
-    std::array<Eigen::Vector3d, number_of_unknowns> pipeAdvectionVectors()
+    std::array<Eigen::Vector3d, number_of_unknowns> pipeAdvectionVectors(
+        Eigen::Vector3d /*elem_direction_vec*/)
         const;
 
     template <int NPoints,
diff --git a/ProcessLib/HeatTransportBHE/BHE/BHE_2U.cpp b/ProcessLib/HeatTransportBHE/BHE/BHE_2U.cpp
index 734058abeabf52fd5e06cd5e07f5884f0ce69172..ce09a6e0482981a4008e7a59da33239077c685fc 100644
--- a/ProcessLib/HeatTransportBHE/BHE/BHE_2U.cpp
+++ b/ProcessLib/HeatTransportBHE/BHE/BHE_2U.cpp
@@ -94,7 +94,7 @@ std::array<double, BHE_2U::number_of_unknowns> BHE_2U::pipeHeatConductions()
 }
 
 std::array<Eigen::Vector3d, BHE_2U::number_of_unknowns>
-BHE_2U::pipeAdvectionVectors() const
+    BHE_2U::pipeAdvectionVectors(Eigen::Vector3d /*elem_direction_vec*/) const
 {
     double const rho_r = refrigerant.density;
     double const Cp_r = refrigerant.specific_heat_capacity;
diff --git a/ProcessLib/HeatTransportBHE/BHE/BHE_2U.h b/ProcessLib/HeatTransportBHE/BHE/BHE_2U.h
index c63ff205d3833d960d6bdde40d045486b69604c4..d3c69e92d8f39ab3c5ecbef1127a6f39d44bce35 100644
--- a/ProcessLib/HeatTransportBHE/BHE/BHE_2U.h
+++ b/ProcessLib/HeatTransportBHE/BHE/BHE_2U.h
@@ -55,7 +55,8 @@ public:
 
     std::array<double, number_of_unknowns> pipeHeatConductions() const;
 
-    std::array<Eigen::Vector3d, number_of_unknowns> pipeAdvectionVectors()
+    std::array<Eigen::Vector3d, number_of_unknowns> pipeAdvectionVectors(
+        Eigen::Vector3d /*elem_direction_vec*/)
         const;
 
     template <int NPoints,
diff --git a/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerBHE-impl.h b/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerBHE-impl.h
index 33e35b1bd518f3f4330c3733e4d25d9df12ebb40..9583cd5d02bbd83f270673a5ffed2b0ed8c6193f 100644
--- a/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerBHE-impl.h
+++ b/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerBHE-impl.h
@@ -57,6 +57,15 @@ HeatTransportBHELocalAssemblerBHE<ShapeFunction, IntegrationMethod, BHEType>::
         _secondary_data.N[ip] = sm.N;
     }
 
+    // calculate the element direction vector
+    const double* coords_1 = e.getNode(0)->getCoords();
+    const double* coords_2 = e.getNode(1)->getCoords();
+    _element_direction_vector.setZero(3);
+    _element_direction_vector(0) = coords_2[0] - coords_1[0];
+    _element_direction_vector(1) = coords_2[1] - coords_1[1];
+    _element_direction_vector(2) = coords_2[2] - coords_1[2];
+    _element_direction_vector /= _element_direction_vector.norm();
+
     _R_matrix.setZero(bhe_unknowns_size, bhe_unknowns_size);
     _R_pi_s_matrix.setZero(bhe_unknowns_size, soil_temperature_size);
     _R_s_matrix.setZero(soil_temperature_size, soil_temperature_size);
@@ -65,12 +74,10 @@ HeatTransportBHELocalAssemblerBHE<ShapeFunction, IntegrationMethod, BHEType>::
     for (int idx_bhe_unknowns = 0; idx_bhe_unknowns < bhe_unknowns;
          idx_bhe_unknowns++)
     {
-        typename ShapeMatricesType::template MatrixType<
-            n_int_points, n_int_points>
+        typename ShapeMatricesType::template MatrixType<n_int_points,
+                                                        n_int_points>
             matBHE_loc_R = ShapeMatricesType::template MatrixType<
-                n_int_points,
-                n_int_points>::Zero(n_int_points,
-                                                n_int_points);
+                n_int_points, n_int_points>::Zero(n_int_points, n_int_points);
         // Loop over Gauss points
         for (unsigned ip = 0; ip < n_integration_points; ip++)
         {
@@ -131,7 +138,7 @@ void HeatTransportBHELocalAssemblerBHE<ShapeFunction, IntegrationMethod,
 
     auto const& pipe_heat_capacities = _bhe.pipeHeatCapacities();
     auto const& pipe_heat_conductions = _bhe.pipeHeatConductions();
-    auto const& pipe_advection_vectors = _bhe.pipeAdvectionVectors();
+    auto const& pipe_advection_vectors = _bhe.pipeAdvectionVectors(_element_direction_vector);
     auto const& cross_section_areas = _bhe.crossSectionAreas();
 
     // the mass and conductance matrix terms
@@ -155,26 +162,22 @@ void HeatTransportBHELocalAssemblerBHE<ShapeFunction, IntegrationMethod,
             auto const& A = cross_section_areas[idx_bhe_unknowns];
 
             int const single_bhe_unknowns_index =
-                bhe_unknowns_index +
-                n_int_points * idx_bhe_unknowns;
+                bhe_unknowns_index + n_int_points * idx_bhe_unknowns;
             // local M
             local_M
-                .template block<n_int_points,
-                                n_int_points>(
+                .template block<n_int_points, n_int_points>(
                     single_bhe_unknowns_index, single_bhe_unknowns_index)
                 .noalias() += N.transpose() * N * mass_coeff * A * w;
 
             // local K
             // laplace part
             local_K
-                .template block<n_int_points,
-                                n_int_points>(
+                .template block<n_int_points, n_int_points>(
                     single_bhe_unknowns_index, single_bhe_unknowns_index)
                 .noalias() += dNdx.transpose() * dNdx * lambda * A * w;
             // advection part
             local_K
-                .template block<n_int_points,
-                                n_int_points>(
+                .template block<n_int_points, n_int_points>(
                     single_bhe_unknowns_index, single_bhe_unknowns_index)
                 .noalias() +=
                 N.transpose() * advection_vector.transpose() * dNdx * A * w;
diff --git a/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerBHE.h b/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerBHE.h
index da02fb4465429a54951a904b4c22c0079ea8765f..d259d3884d2ec00c04b2675743b420fa8c4c6618 100644
--- a/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerBHE.h
+++ b/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerBHE.h
@@ -32,8 +32,7 @@ class HeatTransportBHELocalAssemblerBHE
     static constexpr int n_int_points = ShapeFunction::NPOINTS;
     static constexpr int soil_temperature_size = ShapeFunction::NPOINTS;
     static constexpr int soil_temperature_index = 0;
-    static constexpr int bhe_unknowns_size =
-        n_int_points * bhe_unknowns;
+    static constexpr int bhe_unknowns_size = n_int_points * bhe_unknowns;
     static constexpr int bhe_unknowns_index = ShapeFunction::NPOINTS;
     static constexpr int local_matrix_size =
         soil_temperature_size + bhe_unknowns_size;
@@ -91,6 +90,8 @@ private:
 
     SecondaryData<typename ShapeMatrices::ShapeType> _secondary_data;
 
+    Eigen::Vector3d _element_direction_vector;
+
     typename ShapeMatricesType::template MatrixType<bhe_unknowns_size,
                                                     bhe_unknowns_size>
         _R_matrix;