diff --git a/examples/script_design_capacity.py b/examples/script_design_capacity.py
index 4279abba18a266ce9693adba48cf371fb0f27628..d17120e454f2871001ee11c735859c6e05b4c8f8 100644
--- a/examples/script_design_capacity.py
+++ b/examples/script_design_capacity.py
@@ -4,26 +4,35 @@
 # TIP: this script requires the data found in the hyhetra-pipedata and hyhetra-
 # -fluiddata repositories
 
-pipedata_files = ['pipes/isoplus_single_disconti_s1.csv']
-        
+# pipedata_files = ['pipes/isoplus_single_disconti_s1.csv'] 
+pipedata_files = [
+    'pipes/enerpipe_caldopex_single.csv',
+    'pipes/enerpipe_caldopex_twin.csv',
+    'pipes/isoplus_single_disconti_s1.csv',
+    'pipes/isoplus_single_disconti_s2.csv',
+    'pipes/isoplus_single_disconti_s3.csv',
+    'pipes/isoplus_twin_disconti_s1.csv',
+    'pipes/isoplus_twin_disconti_s2.csv',
+    'pipes/isoplus_twin_disconti_s3.csv',
+    ]
 fluiddata_file = 'fluids/incropera2006_saturated_water.csv'
 
 # other files
-
 file_best2018 = 'papers/Best2018_design_flow_speeds.csv'
-
 file_sven2013 = 'papers/Sven2013_design_flow_speeds.csv'
-
 file_nussbaumer2016 = 'papers/Nussbaumer2016_design_flow_speeds.csv'
-
 file_euroheat2021 = 'papers/EuroHeat2021_design_flow_speeds.csv'
 
 # import libraries
+import math
 import numpy as np
 import matplotlib.pyplot as plt
 # topupheat
 import topupheat.pipes.fic as fic
+import topupheat.pipes.trenches as _tre
 from topupheat.pipes.single import StandardisedPipe, StandardisedPipeDatabase
+from topupheat.pipes.twin import StandardisedTwinPipe
+from topupheat.pipes.trenches import TwinPipeTrench
 import topupheat.pipes.utils as utils
 from topupheat.common import formulas as core
 
@@ -71,21 +80,21 @@ list_pipe_tuples = [(20, 1, 'isoplus', 'DRE-20-STD'),
                     (50, 1, 'isoplus', 'DRE-50-STD'),
                     (65, 1, 'isoplus', 'DRE-65-STD'),
                     (80, 1, 'isoplus', 'DRE-80-STD'),
-                    (100, 1, 'isoplus', 'DRE-100-ST'),
-                    (125, 1, 'isoplus', 'DRE-125-ST'),
-                    (150, 1, 'isoplus', 'DRE-150-ST'),
-                    (200, 1, 'isoplus', 'DRE-200-ST'),
-                    (250, 1, 'isoplus', 'DRE-250-ST'),
-                    (300, 1, 'isoplus', 'DRE-300-ST'),
-                    (350, 1, 'isoplus', 'DRE-350-ST'),
-                    (400, 1, 'isoplus', 'DRE-400-ST'),
-                    (450, 1, 'isoplus', 'DRE-450-ST'),
-                    (500, 1, 'isoplus', 'DRE-500-ST'),
-                    (600, 1, 'isoplus', 'DRE-600-ST'),
-                    (700, 1, 'isoplus', 'DRE-700-ST'),
-                    (800, 1, 'isoplus', 'DRE-800-ST'),
-                    (900, 1, 'isoplus', 'DRE-900-ST'),
-                    (1000, 1, 'isoplus', 'DRE-1000-S')]
+                    (100, 1, 'isoplus', 'DRE-100-STD'),
+                    (125, 1, 'isoplus', 'DRE-125-STD'),
+                    (150, 1, 'isoplus', 'DRE-150-STD'),
+                    (200, 1, 'isoplus', 'DRE-200-STD'),
+                    (250, 1, 'isoplus', 'DRE-250-STD'),
+                    (300, 1, 'isoplus', 'DRE-300-STD'),
+                    (350, 1, 'isoplus', 'DRE-350-STD'),
+                    (400, 1, 'isoplus', 'DRE-400-STD'),
+                    (450, 1, 'isoplus', 'DRE-450-STD'),
+                    (500, 1, 'isoplus', 'DRE-500-STD'),
+                    (600, 1, 'isoplus', 'DRE-600-STD'),
+                    (700, 1, 'isoplus', 'DRE-700-STD'),
+                    (800, 1, 'isoplus', 'DRE-800-STD'),
+                    (900, 1, 'isoplus', 'DRE-900-STD'),
+                    (1000, 1, 'isoplus', 'DRE-1000-STD')]
 
 #******************************************************************************
 #******************************************************************************
@@ -95,12 +104,10 @@ list_pipe_tuples = [(20, 1, 'isoplus', 'DRE-20-STD'),
 #******************************************************************************
 
 # pipe data
-
 pipedb = StandardisedPipeDatabase(source=pipedata_files)
 
 # water data
-
-waterdb = fic.FluidDatabase(fluid='water', phase='l', source=fluiddata_file)
+fluiddb = fic.FluidDatabase(fluid='water', phase='l', source=fluiddata_file)
 
 #******************************************************************************
 #******************************************************************************
@@ -111,19 +118,19 @@ dh_flow = fic.Fluid(
     phase='l',
     temperature=dh_flow_temperature,
     pressure=1e5,
-    db=waterdb)
+    db=fluiddb)
 
 dh_return = fic.Fluid(
     phase='l',
     temperature=dh_return_temperature,
     pressure=1e5,
-    db=waterdb)
+    db=fluiddb)
 
 fluid_bulk = fic.Fluid(
     phase='l',
     temperature=temperature_fluid_bulk,
     pressure=1e5,
-    db=waterdb)
+    db=fluiddb)
             
 #******************************************************************************
 #******************************************************************************
@@ -542,6 +549,7 @@ for i, max_specific_pressure_loss in enumerate(list_max_specific_pressure_loss):
     ax.plot(list_pipe_diameters,
             [spl for spl in dict_specific_pressure_loss[
                 (max_specific_pressure_loss,list_pipe_roughness[i])]],
+            list_linestyle_our[i],
             label='MSP='+str(max_specific_pressure_loss)+' Pa/m')
 
     # ax.plot(data_best2018['x']*1e3,
@@ -570,7 +578,201 @@ plt.show()
 
 #******************************************************************************
 #******************************************************************************
-#******************************************************************************
-#******************************************************************************
-#******************************************************************************
-#******************************************************************************
\ No newline at end of file
+
+# ISOPLUS validation
+
+# @ 30K
+
+isoplus_twin_mlfs = {
+    21.7: [0.5, 1.1], # 20
+    27.3: [0.5, 1.1], # 25
+    36.0: [0.6, 1.2], # 32
+    41.9: [0.6, 1.2], # 40
+    53.9: [0.7, 1.4], # 50
+    69.7: [0.7, 1.4], # 65
+    82.5: [0.8, 1.6], # 80
+    107.1: [0.8, 1.6], # 100
+    132.5: [1.0, 1.8], # 125
+    160.3: [1.2, 2.1], # 150
+    210.1: [1.4, 2.4], # 200
+    }
+
+# volumetric flow rate m3/hour
+isoplus_twin_vfr = {
+    21.7: [0.703, 1.547], # 20
+    27.3: [1.148, 2.526], # 25
+    36.0: [2.348, 4.695], # 32
+    41.9: [3.151, 6.303], # 40
+    53.9: [5.879, 11.757], # 50
+    69.7: [9.781, 19.563], # 65
+    82.5: [15.395, 30.791], # 80
+    107.1: [25.945, 51.891], # 100
+    132.5: [49.639, 89.350], # 125
+    160.3: [87.185, 152.573], # 150
+    210.1: [174.732, 299.541], # 200
+    }
+
+# capacity in kW
+isoplus_twin_capacity = {
+    21.7: [25, 54],         # 20
+    27.3: [40, 88],         # 25
+    36.0: [82, 164],        # 32
+    41.9: [110, 220],       # 40
+    53.9: [205, 410],       # 50
+    69.7: [341, 683],       # 65
+    82.5: [537, 1074],      # 80
+    107.1: [905, 1811],     # 100
+    132.5: [1732, 3118],    # 125
+    160.3: [3042, 5324],    # 150
+    210.1: [6097, 10451],   # 200
+    }
+        
+# *****************************************************************************
+# *****************************************************************************
+
+# specifications
+
+# district heating details
+# dh_flow_temperature = 273.15+115 # K  # leads to problems in friction factor calculation
+# dh_return_temperature = 273.15+85 # K # leads to problems in friction factor calculation
+dh_flow_temperature = 273.15+70 # K
+dh_return_temperature = 273.15+40 # K
+
+# temperature fluid bulk
+temperature_fluid_bulk = 0.5*(dh_return_temperature+dh_flow_temperature)
+
+# ground temperature
+ground_temperature = 273.15+10 # K
+
+# pipe absolute/effective roughness
+pipe_e_eff = 0.01/1000 # m
+
+# pipe length
+pipe_length = 1000
+
+# pipe depth
+pipe_depth = 0.8
+
+# soil thermal conductivity
+soil_k = 1 #0.5
+    
+# pipe distance
+pipe_distance = 0.1 # plus diameter
+ground_air_heat_transfer_coefficient = math.inf
+
+# single pipe method
+# single_pipe_trench_config = (_tre.SHT_METHOD_DIRECT, _tre.TRGTPT_KRISCHER1936)
+# single_pipe_trench_config = (_tre.SHT_METHOD_SYM_ASYM, _tre.TRGTPT_MULTIPOLE_ZERO_ORDER)
+# single_pipe_trench_config = (_tre.SHT_METHOD_SYM_ASYM, _tre.TRGTPT_MULTIPOLE_FIRST_ORDER)
+# single_pipe_trench_config = (_tre.SHT_METHOD_DIRECT, _tre.TRGTPT_MULTIPOLE_ZERO_ORDER)
+single_pipe_trench_config = (_tre.SHT_METHOD_DIRECT, _tre.TRGTPT_MULTIPOLE_FIRST_ORDER)
+
+# twin pipe method
+# twin_pipe_trench_config = (_tre.SHT_METHOD_SYM_ASYM, _tre.TRGTPT_MULTIPOLE_ZERO_ORDER)
+# twin_pipe_trench_config = (_tre.SHT_METHOD_SYM_ASYM, _tre.TRGTPT_MULTIPOLE_FIRST_ORDER)
+twin_pipe_trench_config = (_tre.SHT_METHOD_SYM_ASYM, _tre.TRGTPT_TWO_MODEL_APPROX)
+        
+# *****************************************************************************
+# *****************************************************************************
+
+pipe_tuples = [
+    (20, 1, 'isoplus', 'DRD-20-STD'), 
+    (25, 1, 'isoplus', 'DRD-25-STD'), 
+    (32, 1, 'isoplus', 'DRD-32-STD'),
+    (40, 1, 'isoplus', 'DRD-40-STD'),
+    (50, 1, 'isoplus', 'DRD-50-STD'), 
+    (65, 1, 'isoplus', 'DRD-65-STD'), 
+    (80, 1, 'isoplus', 'DRD-80-STD'),
+    (100, 1, 'isoplus', 'DRD-100-STD'),
+    (125, 1, 'isoplus', 'DRD-125-STD'),
+    (150, 1, 'isoplus', 'DRD-150-STD'), 
+    (200, 1, 'isoplus', 'DRD-200-STD'),
+    ]
+
+list_twin_pipes = [
+    StandardisedTwinPipe(
+        pipe_tuple=pipe_tuple,
+        e_eff=pipe_e_eff, # override pipe absolute roughness
+        length=pipe_length, # override the pipe length
+        db=pipedb) 
+    for pipe_tuple in pipe_tuples
+    ]
+            
+# *****************************************************************************
+# *****************************************************************************
+
+twin_trench = TwinPipeTrench(
+    pipe_center_depth=[pipe_depth+pipe.d_cas/2 for pipe in list_twin_pipes],   
+    fluid_db=fluiddb, 
+    phase=fluiddb.fluid_LIQUID, 
+    pressure=[1e5 for i in range(len(list_twin_pipes))], 
+    supply_temperature=[dh_flow_temperature for i in range(len(list_twin_pipes))], 
+    return_temperature=[dh_return_temperature for i in range(len(list_twin_pipes))], 
+    max_specific_pressure_loss=[100 for i in range(len(list_twin_pipes))], 
+    supply_pipe=list_twin_pipes
+    )
+
+twin_trench2 = TwinPipeTrench(
+    pipe_center_depth=[pipe_depth+pipe.d_cas/2 for pipe in list_twin_pipes],   
+    fluid_db=fluiddb, 
+    phase=fluiddb.fluid_LIQUID, 
+    pressure=[1e5 for i in range(len(list_twin_pipes))], 
+    supply_temperature=[dh_flow_temperature for i in range(len(list_twin_pipes))], 
+    return_temperature=[dh_return_temperature for i in range(len(list_twin_pipes))], 
+    max_specific_pressure_loss=[200 for i in range(len(list_twin_pipes))], 
+    supply_pipe=list_twin_pipes
+    )
+
+# pipe capacity vs flow rate
+fig, ax = plt.subplots()
+
+# datasheet
+level_label_datasheet = ['low speed', 'high speed']
+level_linestyle_datasheet = ['rd-','bd-']
+for level in range(2):
+    
+    ax.semilogy([d for d in isoplus_twin_capacity.keys()],
+            [_cap[level] for d, _cap in isoplus_twin_capacity.items()],
+            level_linestyle_datasheet[level],
+            label=level_label_datasheet[level])
+    
+# model
+level_label_model = ['low speed (model)', 'high speed (model)']
+level_linestyle_model = ['ro--','bo--']
+fluid_density = twin_trench.supply_fluid[0].mass_density #1000 # Kg/m3
+for level in range(2):
+    # convert vfr to m3/s
+    m = [_v[level]*fluid_density/3600 for k, _v in isoplus_twin_vfr.items()]
+    trans_cap = []
+    for i, _m in enumerate(m):
+        _cap = twin_trench.heat_transfer_rate(mass_flow_rate=_m)
+        trans_cap.append(_cap[i]/1e3)
+    
+    ax.semilogy([pipe.d_int*1e3 for pipe in twin_trench.supply_pipe],
+            trans_cap,
+            level_linestyle_model[level],
+            markersize=15,
+            markerfacecolor='none',
+            label=level_label_model[level])
+# rated capacity @ 100 Pa/m
+cap = [_cap/1e3 for _cap in twin_trench.rated_heat_capacity()]
+ax.semilogy([pipe.d_int*1e3 for pipe in twin_trench.supply_pipe],
+        cap,
+        'k--',
+        label=str(twin_trench.max_specific_pressure_loss[0])+' Pa/m (model)')
+# rated capacity @ 200 Pa/m
+cap2 = [_cap/1e3 for _cap in twin_trench2.rated_heat_capacity()]
+ax.semilogy([pipe.d_int*1e3 for pipe in twin_trench2.supply_pipe],
+        cap2,
+        'k-',
+        label=str(twin_trench2.max_specific_pressure_loss[0])+' Pa/m (model)')
+    
+ax.set(xlabel='Pipe internal diameter [mm]', 
+       ylabel='Transmittable Capacity [kW]')
+
+ax.legend()
+
+ax.grid()
+ax.set(xlim=(10, 220),ylim=(0, 12e3))
+# fig.savefig("test.png")
+plt.show()
diff --git a/examples/script_moody_diagram.py b/examples/script_moody_diagram.py
index 1b17462b33b4268986eca7d7a209f8330c224036..df177a97e9b4b2e48327b70087048eb58e8d9177 100644
--- a/examples/script_moody_diagram.py
+++ b/examples/script_moody_diagram.py
@@ -26,10 +26,10 @@ list_pipe_tuples = [(20, 1, 'isoplus', 'DRE-20-STD'),
                     (50, 1, 'isoplus', 'DRE-50-STD'),
                     (65, 1, 'isoplus', 'DRE-65-STD'),
                     (80, 1, 'isoplus', 'DRE-80-STD'),
-                    (100, 1, 'isoplus', 'DRE-100-ST'),
-                    (125, 1, 'isoplus', 'DRE-125-ST'),
-                    (150, 1, 'isoplus', 'DRE-150-ST'),
-                    (200, 1, 'isoplus', 'DRE-200-ST')]
+                    (100, 1, 'isoplus', 'DRE-100-STD'),
+                    (125, 1, 'isoplus', 'DRE-125-STD'),
+                    (150, 1, 'isoplus', 'DRE-150-STD'),
+                    (200, 1, 'isoplus', 'DRE-200-STD')]
 
 # relative pipe roughness
 
diff --git a/examples/script_nikuradse_experiments.py b/examples/script_nikuradse_experiments.py
index c77a9ecaaa3577b6d2e1831ae091657ccd9bd21c..b38d29061e06f6877018c44a095e8e3e16356b8e 100644
--- a/examples/script_nikuradse_experiments.py
+++ b/examples/script_nikuradse_experiments.py
@@ -38,10 +38,10 @@ list_pipe_tuples = [(20, 1, 'isoplus', 'DRE-20-STD'),
                     (50, 1, 'isoplus', 'DRE-50-STD'),
                     (65, 1, 'isoplus', 'DRE-65-STD'),
                     (80, 1, 'isoplus', 'DRE-80-STD'),
-                    (100, 1, 'isoplus', 'DRE-100-ST'),
-                    (125, 1, 'isoplus', 'DRE-125-ST'),
-                    (150, 1, 'isoplus', 'DRE-150-ST'),
-                    (200, 1, 'isoplus', 'DRE-200-ST')]
+                    (100, 1, 'isoplus', 'DRE-100-STD'),
+                    (125, 1, 'isoplus', 'DRE-125-STD'),
+                    (150, 1, 'isoplus', 'DRE-150-STD'),
+                    (200, 1, 'isoplus', 'DRE-200-STD')]
 
 # relative pipe roughness
 
diff --git a/examples/script_pipe_losses.py b/examples/script_pipe_losses.py
index ba7a385cd4066f17b9c26a559aa7748cdf24ca1e..1f2478bfa9448a305ca49aef1989da81a4aa96bb 100644
--- a/examples/script_pipe_losses.py
+++ b/examples/script_pipe_losses.py
@@ -10,6 +10,7 @@ from topupheat.pipes.single import StandardisedPipe, StandardisedPipeDatabase
 from topupheat.pipes.twin import StandardisedTwinPipe
 from topupheat.pipes.trenches import SupplyReturnPipeTrench, TwinPipeTrench
 import topupheat.pipes.trenches as _tre
+from topupheat.pipes.utils import minimum_assembling_distance_isoplus
 
 # *****************************************************************************
 # *****************************************************************************
@@ -25,65 +26,6 @@ pipedata_files = [
     'pipes/isoplus_twin_disconti_s2.csv',
     'pipes/isoplus_twin_disconti_s3.csv',
     ]
-
-def minimum_assembling_distance(pipe_outer_diameter: float):
-    
-    if pipe_outer_diameter <= 0.065:
-        return 0.100
-    elif pipe_outer_diameter <= 0.075:
-        return 0.100    
-    elif pipe_outer_diameter <= 0.090:
-        return 0.100    
-    elif pipe_outer_diameter <= 0.110:
-        return 0.150    
-    elif pipe_outer_diameter <= 0.125:
-        return 0.150    
-    elif pipe_outer_diameter <= 0.140:
-        return 0.150
-    elif pipe_outer_diameter <= 0.160:
-        return 0.200    
-    elif pipe_outer_diameter <= 0.180:
-        return 0.200    
-    elif pipe_outer_diameter <= 0.200:
-        return 0.200    
-    elif pipe_outer_diameter <= 0.225:
-        return 0.200
-    elif pipe_outer_diameter <= 0.250:
-        return 0.200    
-    elif pipe_outer_diameter <= 0.280:
-        return 0.300
-    elif pipe_outer_diameter <= 0.315:
-        return 0.300
-    elif pipe_outer_diameter <= 0.355:
-        return 0.300   
-    elif pipe_outer_diameter <= 0.400:
-        return 0.400    
-    elif pipe_outer_diameter <= 0.450:
-        return 0.400
-    elif pipe_outer_diameter <= 0.500:
-        return 0.400    
-    elif pipe_outer_diameter <= 0.560:
-        return 0.500    
-    elif pipe_outer_diameter <= 0.630:
-        return 0.500    
-    elif pipe_outer_diameter <= 0.670:
-        return 0.600
-    elif pipe_outer_diameter <= 0.710:
-        return 0.600
-    elif pipe_outer_diameter <= 0.800:
-        return 0.700
-    elif pipe_outer_diameter <= 0.900:
-        return 0.700
-    elif pipe_outer_diameter <= 1.000:
-        return 0.800
-    elif pipe_outer_diameter <= 1.100:
-        return 0.800
-    elif pipe_outer_diameter <= 1.200:
-        return 0.900
-    elif pipe_outer_diameter <= 1.300:
-        return 0.900
-    else:
-        raise NotImplementedError
         
 # *****************************************************************************
 # *****************************************************************************
@@ -94,8 +36,10 @@ def minimum_assembling_distance(pipe_outer_diameter: float):
 list_max_specific_pressure_loss = [80] # Pa/m
 
 # district heating details
-dh_flow_temperature = 273.15+115 # K
-dh_return_temperature = 273.15+85 # K
+# dh_flow_temperature = 273.15+115 # K  # leads to problems in friction factor calculation
+# dh_return_temperature = 273.15+85 # K # leads to problems in friction factor calculation
+dh_flow_temperature = 273.15+70 # K
+dh_return_temperature = 273.15+40 # K
 
 # temperature fluid bulk
 temperature_fluid_bulk = 0.5*(dh_return_temperature+dh_flow_temperature)
@@ -112,7 +56,7 @@ pipe_length = 1000
 # pipe depth
 pipe_depth = 0.8
 
-# pipe soil
+# soil thermal conductivity
 soil_k = 1 #0.5
     
 # pipe distance
@@ -399,8 +343,7 @@ u_pipes = {
     (150, 3, 'isoplus', 'DRD-150-TWIRE'): 0.2379
     }
 
-# pipe speeds
-
+# pipe mass flow in tones per hour
 m_pipes = { # 60-80 Pa
     21.7: [0.4, 0.5], # 20
     27.3: [0.8, 1.0],
@@ -632,7 +575,7 @@ for max_specific_pressure_loss in list_max_specific_pressure_loss:
         single_trench = SupplyReturnPipeTrench(
             pipe_center_depth=[pipe_depth+pipe.d_cas/2 for pipe in list_single_pipes],  
             # pipe_center_distance=[pipe_distance+pipe.d_cas for pipe in list_single_pipes], # minimum_assembling_distance
-            pipe_center_distance=[minimum_assembling_distance(pipe.d_cas)+pipe.d_cas for pipe in list_single_pipes],
+            pipe_center_distance=[minimum_assembling_distance_isoplus(pipe.d_cas)+pipe.d_cas for pipe in list_single_pipes],
             fluid_db=fluiddb, 
             phase=fluiddb.fluid_LIQUID, 
             pressure=[1e5 for i in range(len(list_single_pipes))], 
@@ -753,29 +696,25 @@ twin_pipe_tuples_per_tech = {
 for max_specific_pressure_loss in list_max_specific_pressure_loss:
     # single pipe
     for tech_index, _single_pipe_tuples in single_pipe_tuples_per_tech.items():
-        
-        _indices = [list_single_pipe_tuples.index(pipe_tuple) for pipe_tuple in _single_pipe_tuples]
-        
+        _indices = [list_single_pipe_tuples.index(pipe_tuple) for pipe_tuple in _single_pipe_tuples]        
         ax.semilogy(
             [1000*list_single_pipes[index].d_int for index in _indices], 
             [2*u_pipes[pipe_tuple]*(temperature_fluid_bulk-ground_temperature)
              for pipe_tuple in _single_pipe_tuples],
             style_category_single_official[tech_index],
             markersize=10,
-            label='isoplus, 2x single, tech '+str(model_string_ids[tech_index][1:])+' (datasheet)'
+            label='isoplus, 2x single, '+str(model_string_ids[tech_index][1:])+' series (datasheet)'
             )
     # twin pipe  
-    for tech_index, _twin_pipe_tuples in twin_pipe_tuples_per_tech.items():
-        
-        _indices = [list_twin_pipe_tuples.index(pipe_tuple) for pipe_tuple in _twin_pipe_tuples]
-        
+    for tech_index, _twin_pipe_tuples in twin_pipe_tuples_per_tech.items():        
+        _indices = [list_twin_pipe_tuples.index(pipe_tuple) for pipe_tuple in _twin_pipe_tuples]        
         ax.semilogy(
             [1000*list_twin_pipes[index].d_int for index in _indices], 
             [u_pipes[pipe_tuple]*(temperature_fluid_bulk-ground_temperature)
              for pipe_tuple in _twin_pipe_tuples],
             style_category_twin_official[tech_index],
             markersize=10,
-            label='isoplus, 1x twin, tech '+str(model_string_ids[tech_index][1:])+' (datasheet)'
+            label='isoplus, 1x twin, '+str(model_string_ids[tech_index][1:])+' series (datasheet)'
             )
         
 # model, single
@@ -797,26 +736,12 @@ for tech_index, _twin_pipe_tuples in twin_pipe_tuples_per_tech.items():
             markersize=10,
             label='isoplus, 1x twin (model)')
 
-# ax.plot([1000*pipe.d_int for pipe in list_single_pipes],
-#         [dict_single_heat_transfer_rate[max_specific_pressure_loss][i]
-#          for i, pipe_tuple in enumerate(list_single_pipe_tuples)],
-#         'gx',
-#         markersize=10,
-#         label='isoplus, 2x single (model)')
-
-# ax.plot([1000*pipe.d_int for pipe in list_twin_pipes],
-#         [dict_twin_heat_transfer_rate[max_specific_pressure_loss][i]
-#          for i, pipe_tuple in enumerate(list_twin_pipe_tuples)],
-#         'k^',
-#         markersize=10,
-#         label='isoplus, 1x twin (model)')
-
 ax.set(xlabel='Nominal diameter (DN)', 
        ylabel='Specific heat losses [W/m]',
        title='Specific heat losses as a function of pipe diameter')
 
-ax.grid()
-
+ax.set(xlim=(0, 1000),ylim=(5, 100))
+ax.grid(which='both')
 ax.legend()
 
 # fig.savefig("test.png")
@@ -830,7 +755,7 @@ plt.show()
 fig, ax = plt.subplots()
 fig.set_size_inches(15,10)
 
-cpdt = 4187*30*(3600/1000)/1000 # m in tons per hour; cpdt has to be in J/Kg
+cpdt = 4187*(dh_flow_temperature-dh_return_temperature)*(1000/3600)/1000 # m in tons per hour; cpdt has to be in J/Kg
         
 max_specific_pressure_loss = 80
         
@@ -870,6 +795,10 @@ for max_specific_pressure_loss in list_max_specific_pressure_loss:
             markersize=10,
             label='isoplus, 2x single (model)')
         break
+    
+    # [dict_single_trench[max_specific_pressure_loss].nmls_supply
+    #  for pipe_tuple in _single_pipe_tuples]
+    
 
 # ax.semilogy([1000*pipe.d_int for pipe in list_twin_pipes],
 #         [dict_twin_capacity[max_specific_pressure_loss][i]/1e3
diff --git a/src/topupheat/pipes/system.py b/src/topupheat/pipes/system.py
index c446094533058bf95dc8f4d1cc4761880e9d72fb..04032d1562c57b67e3f1d83d76d5a3eb9e3ec30f 100644
--- a/src/topupheat/pipes/system.py
+++ b/src/topupheat/pipes/system.py
@@ -569,7 +569,6 @@ class SupplyReturnPipeSystem:
             unit_conversion_factor: float = 1.0
             ):
         
-        # length >> options
         # specific_heat_transfer_coefficient >> options
         # temperature_surroundings >> temporal
         # determine if the inputs are for temporal vector/normal mode
@@ -585,7 +584,8 @@ class SupplyReturnPipeSystem:
         if ((not self.vector_mode and 
              not isinstance(specific_heat_transfer_coefficient, Real)) or
             (self.vector_mode and 
-             len(self.supply_temperature) < 1
+             len(self.supply_temperature) < 1 or
+             len(specific_heat_transfer_coefficient) != self.number_options()
              )
             ):
             raise TypeError('Inconsistent inputs.')
@@ -646,7 +646,7 @@ class SupplyReturnPipeSystem:
         return out
     
     # *************************************************************************
-        
+    # TODO: reconsider the implementation in vector mode
     def simplified_heat_transfer_surroundings(
             self,
             length: float or list,
@@ -655,6 +655,40 @@ class SupplyReturnPipeSystem:
             time_interval_duration: float or list,
             unit_conversion_factor: float = 1.0
             ):
+        """
+        Calculates the heat transfer with the surroundings using a constant
+        specific heat transfer coefficient.
+
+        Parameters
+        ----------
+        length : float or list
+            The length of each segment.
+        specific_heat_transfer_coefficient : float or list
+            The heat transfer coefficient per unit length.
+        temperature_surroundings : float or list
+            The temperature of the surroundings.
+        time_interval_duration : float or list
+            The duration of each time step under consideration.
+        unit_conversion_factor : float, optional
+            A factor to adjust the final units. The default is 1.0.
+
+        Raises
+        ------
+        TypeError
+            This error is raised if the combination of inputs is incompatible.
+
+        Returns
+        -------
+        out : float, list or list of lists
+            The heat transfer rate based on the inputs. If single inputs are
+            provided and a single option is under consideration, then a float
+            is returned. If there are multiple options and single inputs are
+            provided, then a list with one heat transfer rate per option is
+            returned. If there are multiple options and multiple inputs, then
+            a list of lists is returned: one list with one list per option, 
+            each of which with one heat transfer rate per input.
+
+        """
         
         # length >> options
         # specific_heat_transfer_coefficient >> options
@@ -669,7 +703,7 @@ class SupplyReturnPipeSystem:
                type(time_interval_duration) == tuple)):
             temporal_vector_mode = True
         else:
-            raise ValueError('Incompatible inputs.')
+            raise TypeError('Incompatible inputs.')
         
         # confirm inputs are consistent with vector/normal mode
         if ((not self.vector_mode and 
@@ -678,7 +712,8 @@ class SupplyReturnPipeSystem:
             (self.vector_mode and 
              (type(length) != list and 
               type(length) != tuple or
-              len(length) != len(self.supply_temperature))
+              len(length) != self.number_options() or
+              len(specific_heat_transfer_coefficient) != len(length))
              )
             ):
             raise TypeError('Inconsistent inputs.')
diff --git a/tests/test_trenches.py b/tests/test_trenches.py
index d47da40f18c334219a6a9688cc2dcd97aa479d57..930873848d08ec86e070c5b534d428015033ad84 100644
--- a/tests/test_trenches.py
+++ b/tests/test_trenches.py
@@ -1,4 +1,4 @@
-
+import math
 import numpy as np
 # from numbers import Real
 from src.topupheat.common.fluids import FluidDatabase, Fluid
@@ -9,6 +9,7 @@ from src.topupheat.pipes.twin import SymmetricalTwinPipe, StandardisedTwinPipe
 from math import isclose, inf, ceil
 from src.topupheat.common.factors import shape_factor_buried_horizontal_isothermal_cylinder
 from src.topupheat.common.factors import thermal_resistance_from_shape_factor
+from topupheat.pipes.utils import minimum_assembling_distance_isoplus
 
 # *****************************************************************************
 # *****************************************************************************
@@ -2704,9 +2705,170 @@ class TestPipeTrench:
             #     )
             
             # assert type(trench.printable_description()) == str
+                    
+    # *************************************************************************
+    # *************************************************************************
+    
+    def test_simplified_heat_transfer_twin_pipes(self):
             
-    # TODO: test the simplified specific heat transfer method using trench objects
+        # maximum specific pressure drop
+        max_specific_pressure_loss = 80 # Pa/m
+
+        # district heating details
+        dh_flow_temperature = 273.15+70 # K
+        dh_return_temperature = 273.15+40 # K
+        
+        # ground temperature
+        ground_temperature = 273.15+10 # K
+
+        # pipe absolute/effective roughness
+        pipe_e_eff = 0.01/1000 # m
+
+        # pipe length
+        pipe_length = 1000
+
+        # pipe depth
+        pipe_depth = 0.8
+        
+        # pipe data
+        pipedata_files = [
+            'tests/data/isoplus_twin_disconti_s1.csv',
+            # 'tests/data/isoplus_twin_disconti_s2.csv',
+            # 'tests/data/isoplus_twin_disconti_s3.csv',
+            ]
+        
+        # pipe data
+        pipedb = StandardisedPipeDatabase(source=pipedata_files)
+
+        # fluid data   
+        fluiddata_file = 'tests/data/incropera2006_saturated_water.csv'
+        fluiddb = fic.FluidDatabase(fluid='water', phase='l', source=fluiddata_file)
+            
+        # *****************************************************************************
+        # *****************************************************************************
+
+        # pipe tuples
+        list_pipe_tuples = [
+            pipe_tuple
+            for pipe_tuple in pipedb.pipe_tuples
+            ]
+        # pipe objects
+        list_pipes = [
+            StandardisedTwinPipe(
+                pipe_tuple=pipe_tuple,
+                e_eff=pipe_e_eff, # override pipe absolute roughness
+                length=pipe_length, # override the pipe length
+                db=pipedb)
+            for pipe_tuple in list_pipe_tuples
+            ]
+        
+        # single pipe trench    
+        single_trench = trenches.TwinPipeTrench(
+            pipe_center_depth=[pipe_depth+pipe.d_cas/2 for pipe in list_pipes],  
+            fluid_db=fluiddb, 
+            phase=fluiddb.fluid_LIQUID, 
+            pressure=[1e5 for i in range(len(list_pipes))], 
+            supply_temperature=[dh_flow_temperature for i in range(len(list_pipes))], 
+            return_temperature=[dh_return_temperature for i in range(len(list_pipes))], 
+            max_specific_pressure_loss=[max_specific_pressure_loss for i in range(len(list_pipes))], 
+            supply_pipe=list_pipes
+            )
+        assert single_trench.vector_mode
+        assert single_trench.number_options() == len(list_pipes)
+        # one u per option
+        u = [i+1 for i in range(len(list_pipe_tuples))]
+        q = single_trench.simplified_specific_heat_transfer_surroundings(
+            specific_heat_transfer_coefficient=u, 
+            temperature_surroundings=ground_temperature, 
+            )
+        for i, q_i in enumerate(q):
+            assert math.isclose(
+                q_i, 
+                u[i]*(dh_flow_temperature/2+dh_return_temperature/2-ground_temperature),
+                abs_tol=1e-3
+                )
+        # TODO: multiple time intervals
+            
+    # *************************************************************************
+    # *************************************************************************
+     
+    def test_simplified_heat_transfer_two_single_pipes(self):
+        
+        # maximum specific pressure drop
+        max_specific_pressure_loss = 80 # Pa/m
+
+        # district heating details
+        dh_flow_temperature = 273.15+70 # K
+        dh_return_temperature = 273.15+40 # K
+        # ground temperature
+        ground_temperature = 273.15+10 # K
+
+        # pipe absolute/effective roughness
+        pipe_e_eff = 0.01/1000 # m
+
+        # pipe length
+        pipe_length = 1000
+
+        # pipe depth
+        pipe_depth = 0.8
+        
+        # pipe data
+        pipedata_files = [
+            'tests/data/isoplus_single_disconti_s1.csv',
+            # 'tests/data/isoplus_single_disconti_s2.csv',
+            # 'tests/data/isoplus_single_disconti_s3.csv',
+            ]
+        
+        # pipe data
+        pipedb = StandardisedPipeDatabase(source=pipedata_files)
+
+        # fluid data   
+        fluiddata_file = 'tests/data/incropera2006_saturated_water.csv'
+        fluiddb = fic.FluidDatabase(fluid='water', phase='l', source=fluiddata_file)
     
+        # pipe tuples
+        list_pipe_tuples = [
+            pipe_tuple
+            for pipe_tuple in pipedb.pipe_tuples
+            ]
+        
+        # pipe objects
+        list_pipes = [
+            StandardisedPipe(
+                pipe_tuple=pipe_tuple,
+                e_eff=pipe_e_eff, # override pipe absolute roughness
+                length=pipe_length, # override the pipe length
+                db=pipedb)
+            for pipe_tuple in list_pipe_tuples
+            ]
+        
+        # single pipe trench    
+        single_trench = trenches.SupplyReturnPipeTrench(
+            pipe_center_depth=[pipe_depth+pipe.d_cas/2 for pipe in list_pipes],  
+            pipe_center_distance=[minimum_assembling_distance_isoplus(pipe.d_cas)+pipe.d_cas for pipe in list_pipes],
+            fluid_db=fluiddb, 
+            phase=fluiddb.fluid_LIQUID, 
+            pressure=[1e5 for i in range(len(list_pipes))], 
+            supply_temperature=[dh_flow_temperature for i in range(len(list_pipes))], 
+            return_temperature=[dh_return_temperature for i in range(len(list_pipes))], 
+            max_specific_pressure_loss=[max_specific_pressure_loss for i in range(len(list_pipes))], 
+            supply_pipe=list_pipes
+            )
+        assert single_trench.vector_mode
+        # one u per option
+        u = [i+1 for i in range(len(list_pipe_tuples))]
+        q = single_trench.simplified_specific_heat_transfer_surroundings(
+            specific_heat_transfer_coefficient=u, 
+            temperature_surroundings=ground_temperature, 
+            )
+        for i, q_i in enumerate(q):
+            assert math.isclose(
+                q_i, 
+                u[i]*(dh_flow_temperature/2+dh_return_temperature/2-ground_temperature),
+                abs_tol=1e-3
+                )
+    
+    # *************************************************************************
     # *************************************************************************
         
     def test_different_single_pipes(self):
@@ -2852,98 +3014,6 @@ class TestPipeTrench:
     # *************************************************************************
     # *************************************************************************
     
-    # TODO: try to do this
-    
-    # def test_u_coefficient_formula_bohm2005_twin(self):
-        
-    #     # *********************************************************************
-        
-    #     # trench
-        
-    #     pipe_distance = 0.0251+0.0889 # m
-        
-    #     pipe_depth = 1.28 # m
-        
-    #     h_gs = 14.6 # W/m2K
-        
-    #     ground_temperature = 8+273.15 # K
-        
-    #     soil_k = 1.5 # W/mK
-        
-    #     supply_temperature = 80+273.15
-        
-    #     return_temperature = 40+273.15
-        
-    #     #**************************************************************************
-        
-    #     # pipes
-        
-    #     twin_pipe = SymmetricalTwinPipe(
-    #         length=1000,
-    #         k=60,
-    #         e_eff=1e-5,
-    #         d_int=0.0889-2*0.0032,
-    #         d_ext=0.0889,
-    #         d_ins=0.2483,
-    #         d_cas=0.2483+2*0.0036,
-    #         k_ins=0.0265,
-    #         k_cas=0.43,
-    #         pipe_center_distance=pipe_distance)
-  
-    #     #**************************************************************************
-        
-    #     u11_true = 0.2517 # W/mK
-        
-    #     u22_true = 0.2534 # W/mK
-        
-    #     u12_true = 0.0784 # W/mK
-        
-    #     q_total_true = 18.08 # W/m
-    
-    #     list_abs_tol_u11 = [1e-9, 1e-9, 1e-9]
-    #     list_abs_tol_u12 = [1e-9, 1e-9, 1e-9]
-        
-    #     # calculations
-        
-    #     list_ground_resistance_models = [
-    #         trenches.TRGTPT_TWO_MODEL_APPROX,
-    #         trenches.TRGTPT_MULTIPOLE_ZERO_ORDER,
-    #         trenches.TRGTPT_MULTIPOLE_FIRST_ORDER
-    #         ]
-        
-    #     for i, model in enumerate(list_ground_resistance_models):
-        
-    #         (rsym,
-    #          rasym) = trenches.ThermalResistanceBuriedTwinPipesWallenten(
-    #             twin_pipe,
-    #             soil_k,
-    #             pipe_depth,
-    #             pipe_distance,
-    #             h_gs,
-    #             model=model)
-            
-    #         (u1122,
-    #          u12) = trenches.HeatTransferCoefficientsFromSymmetricalAsymmetricalResistances(
-    #               rsym, rasym)
-    #         print('hey')
-    #         print(u1122)
-    #         print(u12)
-    #         assert isclose( 
-    #             u1122, 
-    #             u11_true,
-    #             abs_tol=list_abs_tol_u11[i]
-    #             )
-                                                                                          
-    #         assert isclose(
-    #             u12, 
-    #             u12_true,
-    #             abs_tol=list_abs_tol_u12[i]
-    #             )
-        
-    #     # *********************************************************************
-        
-    # *************************************************************************
-    
     def test_district_cooling(self):
         
         # load pipe data