diff --git a/ProcessLib/HeatTransportBHE/Tests.cmake b/ProcessLib/HeatTransportBHE/Tests.cmake
index 4434c63ec81..a2c60ddefd0 100644
--- a/ProcessLib/HeatTransportBHE/Tests.cmake
+++ b/ProcessLib/HeatTransportBHE/Tests.cmake
@@ -160,3 +160,7 @@ AddTest(
     BHE_1P_newton_ts_10_t_600.000000.vtu BHE_1P_newton_ts_10_t_600.000000.vtu temperature_BHE1 temperature_BHE1 1e-12 1e-14
     BHE_1P_newton_ts_10_t_600.000000.vtu BHE_1P_newton_ts_10_t_600.000000.vtu temperature_soil temperature_soil 1e-12 1e-13
 )
+
+if(NOT OGS_USE_PETSC)
+    NotebookTest(NOTEBOOKFILE Parabolic/T/BHE_1P/pipe_flow_ebhe.md RUNTIME 600)
+endif()
diff --git a/Tests/Data/Parabolic/T/BHE_1P/BHE_1P.prj b/Tests/Data/Parabolic/T/BHE_1P/BHE_1P.prj
index 5f0776157f8..c0a197093dd 100644
--- a/Tests/Data/Parabolic/T/BHE_1P/BHE_1P.prj
+++ b/Tests/Data/Parabolic/T/BHE_1P/BHE_1P.prj
@@ -144,9 +144,25 @@
                     <t_end> 600 </t_end>
                     <timesteps>
                         <pair>
-                            <repeat>10</repeat>
+                            <repeat>60</repeat>
                             <delta_t>60</delta_t>
                         </pair>
+                        <pair>
+                            <repeat>60</repeat>
+                            <delta_t>360</delta_t>
+                        </pair>
+                        <pair>
+                            <repeat>41</repeat>
+                            <delta_t>3600</delta_t>
+                        </pair>
+                        <pair>
+                            <repeat>16</repeat>
+                            <delta_t>43200</delta_t>
+                        </pair>
+                        <pair>
+                            <repeat>28</repeat>
+                            <delta_t>86400</delta_t>
+                        </pair>
                     </timesteps>
                 </time_stepping>
             </process>
diff --git a/web/content/docs/benchmarks/heat-transport-bhe/pipe_flow_EBHE/pipe_flow_3d_model.png b/Tests/Data/Parabolic/T/BHE_1P/figures/pipe_flow_3d_model.png
similarity index 100%
rename from web/content/docs/benchmarks/heat-transport-bhe/pipe_flow_EBHE/pipe_flow_3d_model.png
rename to Tests/Data/Parabolic/T/BHE_1P/figures/pipe_flow_3d_model.png
diff --git a/Tests/Data/Parabolic/T/BHE_1P/pipe_flow_ebhe.md b/Tests/Data/Parabolic/T/BHE_1P/pipe_flow_ebhe.md
new file mode 100644
index 00000000000..bfd9bba7c02
--- /dev/null
+++ b/Tests/Data/Parabolic/T/BHE_1P/pipe_flow_ebhe.md
@@ -0,0 +1,425 @@
++++
+title = "Wellbore Heat Transport - EUBHE"
+date = "2023-08-03"
+author = "Chaofan Chen, Haibing Shao, Frieder Loer"
+image = "figures/pipe_flow_3d_model.png"
+web_subsection = "heat-transport-bhe"
+project = ["Parabolic/T/BHE_1P/BHE_1P.prj"]
++++
+
+
+## Problem description
+
+(Ramey et al. (1962)) proposed the analytical solution concerning the wellbore heat transmission, which can be used to quantify the fluid temperature change in the wellbore.  
+In order to verify the single pipe flow model in the OGS, the numerical results were compared with Ramey's analytical solution (Ramey et al. (1962)).  
+The detailed calculation of the Ramey's analytical solution is given below.  
+A detailed analysis of an enhanced U-tube borehole heat exchanger (EUBHE) can be found in (Chen, C. et al. (2021)).
+
+## Model Setup
+
+In this benchmark, the length of the wellbore is 30 m as shown in Figure 1 and the cold water is injected into the inlet point of the wellbore with the constant temperature of 20°C.  
+The initial temperature of the fluid and grout in the wellbore is 20 $^{\circ}$C, and temperature of the surrounding rock is 55 $^{\circ}$C.  
+The wellbore and pipe diameter are 0.28 m and 0.25826 m, respectively.  
+And the flow rate is 0.0002 $m^3/s$.
+
+![Single pipe flow model](./figures/pipe_flow_3d_model.png)
+
+Figure 1: Single pipe flow model
+
+```python
+# Import required libraries
+
+import numpy as np
+import matplotlib.pyplot as plt
+import pandas as pd
+```
+
+## Ramey's analytical solution
+
+The equations are followed by their implementation in python.
+
+```python
+# Define constants
+
+T_d = 55 # Temperature of undistributed formation
+T_i = 20 # Temperature of inlet fluid
+Z = np.insert(np.linspace(0.3,30,100), 0, 1e-10) # Depth below surface, start with very small value to prevent division by zero
+q = 0.0002 # Flow rate of fluid in wellbore, m3/s
+rho_f = 1000 # Density of fluid, unit kg/m3
+c_p_f = 4190 # Specific heat capacity of fluid in wellbore
+mu_f = 1.14e-3  # Unit kg/m/s;
+lambda_f = 0.59  # Unit W/m/K
+rho_re = 1800 # Density of geothermal reservoir
+c_p_re = 1778 # Specific heat capacity of geothermal reservoir
+lambda_re = 2.78018 # Heat conductivity coefficient of geothermal reservoir
+lambda_g = 0.73 # Heat conductivity coefficient of grout and pipe
+lambda_pi = 1.3
+r_pi = 0.12913 # Inner radius of pipe and wellbore
+r_b = 0.14
+t_pi = 0.00587 # Thickness of the pipe
+t = 86400 * 30 # Operation time
+```
+
+In Ramey's analytical solution (Ramey et al. (1962)), the outlet temperature of the pipe inside the wellbore can be calculated by
+
+$$
+    T_o(t) = T_{s} + (T_i(t) - T_{s})\exp(-\Delta z/X)
+$$
+
+```python
+def temp(T_d, T_i, delta_z, X):
+    return T_d + (T_i - T_d) * np.exp(-delta_z / X)
+```
+
+The coefficient $X$ is determined by
+
+$$
+    X = \frac{q\rho_fc_{p,f}(\lambda_{s}+r_pUf(t))}{2\pi r_pU \lambda_{s}}
+$$
+
+where $q$ is the flow rate of the fluid in the wellbore.
+
+```python
+def coefficient_x(q, rho_f, c_p_f, lambda_re, r_pi, U, f_t):
+    return (q * rho_f * c_p_f) * (lambda_re + r_pi * U * f_t) / (2 * np.pi * r_pi * U * lambda_re)
+```
+
+With dimensionless time
+
+$$
+    t_D = \frac{\lambda_{s}t}{(\rho_{s}c_{p,s}r_b^{2})}
+$$
+
+```python
+def dimless_t(lambda_re, delta_t, rho_re, c_p_re, r_b):
+    return lambda_re * delta_t / (rho_re * c_p_re * r_b * r_b)
+```
+
+the time function $f(t)$ can be calculated by
+
+<!-- markdownlint-disable reference-links-images -->
+
+$$
+    f(t) = [0.4063+0.5\ln(t_D)][1+\frac{0.6}{t_D}], t_D > 1.5
+$$
+
+<!-- markdownlint-enable reference-links-images -->
+
+$$
+    f(t) = 1.1281\sqrt{t_D}(1-0.3\sqrt{t_D}), t_D \leqslant 1.5
+$$
+
+```python
+def time_function(t_D):
+    if t_D > 1.5:
+                f_t = (0.4063 + 0.5 * np.log(t_D)) * (1 + 0.6 / t_D)
+    else:
+        f_t = 1.1281 * np.sqrt(t_D) * (1 - 0.3 * np.sqrt(t_D))
+    
+    return f_t
+
+```
+
+The Prandtl and Reynolds number can be calculated as follows
+
+$$
+    \mathrm{Pr} = \frac{\mu_f c_{p,f}}{\lambda_f}
+$$
+
+$$
+    \mathrm{Re} = \frac{\rho_f v d_{pi}}{\mu_f}
+$$
+
+where, $\mu_f, \rho_f$ and $\lambda_f$ is the fluid viscosity, density and thermal conductivity.
+
+```python
+# Calculate the convective film resistance inside pipe
+v = q / (np.pi * r_pi * r_pi)
+
+# Reynolds number
+Pr = mu_f * c_p_f / lambda_f
+Re = rho_f * v * (2 * r_pi) / mu_f
+
+```
+
+The Nusselt number can be determined by the following equation (Diersch, (2011)):
+
+$$
+    \mathrm{Nu} = 4.364,\  \mathrm{Re} < 2300
+$$
+
+$$
+    \mathrm{Nu} = \frac{(\xi_{k}/8)\  \mathrm{Re}_{k}\  \mathrm{Pr}}{1+12.7\sqrt{{\xi_k}/8}(\mathrm{Pr}^{2/3}-1)} [ 1+(\frac{{d_k}^{i}}{L})^{2/3}], Re \geq 10^4
+$$
+
+$$
+
+$$
+    \mathrm{Nu} = (1-\gamma_{k})\  4.364 + \gamma_{k} ( \frac{(0.0308/8)10^{4}\mathrm{Pr}}{1+12.7\  \sqrt{0.0308/8}(\mathrm{Pr}^{2/3}-1)} \left[ 1+\left(\frac{d_{k}^{i}}{L}\right)^{2/3} \right] ), 2 300 \leq \mathrm{Re} < 10^{4}
+$$
+
+with
+
+$$
+    \xi_{k}=(1.8\  \log_{10}\mathrm{Re}_{k}-1.5)^{-2}
+$$
+
+and
+
+$$
+    \gamma_{k} = \frac{\mathrm{Re}_{k}-2300}{10^{4}-2300}
+$$
+
+```python
+# Estimate convective film coefficient
+xi = (1.8 * np.log(Re) - 1.5)**(-2)
+Nu_l = (xi / 8.0 * Re * Pr) / (1.0 + 12.7 * np.sqrt(xi / 8.0) * ((Pr)**(2.0 / 3.0) - 1.0)) * (1.0 + (r_pi / Z)**(2.0 / 3.0))
+Nu_lam = 4.364
+gamma = (Re - 2300) / (10000 - 2300)
+
+Nu_m = (1.0 - gamma) * 4.364 + gamma * ((0.0308 / 8.0 * 1.0e4 * Pr) / (1.0 + 12.7 * np.sqrt(0.0308 / 8.0) * ((Pr)**(2.0 / 3.0) - 1.0)) * (1.0 + (r_pi / Z)**(2.0 / 3.0)))
+
+# Evaluate Nusselt number; based on the value of Re choose appropriate correlation
+if Re < 2300:
+    Nu_p = Nu_lam
+elif Re > 10000:
+    Nu_p = Nu_l
+else:
+    Nu_p = Nu_m
+```
+
+and the overall heat transfer coefficient $U$ is written as follows,
+
+$$
+    U = [\frac{r_{pi}+t_{pi}}{r_{pi}h}+(r_{pi}+t_{pi})(\frac{\ln{\frac{r_{pi}+t_{pi}}{r_{pi}}}}{\lambda_{pi}}+\frac{\ln{\frac{r_b}{r_{pi}+t_{pi}}}}{\lambda_{grout}})]^{-1}
+$$
+
+with
+
+$$
+    h = \frac{\lambda_f Nu}{2r_{pi}}
+$$
+
+```python
+h = lambda_f * Nu_p / (2 * r_pi)  # Unit: W/m2/K
+
+U = 1 / (((r_pi + t_pi) / (r_pi * h)) + (r_pi + t_pi) * (np.log((r_pi + t_pi) / r_pi) / lambda_pi + np.log(r_b / (r_pi + t_pi)) / lambda_g))
+```
+
+The friction factor $f$, is evaluated by Churchill correlation (Churchill et al. (1977)),
+
+$$
+    f = \frac{1}{(\frac{1}{[((\frac{8}{Re})^{10}+(\frac{Re}{36500})^{20})]^{1/2}}+[2.21(\ln{\frac{Re}{7}})]^{10})^{1/5}}
+$$
+
+```python
+# Churchill correlation for friction factor
+f = 1 / ((1 / (np.sqrt((8 / Re) ** 10 + (Re / 36500) ** 20))) + (2.21 * np.log(Re / 7)) ** 10) ** (1 / 5)
+```
+
+With these equations the analytical solution can be computed:
+
+```python
+# Compute solution
+
+time_range = range(0, t+1, 43200)
+
+# Create a dataframe to store results for different z- and t-steps
+row_names = Z
+column_names = time_range
+
+
+df = pd.DataFrame(index=row_names, columns=column_names) 
+
+
+for delta_z in Z:
+    T_0 = []
+    for delta_t in time_range:
+        
+        # Compute dimensionless time
+        t_D = dimless_t(lambda_re, delta_t, rho_re, c_p_re, r_b)
+
+        # Compute time function
+        f_t = time_function(t_D)
+
+        # Compute X
+        X = coefficient_x(q, rho_f, c_p_f, lambda_re, r_pi, U, f_t)
+
+        # Compute outlet temperature
+        T_0.append(temp(T_d, T_i, delta_z, X))
+
+    df.loc[delta_z] = T_0
+```
+
+## Numerical solution
+
+```python
+# Create output path 
+import os
+
+out_dir = os.environ.get('OGS_TESTRUNNER_OUT_DIR', '_out')
+if not os.path.exists(out_dir):
+    os.makedirs(out_dir)
+```
+
+```python
+from ogs6py.ogs import OGS
+
+# Create an OGS-object
+# Pass the project file and set an output file path
+model=OGS(INPUT_FILE="BHE_1P.prj", PROJECT_FILE=f"{out_dir}/BHE_1P.prj")
+
+# set end time to 30 days (in seconds)
+model.replace_text(30 * 24 * 60 * 60, xpath="./time_stepping/t_end/value")
+
+# Write to the output file 
+model.write_input()
+
+# Run OGS
+model.run_model(logfile=os.path.join(out_dir, "log.txt"), args=f"-o {out_dir} -m .")
+```
+
+## Results and discussion
+
+The outlet temperature change over time was compared against the analytical solution and presented in Figure 2.  
+After 30 days, the fluid temperature distribution in the wellbore is shown in Figure 3.  
+The maximum relative error between the numerical model and Ramey's analytical solution is less than 0.15%.
+
+In numerical model, the outlet temperature at beginning stage is affected by the initial temperature in the pipe inside the wellbore.  
+The initial fluid temperature set in the benchmark means there is water with 20°C filled in the wellbore already before injecting water into the wellbore.  
+But in the analytical solution, no initial temperature is set and the temperature keeps equilibrium state at every moment.  
+The impact of initial temperature condition in numerical model is decreasing with the increasing of the operational time as shown in Figure 2.
+
+```python
+import vtuIO
+import numpy as np
+import os
+
+out_dir = os.environ.get('OGS_TESTRUNNER_OUT_DIR', '_out')
+
+# Load the result
+pvdfile = vtuIO.PVDIO(f"{out_dir}/BHE_1P.pvd", dim=3)
+
+# Get point field names
+fields = pvdfile.get_point_field_names()
+print(fields)
+
+# Get all written timesteps
+time = pvdfile.timesteps
+
+# Define observation points
+observation_points = {'outlet':(0.0,10.0,-30.0)}
+
+# Read the time series of point field 'temperature_BHE1'
+temperature_BHE1 = pvdfile.read_time_series('temperature_BHE1', observation_points)
+```
+
+```python
+# Compute analytical solution again but at timesteps of numerical solution for error calculation
+for delta_z in Z:
+    T_0 = []
+    for delta_t in time:
+        
+        # Compute dimensionless time
+        t_D = dimless_t(lambda_re, delta_t, rho_re, c_p_re, r_b)
+
+        # Compute time function
+        f_t = time_function(t_D)
+
+        # Compute X
+        X = coefficient_x(q, rho_f, c_p_f, lambda_re, r_pi, U, f_t)
+
+        # Compute outlet temperature
+        T_0.append(temp(T_d, T_i, delta_z, X))
+```
+
+```python
+# Plot numerical vs. analytical result
+fig, ax1 = plt.subplots(figsize=(10, 8))
+
+x_pos = np.arange(0, 2592000+60*60*24*5, 60*60*24*5)
+x_ticks = np.arange(0,35,5)
+
+
+ax1.plot(time_range, df.iloc[-1, :], 'k.', markersize=10, markerfacecolor='none', label = 'Ramey\'s analytical solution')
+ax1.plot(time, temperature_BHE1['outlet'][:,0], 'r', label = 'OGS')
+ax1.set_ylabel('Outlet temperature (°C)', fontsize=20)
+ax1.set_xlabel('Time (d)', fontsize=20)
+ax1.set_xticks(x_pos, x_ticks)
+ax1.spines['left'].set_linewidth(2)
+ax1.spines['bottom'].set_linewidth(2)
+ax1.tick_params(axis='both', labelsize=16)
+
+color = 'blue'
+ax2 = ax1.twinx()
+ax2.plot(time ,temperature_BHE1['outlet'][:,0]-T_0, color=color, label = 'Absolute error')
+ax2.set_ylabel('Absolute error (°C)', color=color, fontsize = 20)
+ax2.tick_params(axis='y', labelcolor=color)
+ax2.spines['right'].set_color('blue')
+ax2.spines['right'].set_linewidth(2)
+ax2.tick_params(axis='y', labelsize=16)
+
+plt.figlegend(fontsize=18, loc = (0.4, 0.1), frameon=False)
+fig.tight_layout()
+
+plt.show()
+```
+
+Figure 2: Comparison with analytical solution results
+
+```python
+# Extract values along a line in the domain
+
+# Space axis for plotting
+length = np.linspace(0, -30, 101)
+
+# Draws a line through the domain for sampling results
+x_axis=[(0,10,i) for i in length]
+
+# At timestep i
+i = -1
+temp = pvdfile.read_set_data(time[i], "temperature_BHE1", pointsetarray=x_axis)
+
+# Plot
+fig, ax1 = plt.subplots(figsize=(10, 8))
+ax1.set_xlabel('Distance (m)', fontsize = 20)
+ax1.set_ylabel('Temperature (°C)', fontsize = 20)
+ax1.plot(np.linspace(0,30,101), temp[:,0], 'r', linewidth = 2, label = 'OGS')
+ax1.plot(Z,df.iloc[:, -1],'k.', markersize=10, markerfacecolor='none', label = 'Ramey\'s analytical solution')
+ax1.set_ylim(19,25)
+ax1.set_xlim(0,30)
+ax1.spines['left'].set_linewidth(2)
+ax1.spines['bottom'].set_linewidth(2)
+ax1.tick_params(axis='both', labelsize=16)
+
+color = 'blue'
+ax2 = ax1.twinx()
+ax2.set_ylabel('Absolute error (°C)', color=color, fontsize = 20)
+ax2.plot(np.linspace(0,30,101), temp[:,0]-df.iloc[:, -1], linewidth = 2, color=color, label = 'Absolute error')
+ax2.tick_params(axis='y', labelcolor=color)
+ax2.set_ylim(0,0.04)
+ax2.spines['right'].set_color('blue')
+ax2.spines['right'].set_linewidth(2)
+ax2.tick_params(axis='y', labelsize=16)
+
+fig.tight_layout()
+plt.figlegend(fontsize=18, loc = (0.4, 0.1), frameon=False)
+plt.show()
+
+```
+
+Figure 3: Distributed temperature of fluid and absolute error.
+
+## References
+
+<!-- vale off -->
+
+[1] Ramey Jr, H. J. (1962). Wellbore heat transmission. Journal of petroleum Technology, 14(04), 427-435.
+
+[2] Diersch, H-JG, et al. (2011). "Finite element modeling of borehole heat exchanger systems: Part 1. Fundamentals." Computers & Geosciences 37.8, 1122-1135.
+
+[3] Chaofan Chen, Wanlong Cai, Dmitri Naumov, Kun Tu, Hongwei Zhou, Yuping Zhang, Olaf Kolditz, Haibing Shao (2021). Numerical investigation on the capacity and efficiency of a deep enhanced U-tube borehole heat exchanger system for building heating. Renewable Energy, 169, 557-572.
+
+[4] Churchill, S. W. (1977). Comprehensive correlating equations for heat, mass and momentum transfer in fully developed flow in smooth tubes. Industrial & Engineering Chemistry Fundamentals, 16(1), 109-116.
+
+<!-- vale on -->
diff --git a/web/content/docs/benchmarks/heat-transport-bhe/pipe_flow_EBHE/Analytical_wellbore_heat_transport.zip b/web/content/docs/benchmarks/heat-transport-bhe/pipe_flow_EBHE/Analytical_wellbore_heat_transport.zip
deleted file mode 100644
index c0082f1983d..00000000000
Binary files a/web/content/docs/benchmarks/heat-transport-bhe/pipe_flow_EBHE/Analytical_wellbore_heat_transport.zip and /dev/null differ
diff --git a/web/content/docs/benchmarks/heat-transport-bhe/pipe_flow_EBHE/T_out_comparison.png b/web/content/docs/benchmarks/heat-transport-bhe/pipe_flow_EBHE/T_out_comparison.png
deleted file mode 100644
index a8f75f4fba8..00000000000
Binary files a/web/content/docs/benchmarks/heat-transport-bhe/pipe_flow_EBHE/T_out_comparison.png and /dev/null differ
diff --git a/web/content/docs/benchmarks/heat-transport-bhe/pipe_flow_EBHE/absolute_error_fluid_T_30d.png b/web/content/docs/benchmarks/heat-transport-bhe/pipe_flow_EBHE/absolute_error_fluid_T_30d.png
deleted file mode 100644
index b188b0ab1a0..00000000000
Binary files a/web/content/docs/benchmarks/heat-transport-bhe/pipe_flow_EBHE/absolute_error_fluid_T_30d.png and /dev/null differ