diff --git a/Tests/Data/Parabolic/T/3D_2U_BHE/3D_2U_BHE.prj b/Tests/Data/Parabolic/T/3D_2U_BHE/3D_2U_BHE.prj
index f12e3b2499f44e1ad0cfc20d0d17a4fd9adf7803..7136744786fdc70e9cca5d11f17f37d162b31fba 100644
--- a/Tests/Data/Parabolic/T/3D_2U_BHE/3D_2U_BHE.prj
+++ b/Tests/Data/Parabolic/T/3D_2U_BHE/3D_2U_BHE.prj
@@ -1,6 +1,9 @@
<?xml version="1.0" encoding="ISO-8859-1"?>
<OpenGeoSysProject>
<mesh>3D_2U_BHE.vtu</mesh>
+ <!-- use the following mesh for benchmark with FEFLOW model
+ <mesh>3D_2U_BHE_benchmark.vtu</mesh>
+ -->
<geometry>3D_2U_BHE.gml</geometry>
<processes>
<process>
diff --git a/Tests/Data/Parabolic/T/3D_2U_BHE/3D_2U_BHE_powerBC.prj b/Tests/Data/Parabolic/T/3D_2U_BHE/3D_2U_BHE_powerBC.prj
new file mode 100644
index 0000000000000000000000000000000000000000..78c0d7764d6a8057049ce2858f4bdb91806b24c5
--- /dev/null
+++ b/Tests/Data/Parabolic/T/3D_2U_BHE/3D_2U_BHE_powerBC.prj
@@ -0,0 +1,231 @@
+<?xml version="1.0" encoding="ISO-8859-1"?>
+<OpenGeoSysProject>
+ <mesh>3D_2U_BHE_benchmark.vtu</mesh>
+ <geometry>3D_2U_BHE.gml</geometry>
+ <processes>
+ <process>
+ <name>HeatTransportBHE</name>
+ <type>HEAT_TRANSPORT_BHE</type>
+ <integration_order>2</integration_order>
+ <process_variables>
+ <process_variable>temperature_soil</process_variable>
+ <process_variable>temperature_BHE1</process_variable>
+ </process_variables>
+ <borehole_heat_exchangers>
+ <borehole_heat_exchanger>
+ <type>2U</type>
+ <flow_and_temperature_control>
+ <!-- if power type boundary conditions are selected, the related
+ power value should be set with half of the presumed BHE power-->
+ <type>PowerCurveConstantFlow</type>
+ <power_curve>heating_load</power_curve>
+ <flow_rate>2.0e-4</flow_rate>
+ </flow_and_temperature_control>
+ <borehole>
+ <length>18.0</length>
+ <diameter>0.13</diameter>
+ </borehole>
+ <grout>
+ <density>2190.0</density>
+ <porosity>0.0</porosity>
+ <specific_heat_capacity>1141.55</specific_heat_capacity>
+ <thermal_conductivity>1</thermal_conductivity>
+ </grout>
+ <pipes>
+ <inlet>
+ <diameter> 0.0378</diameter>
+ <wall_thickness>0.0029</wall_thickness>
+ <wall_thermal_conductivity>0.42</wall_thermal_conductivity>
+ </inlet>
+ <outlet>
+ <diameter>0.0378</diameter>
+ <wall_thickness>0.0029</wall_thickness>
+ <wall_thermal_conductivity>0.42</wall_thermal_conductivity>
+ </outlet>
+ <distance_between_pipes>0.053</distance_between_pipes>
+ <longitudinal_dispersion_length>0.001</longitudinal_dispersion_length>
+ </pipes>
+ <refrigerant>
+ <density>1052</density>
+ <viscosity>0.0003</viscosity>
+ <specific_heat_capacity>3802.28</specific_heat_capacity>
+ <thermal_conductivity>0.48</thermal_conductivity>
+ <reference_temperature>22</reference_temperature>
+ </refrigerant>
+ </borehole_heat_exchanger>
+ </borehole_heat_exchangers>
+ </process>
+ </processes>
+ <media>
+ <medium id="0">
+ <phases>
+ <phase>
+ <type>AqueousLiquid</type>
+ <properties>
+ <property>
+ <name>phase_velocity</name>
+ <type>Constant</type>
+ <value>0 0 0</value>
+ </property>
+ <property>
+ <name>specific_heat_capacity</name>
+ <type>Constant</type>
+ <value>4068</value>
+ </property>
+ <property>
+ <name>density</name>
+ <type>Constant</type>
+ <value>992.92</value>
+ </property>
+ </properties>
+ </phase>
+ <phase>
+ <type>Solid</type>
+ <properties>
+ <property>
+ <name>specific_heat_capacity</name>
+ <type>Constant</type>
+ <value>1778</value>
+ </property>
+ <property>
+ <name>density</name>
+ <type>Constant</type>
+ <value>1800</value>
+ </property>
+ </properties>
+ </phase>
+ <phase>
+ <type>Gas</type>
+ <properties>
+ <property>
+ <name>specific_heat_capacity</name>
+ <type>Constant</type>
+ <value>1000</value>
+ </property>
+ <property>
+ <name>density</name>
+ <type>Constant</type>
+ <value>2500</value>
+ </property>
+ </properties>
+ </phase>
+ </phases>
+ <properties>
+ <property>
+ <name>porosity</name>
+ <type>Constant</type>
+ <value>0</value>
+ </property>
+ <property>
+ <name>thermal_conductivity</name>
+ <type>Constant</type>
+ <value>2.78018</value>
+ </property>
+ <property>
+ <name>thermal_longitudinal_dispersivity</name>
+ <type>Constant</type>
+ <value>0</value>
+ </property>
+ <property>
+ <name>thermal_transversal_dispersivity</name>
+ <type>Constant</type>
+ <value>0</value>
+ </property>
+ </properties>
+ </medium>
+ </media>
+ <time_loop>
+ <processes>
+ <process ref="HeatTransportBHE">
+ <nonlinear_solver>basic_picard</nonlinear_solver>
+ <convergence_criterion>
+ <type>DeltaX</type>
+ <norm_type>NORM2</norm_type>
+ <reltol>1e-6</reltol>
+ </convergence_criterion>
+ <time_discretization>
+ <type>BackwardEuler</type>
+ </time_discretization>
+ <time_stepping>
+ <type>FixedTimeStepping</type>
+ <t_initial> 0.0 </t_initial>
+ <t_end> 3600 </t_end>
+ <timesteps>
+ <pair><repeat>60</repeat><delta_t>60</delta_t></pair>
+ </timesteps>
+ </time_stepping>
+ </process>
+ </processes>
+ <output>
+ <type>VTK</type>
+ <prefix>3D_2U_BHE_powerBC</prefix>
+ <timesteps>
+ <pair><repeat> 1</repeat><each_steps> 1 </each_steps></pair>
+ </timesteps>
+ <variables>
+ <variable>temperature_soil</variable>
+ <variable>temperature_BHE1</variable>
+ </variables>
+ <suffix>_ts_{:timestep}_t_{:time}</suffix>
+ </output>
+ </time_loop>
+ <parameters>
+ <parameter>
+ <name>T0</name>
+ <type>Constant</type>
+ <value>22</value>
+ </parameter>
+ <parameter>
+ <name>T0_BHE1</name>
+ <type>Constant</type>
+ <values>22 22 22 22 22 22 22 22</values>
+ </parameter>
+ </parameters>
+ <process_variables>
+ <process_variable>
+ <name>temperature_soil</name>
+ <components>1</components>
+ <order>1</order>
+ <initial_condition>T0</initial_condition>
+ <boundary_conditions>
+ </boundary_conditions>
+ </process_variable>
+ <process_variable>
+ <name>temperature_BHE1</name>
+ <components>8</components>
+ <order>1</order>
+ <initial_condition>T0_BHE1</initial_condition>
+ </process_variable>
+ </process_variables>
+ <nonlinear_solvers>
+ <nonlinear_solver>
+ <name>basic_picard</name>
+ <type>Picard</type>
+ <max_iter>100</max_iter>
+ <linear_solver>general_linear_solver</linear_solver>
+ </nonlinear_solver>
+ </nonlinear_solvers>
+ <linear_solvers>
+ <linear_solver>
+ <name>general_linear_solver</name>
+ <lis>-i cg -p jacobi -tol 1e-16 -maxiter 10000</lis>
+ <eigen>
+ <solver_type>BiCGSTAB</solver_type>
+ <precon_type>ILUT</precon_type>
+ <max_iteration_step>1000</max_iteration_step>
+ <error_tolerance>1e-16</error_tolerance>
+ </eigen>
+ <petsc>
+ <prefix>gw</prefix>
+ <parameters>-gw_ksp_type cg -gw_pc_type bjacobi -gw_ksp_rtol 1e-16 -gw_ksp_max_it 10000</parameters>
+ </petsc>
+ </linear_solver>
+ </linear_solvers>
+ <curves>
+ <curve>
+ <name>heating_load</name>
+ <coords>0 3600</coords>
+ <values>-157.86 -157.86</values>
+ </curve>
+ </curves>
+</OpenGeoSysProject>
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 c76803e626bddd9ce1a930c9b780396971b67983..738d656b3f54b84d6f5e31e445f2841fe68e0708 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
@@ -1,5 +1,5 @@
+++
-author = "Chaofan Chen, Haibing Shao"
+author = "Chaofan Chen, Shuang Chen, Haibing Shao"
date = "2018-09-25T13:44:00+01:00"
title = "3D 2U BHE"
weight = 123
@@ -32,20 +32,54 @@ The numerical model was established based on dual continuum method developed by
## OGS Input Files
+For this benchmark, Two different scenarios were carried out by applying two different boundary conditions imposed on the BHE.
+
+* Fixed inflow boundary condition
+
The detailed input parameters can be seen from the 3D_2U_BHE.prj file. The inflow temperature of the BHE, which was imposed as boundary condition of the BHE is shown in Figure 1. All the initial temperatures are set as 22 $^{\circ}$C. The flow rate within each U-pipe is set to $2.0\times10^{-4}$ $\mathrm{m^{3} s^{-1}}$ during the whole simulation time.
{{< img src="../3D_2U_BHE_figures/In_out_temperature_comparison.png" width="200">}}
Figure 1: Inflow temperature curve and outflow temperature comparison
+* Fixed power boundary condition
+
+The detailed input parameters can be seen from the 3D_2U_BHE_powerBC.prj file.
+A -315.72 $W$ thermal load is imposed on the BHE to extract the heat from the subsurface during the entire simulation.
+All the other parameters adopted in the model is same as the ones used in the scenario with fixed inflow boundary condition.
+
+## FEFLOW Input Files
+For the benchmark a FEFLOW model is presented.
+The mesh used in the OGS model is directly coverted from the FEFLOW model mesh, to ensure that there is no influence to the comparison results from the mesh side.
+Both the FEFLOW and ogs model mesh can be found in the ogs GitLab (<https://gitlab.opengeosys.org/ogs/ogs/-/merge_requests/3426>).
+
+## FEFLOW Input Files
+
+
## Results
-The OGS numerical outflow temperature over time was compared against results of the FEFLOW software as shown in the Figure 1. And the vertical distributed temperature of circulating water was presented in Figure 2 after operation for 3000 s. The comparison figures demonstrate that the OGS numerical results and FEFLOW results can match very well and the biggest absolute error of outflow temperature is 0.19 $^{\circ}$C at the starting up stage, while such error decreases to 0.05 $^{\circ}$C after 3600 s' operation. The maximum relative error of vertical temperature is 0.015 \% after operation for 3000 s.
+The computed resutls from scenario by adopting the fixed inflow boundary condition are illustracted in Figure 1 and Figure 2.
+The OGS numerical outflow temperature over time was compared against results of the FEFLOW software as shown in the Figure 1. And the vertical distributed temperature of circulating water was presented in Figure 2 after operation for 3300 s.
+The comparison figures demonstrate that the OGS numerical results and FEFLOW results can match very well and the biggest absolute error of outflow temperature is 0.20 $^{\circ}$C after 360 s' operation, while such error decreases to 0.037 $^{\circ}$C after 3600 s' operation. The maximum relative error of vertical temperature is 0.019 \% after operation for 3300 s.
{{< img src="../3D_2U_BHE_figures/vertical_temperature_distribution.png" width="200">}}
-Figure 2: Comparison of vertical temperature distribution
+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.
+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.
+
+{{< 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
## References
[1] Diersch, H. J., Bauer, D., Heidemann, W., Rühaak, W., & Schätzl, P. (2011). Finite element modeling of borehole heat exchanger systems: Part 1. Fundamentals. Computers & Geosciences, 37(8), 1122-1135.
+
+[2] FEFLOW online documentation. URL: <http://www.feflow.info/html/help72/feflow/14_References/GUI/Dialogs/bhe_editor.html>.
\ No newline at end of file
diff --git a/web/content/docs/benchmarks/heat-transport-bhe/3D_2U_BHE_figures/In_out_temperature_comparison.png b/web/content/docs/benchmarks/heat-transport-bhe/3D_2U_BHE_figures/In_out_temperature_comparison.png
index 062633fb8b05d3090bfc614aaef7d7e9f9588d61..1e33ca2a2ba571d1517a0509da3668285142cd44 100644
Binary files a/web/content/docs/benchmarks/heat-transport-bhe/3D_2U_BHE_figures/In_out_temperature_comparison.png and b/web/content/docs/benchmarks/heat-transport-bhe/3D_2U_BHE_figures/In_out_temperature_comparison.png differ
diff --git a/web/content/docs/benchmarks/heat-transport-bhe/3D_2U_BHE_figures/vertical_temperature_distribution.png b/web/content/docs/benchmarks/heat-transport-bhe/3D_2U_BHE_figures/vertical_temperature_distribution.png
index 7132ea037d3f8d9826119db8dfa1f4c5d3de2baa..a14cf42a3382f74c0b21a032600fed4ebca7f3ee 100644
Binary files a/web/content/docs/benchmarks/heat-transport-bhe/3D_2U_BHE_figures/vertical_temperature_distribution.png and b/web/content/docs/benchmarks/heat-transport-bhe/3D_2U_BHE_figures/vertical_temperature_distribution.png differ
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
new file mode 100644
index 0000000000000000000000000000000000000000..fb365219aad9fd0345e720be395df8694a3e5052
Binary files /dev/null and b/web/content/docs/benchmarks/heat-transport-bhe/3D_2U_BHE_figures/vertical_temperature_distribution_powerBC.png differ