diff --git a/ProcessLib/HeatTransportBHE/BHE/BHE_1U.cpp b/ProcessLib/HeatTransportBHE/BHE/BHE_1U.cpp
index 4a6272803f63b23b09dc726e9596346d7645cf1f..2ef82cbe6b649f3039adc7a35e9d1c90d51b5ab8 100644
--- a/ProcessLib/HeatTransportBHE/BHE/BHE_1U.cpp
+++ b/ProcessLib/HeatTransportBHE/BHE/BHE_1U.cpp
@@ -124,8 +124,9 @@ double compute_R_gg(double const chi, double const R_gs, double const R_ar,
 /// Check if constraints regarding negative thermal resistances are violated
 /// apply correction procedure.
 /// Section (1.5.5) in FEFLOW White Papers Vol V.
-std::tuple<double, double, double> thermalResistancesGroutSoil(
-    double const chi, double const R_ar, double const R_g)
+std::array<double, 3> thermalResistancesGroutSoil(double const chi,
+                                                  double const R_ar,
+                                                  double const R_g)
 {
     double R_gs = compute_R_gs(chi, R_g);
     double R_gg =
diff --git a/ProcessLib/HeatTransportBHE/BHE/BHE_2U.cpp b/ProcessLib/HeatTransportBHE/BHE/BHE_2U.cpp
index 0fbc897fb0ee9b151415319da2335589c3993317..99555a2ea541cb02815d5efbc62664e8c1db6ba5 100644
--- a/ProcessLib/HeatTransportBHE/BHE/BHE_2U.cpp
+++ b/ProcessLib/HeatTransportBHE/BHE/BHE_2U.cpp
@@ -142,11 +142,10 @@ double compute_R_gg_2U(double const chi, double const R_gs, double const R_ar,
 /// Check if constraints regarding negative thermal resistances are violated
 /// apply correction procedure.
 /// Section (1.5.5) in FEFLOW White Papers Vol V.
-std::tuple<double, double, double, double> thermalResistancesGroutSoil2U(
-    double const chi,
-    double const R_ar_1,
-    double const R_ar_2,
-    double const R_g)
+std::array<double, 4> thermalResistancesGroutSoil2U(double const chi,
+                                                    double const R_ar_1,
+                                                    double const R_ar_2,
+                                                    double const R_g)
 {
     double R_gs = compute_R_gs_2U(chi, R_g);
     double R_gg_1 = compute_R_gg_2U(chi, R_gs, R_ar_1, R_g);
diff --git a/web/content/docs/benchmarks/heat-transport-bhe/3D_2U_BHE.md b/web/content/docs/benchmarks/heat-transport-bhe/3D_2U_BHE.md
index 738d656b3f54b84d6f5e31e445f2841fe68e0708..b806d437db0bcc0026bfe575cf9c0aa98e8bdf0a 100644
--- a/web/content/docs/benchmarks/heat-transport-bhe/3D_2U_BHE.md
+++ b/web/content/docs/benchmarks/heat-transport-bhe/3D_2U_BHE.md
@@ -66,17 +66,18 @@ The comparison figures demonstrate that the OGS numerical results and FEFLOW res
 
 Figure 2: Comparison of vertical temperature distribution from scenario by adopting the fixed inflow boundary condition
 
-Figure 3 shows the the vertical distributed temperature of circulating fluid after operation for 3300 s by adopting the power boundary condition in OGS and FEFLOW models.
-A 0.18 $^{\circ}$C difference is found between the averaged vertical temperature from the two models.
+Figure 3 shows the vertical distributed temperature of circulating fluid after operation for 3300 s by adopting different power boundary conditions in OGS and FEFLOW models.
+Compared to the resutls computed from the OGS model with using a fixed power boundary condtion (illustrated as the blue and green line), A 0.18 $^{\circ}$C difference is found for the averaged vertical temperature from the FEFLOW model (illustrated as the dotted line).
 The reason to the results difference is due to the different power boundary condition type adopted in the two software.
 In FEFLOW the power boundary condtion is based on the outlet temperature calculated from the last time step (non-iterative).
 Compared to it, the default power boundary condition adopted in the OGS `Heat_Transport_BHE` process is based on the outlet temperature calculated from the current time step (with-iterative).
 Besides, by setting python bindings, the current OGS `Heat_Transport_BHE` process is capable to adopt the power boundary condition type used in the FEFLOW software.
-In this way, the computed vertical distributed circulating fluid temperature difference between the OGS and FEFLOW models is becoming much closer to each other.
+In this way, the computed vertical distributed circulating fluid temperature is updated to the black and red solid line illustrated in the figure 3.
+It shows that in this case, the difference between the OGS and FEFLOW models is becoming much closer to each other, which is about 0.037 $^{\circ}$C.
 
 {{< img src="../3D_2U_BHE_figures/vertical_temperature_distribution_powerBC.png" width="200">}}
 
-Figure 3: Comparison of vertical temperature distribution from scenario by adopting the power boundary condition
+Figure 3: Comparison of vertical temperature distribution from scenarios by adopting the power boundary conditions
 
 ## References
 
diff --git a/web/content/docs/benchmarks/heat-transport-bhe/3D_2U_BHE_figures/vertical_temperature_distribution_powerBC.png b/web/content/docs/benchmarks/heat-transport-bhe/3D_2U_BHE_figures/vertical_temperature_distribution_powerBC.png
index fb365219aad9fd0345e720be395df8694a3e5052..1d1f2bc34b85fec755072784a6c932da26e8dd8b 100644
Binary files a/web/content/docs/benchmarks/heat-transport-bhe/3D_2U_BHE_figures/vertical_temperature_distribution_powerBC.png and b/web/content/docs/benchmarks/heat-transport-bhe/3D_2U_BHE_figures/vertical_temperature_distribution_powerBC.png differ