diff --git a/tests/test_system.py b/tests/test_system.py
index 138ef4c17cd4f28f7b03086abceb38fee8c2465c..9d5c5ef59f816cdfc83689b06d5ed78936cd2389 100644
--- a/tests/test_system.py
+++ b/tests/test_system.py
@@ -1,11 +1,11 @@
-
+import math
 from numbers import Real
 from src.topupheat.common.fluids import FluidDatabase
 from src.topupheat.pipes.system import SupplyReturnPipeSystem
 from src.topupheat.pipes.single import StandardisedPipeDatabase
 from src.topupheat.pipes.single import Pipe, InsulatedPipe,StandardisedPipe
 from src.topupheat.pipes.twin import SymmetricalTwinPipe, StandardisedTwinPipe
-from math import isclose, inf, ceil
+from math import isclose, ceil
 
 # *****************************************************************************
 # *****************************************************************************
@@ -15,517 +15,518 @@ class TestPipeSystem:
     # TODO: these cases
     
     def test_simplified_heat_transfer_surroundings_nonvector_single_step(self):
-        
-        pass
-        
-    def test_simplified_heat_transfer_surroundings_nonvector_multi_step(self):
-        
-        pass
             
-    def test_simplified_heat_transfer_surroundings_vector_single_step(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
         
-        pass
+        # 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 data
+        pipedata_files = [
+            'tests/data/isoplus_twin_disconti_s1.csv',
+            # 'tests/data/isoplus_twin_disconti_s2.csv',
+            # 'tests/data/isoplus_twin_disconti_s3.csv',
+            ]
         
-    def test_simplified_heat_transfer_surroundings_vector_multi_step(self):
+        # pipe data
+        pipedb = StandardisedPipeDatabase(source=pipedata_files)
+
+        # fluid data   
+        fluiddata_file = 'tests/data/incropera2006_saturated_water.csv'
+        fluiddb = 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
+            ]
         
-        pass
-    
-    
-    def test_simplified_heat_transfer_surroundings(self):
+        # two pipe system object
+        two_pipe_system = SupplyReturnPipeSystem(
+            fluid_db=fluiddb, 
+            phase=fluiddb.fluid_LIQUID, 
+            pressure=1e5, 
+            supply_temperature=dh_flow_temperature, 
+            return_temperature=dh_return_temperature, 
+            max_specific_pressure_loss=max_specific_pressure_loss, 
+            supply_pipe=list_pipes[0]
+            )
+        assert not two_pipe_system.vector_mode
+        assert two_pipe_system.number_options() == 1
+        
+        # non-vector mode, one u value per option, and one temperature
+        u = 1
+        # q should be a number
+        q = two_pipe_system.simplified_specific_heat_transfer_surroundings(
+            specific_heat_transfer_coefficient=u, 
+            temperature_surroundings=ground_temperature, 
+            )
+        assert isinstance(q, Real)
+        assert 1 == two_pipe_system.number_options()
+        assert math.isclose(
+            q, 
+            u*(dh_flow_temperature/2+dh_return_temperature/2-ground_temperature),
+            abs_tol=1e-3
+            )
+        # test unit conversion
+        u_factor = 2
+        q = two_pipe_system.simplified_specific_heat_transfer_surroundings(
+            specific_heat_transfer_coefficient=u, 
+            temperature_surroundings=ground_temperature, 
+            unit_conversion_factor=u_factor
+            )
+        assert isinstance(q, Real)
+        assert 1 == two_pipe_system.number_options()
+        assert math.isclose(
+            q, 
+            u*(dh_flow_temperature/2+dh_return_temperature/2-ground_temperature)*u_factor,
+            abs_tol=1e-3
+            )
         
-        # load pipe data
-            
-        singlepipedata_files = ['tests/data/enerpipe_caldopex_single.csv']
-        singlepipedb = StandardisedPipeDatabase(source=singlepipedata_files)
+        # with length and time interval duration
+        length = 3
+        time_interval_duration = 10
+        q = two_pipe_system.simplified_heat_transfer_surroundings(
+            specific_heat_transfer_coefficient=u, 
+            temperature_surroundings=ground_temperature, 
+            time_interval_duration=time_interval_duration,
+            length=length
+            )
+        assert isinstance(q, Real)
+        assert 1 == two_pipe_system.number_options()
+        assert math.isclose(
+            q, 
+            u*(dh_flow_temperature/2+dh_return_temperature/2-ground_temperature)*length*time_interval_duration,
+            abs_tol=1e-3
+            )
+        # use unit conversion
+        u_factor = 2
+        q = two_pipe_system.simplified_heat_transfer_surroundings(
+            specific_heat_transfer_coefficient=u, 
+            temperature_surroundings=ground_temperature, 
+            time_interval_duration=time_interval_duration,
+            length=length,
+            unit_conversion_factor=u_factor
+            )
+        assert isinstance(q, Real)
+        assert 1 == two_pipe_system.number_options()
+        assert math.isclose(
+            q, 
+            u*(dh_flow_temperature/2+dh_return_temperature/2-ground_temperature)*length*time_interval_duration*u_factor,
+            abs_tol=1e-3
+            )
         
-        # water 
+    # *************************************************************************
+    # *************************************************************************
         
-        waterdata_file = 'tests/data/incropera2006_saturated_water.csv'
-        phase = FluidDatabase.fluid_LIQUID
-        fluid_db = FluidDatabase(
-            fluid='fluid',
-            phase=phase,
-            source=waterdata_file
-            )
+    def test_simplified_heat_transfer_surroundings_nonvector_multi_step(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
         
-        supply_temperature = 80+273.15
-        return_temperature = 50+273.15
-        pressure = 1e5
-        max_specific_pressure_loss = 100 # Pa/m
-        temperature_surroundings = 10+273.15
-        
-        pipe_tuples = [
-            (20,2),     # 1
-            (25,2),     # 2
-            (32,2),     # 3
-            (40,2),     # 4 
-            (50,2),     # 5 
-            (65,2),     # 6 twin pipe too! the DN is all that matters here
-            (80,2),     # 7
-            (100,2),    # 8
-            (125,2),    # 9
-            (150,1)     # 11
-            ]
-    
-        pipes = {
-            pipe_tuple[0:2]: StandardisedPipe(
-                pipe_tuple=pipe_tuple,
-                db=singlepipedb)
-            for pipe_tuple in singlepipedb.pipe_tuples
-            }
+        # ground temperature
+        ground_temperature = 273.15+10 # K
 
-        u_coefficients = [
-            0.1392/2,   # 1
-            0.1571/2,   # 2
-            0.1741/2,   # 3
-            0.1662/2,   # 4
-            0.2075/2,   # 5
-            0.2802/2,   # 6
-            0.1747,     # 7
-            0.1992,     # 8
-            0.2771,     # 9 
-            0.3028      # 11
-            ]
-                
-        # specific heat losses extracted graphically from the paper
-        specific_heat_losses = [
-            7.653,      # 1
-            8.6724,     # 2
-            9.5565,     # 3
-            9.1742,     # 4
-            11.3841,    # 5
-            15.4729,    # 6
-            19.2494,    # 7
-            21.9193,    # 8
-            30.4776,    # 9
-            33.3056     # 11
+        # pipe absolute/effective roughness
+        pipe_e_eff = 0.01/1000 # m
+
+        # pipe length
+        pipe_length = 1000
+        
+        # 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',
             ]
         
-        specific_heat_losses_tol = [
-            3.1e-3, # 1 : 0.0030000000000001137
-            0.0320,  # 2 : 0.03190000000000026
-            0.0191,  # 3 : 0.019000000000000128
-            0.0333,  # 4 : 0.033200000000000784
-            0.0284,  # 5 : 0.028399999999999537
-            0.0619, # 6 : 0.06189999999999962
-            0.0325, # 7 : 0.03240000000000265
-            0.0074, # 8 : 0.00730000000000075
-            0.0035, # 9 : 0.003400000000002734
-            0.0025, # 11: 0.002400000000001512
+        # pipe data
+        pipedb = StandardisedPipeDatabase(source=pipedata_files)
+
+        # fluid data   
+        fluiddata_file = 'tests/data/incropera2006_saturated_water.csv'
+        fluiddb = 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
             ]
         
-        # pipe
-        for pipe_tuple, sht_true, sht_tol, u_coef in zip(
-                pipe_tuples, 
-                specific_heat_losses, 
-                specific_heat_losses_tol,
-                u_coefficients
-                ):
-            
-            srps = SupplyReturnPipeSystem(
-                fluid_db=fluid_db,
-                phase=phase, 
-                pressure=pressure, 
-                supply_temperature=supply_temperature, 
-                return_temperature=return_temperature, 
-                supply_pipe=pipes[pipe_tuple],
-                max_specific_pressure_loss=max_specific_pressure_loss # Pa/m
+        # two pipe system object
+        two_pipe_system = SupplyReturnPipeSystem(
+            fluid_db=fluiddb, 
+            phase=fluiddb.fluid_LIQUID, 
+            pressure=1e5, 
+            supply_temperature=dh_flow_temperature, 
+            return_temperature=dh_return_temperature, 
+            max_specific_pressure_loss=max_specific_pressure_loss, 
+            supply_pipe=list_pipes[0]
+            )
+        assert not two_pipe_system.vector_mode
+        assert two_pipe_system.number_options() == 1
+        
+        # non-vector mode, one u value per time interval (and option), multiple temperatures
+        number_time_intervals = 3
+        u = [1+k for k in range(number_time_intervals)]
+        t_surroundings = [ground_temperature+k for k in range(number_time_intervals)]
+        # q should be a list of lists
+        q = two_pipe_system.simplified_specific_heat_transfer_surroundings(
+            specific_heat_transfer_coefficient=u, 
+            temperature_surroundings=t_surroundings, 
+            )
+        assert type(q) == list or type(q) == tuple
+        assert len(q) == number_time_intervals
+        for u_i, q_i, t_i in zip(u, q, t_surroundings):
+            assert isinstance(q_i, Real)
+            assert math.isclose(
+                q_i, 
+                u_i*(dh_flow_temperature/2+dh_return_temperature/2-t_i),
+                abs_tol=1e-3
                 )
-            
-            sht = srps.simplified_heat_transfer_surroundings(
-                length=1, 
-                specific_heat_transfer_coefficient=u_coef, 
-                temperature_surroundings=temperature_surroundings, 
-                time_interval_duration=1
+        # use unit conversion
+        u_factor = 2
+        q = two_pipe_system.simplified_specific_heat_transfer_surroundings(
+            specific_heat_transfer_coefficient=u, 
+            temperature_surroundings=t_surroundings, 
+            unit_conversion_factor=u_factor
+            )
+        assert type(q) == list or type(q) == tuple
+        assert len(q) == number_time_intervals
+        for u_i, q_i, t_i in zip(u, q, t_surroundings):
+            assert isinstance(q_i, Real)
+            assert math.isclose(
+                q_i, 
+                u_i*(dh_flow_temperature/2+dh_return_temperature/2-t_i)*u_factor,
+                abs_tol=1e-3
                 )
+    
+    # *************************************************************************
+    # *************************************************************************
             
-            assert isclose(sht, sht_true, abs_tol=sht_tol)
+    def test_simplified_heat_transfer_surroundings_vector_single_step(self):
         
-            sht = srps.simplified_specific_heat_transfer_surroundings(
-                specific_heat_transfer_coefficient=u_coef, 
-                temperature_surroundings=temperature_surroundings, 
-                )
-            
-            assert isclose(sht, sht_true, abs_tol=sht_tol)
-            
-            number_steps = 2
-        
-            sht = srps.simplified_heat_transfer_surroundings(
-                length=1, 
-                specific_heat_transfer_coefficient=u_coef, 
-                temperature_surroundings=[
-                    temperature_surroundings for i in range(number_steps)], 
-                time_interval_duration=[
-                    1 for i in range(number_steps)], 
-                )
-            assert type(sht) == list
-            assert len(sht) == number_steps
-            for _sht in sht:
-                assert isclose(_sht, sht_true, abs_tol=sht_tol)
-                
-            sht = srps.simplified_specific_heat_transfer_surroundings(
-                specific_heat_transfer_coefficient=u_coef, 
-                temperature_surroundings=[
-                    temperature_surroundings for i in range(number_steps)], 
-                )
-            assert type(sht) == list
-            assert len(sht) == number_steps
-            for _sht in sht:
-                assert isclose(_sht, sht_true, abs_tol=sht_tol)
-                
-            # try using incorrect option-based inputs
-                
-            error_raised = False
-            try:
-                srps.simplified_heat_transfer_surroundings(
-                    length=[1], 
-                    specific_heat_transfer_coefficient=[u_coef], 
-                    temperature_surroundings=temperature_surroundings, 
-                    time_interval_duration=1
-                    )
-            except TypeError:
-                error_raised = True
-            assert error_raised
-                            
-            error_raised = False
-            try:
-                srps.simplified_specific_heat_transfer_surroundings(
-                    specific_heat_transfer_coefficient=[u_coef], 
-                    temperature_surroundings=temperature_surroundings
-                    )
-            except TypeError:
-                error_raised = True
-            assert error_raised
-                
-            # *****************************************************************
-            
-            # try using other units
-            ucf = 2
-            sht = srps.simplified_heat_transfer_surroundings(
-                length=1, 
-                specific_heat_transfer_coefficient=u_coef, 
-                temperature_surroundings=temperature_surroundings, 
-                time_interval_duration=1,
-                unit_conversion_factor=ucf
+        # 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 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 = 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
+            ]
+        number_options = len(list_pipes)
+        # two pipe system
+        two_pipe_system = SupplyReturnPipeSystem(
+            fluid_db=fluiddb, 
+            phase=fluiddb.fluid_LIQUID, 
+            pressure=[1e5 for i in range(number_options)],
+            supply_temperature=[dh_flow_temperature for i in range(number_options)], 
+            return_temperature=[dh_return_temperature for i in range(number_options)],  
+            max_specific_pressure_loss=[max_specific_pressure_loss for i in range(number_options)], 
+            supply_pipe=list_pipes
+            )
+        assert two_pipe_system.vector_mode
+        assert two_pipe_system.number_options() == len(list_pipes)
+        
+        # vector mode, one u value per option, and one temperature
+        u = [i+1 for i in range(len(list_pipe_tuples))]
+        # q should be a list of numbers
+        q = two_pipe_system.simplified_specific_heat_transfer_surroundings(
+            specific_heat_transfer_coefficient=u, 
+            temperature_surroundings=ground_temperature, 
+            )
+        assert type(q) == list or type(q) == tuple
+        assert len(q) == two_pipe_system.number_options()
+        for u_i, q_i in zip(u, q):
+            assert math.isclose(
+                q_i, 
+                u_i*(dh_flow_temperature/2+dh_return_temperature/2-ground_temperature),
+                abs_tol=1e-3
                 )
-            
-            assert isclose(sht, sht_true*ucf, abs_tol=sht_tol*ucf)
-            
-            ucf = 2
-            sht = srps.simplified_specific_heat_transfer_surroundings(
-                specific_heat_transfer_coefficient=u_coef, 
-                temperature_surroundings=temperature_surroundings, 
-                unit_conversion_factor=ucf
+        # use unit conversion
+        u_factor = 2        
+        q = two_pipe_system.simplified_specific_heat_transfer_surroundings(
+            specific_heat_transfer_coefficient=u, 
+            temperature_surroundings=ground_temperature, 
+            unit_conversion_factor=u_factor
+            )
+        assert type(q) == list or type(q) == tuple
+        assert len(q) == two_pipe_system.number_options()
+        for u_i, q_i in zip(u, q):
+            assert math.isclose(
+                q_i, 
+                u_i*(dh_flow_temperature/2+dh_return_temperature/2-ground_temperature)*u_factor,
+                abs_tol=1e-3
                 )
-            
-            assert isclose(sht, sht_true*ucf, abs_tol=sht_tol*ucf)
-            
-            number_steps = 2
-        
-            sht = srps.simplified_heat_transfer_surroundings(
-                length=1, 
-                specific_heat_transfer_coefficient=u_coef, 
-                temperature_surroundings=[
-                    temperature_surroundings for i in range(number_steps)], 
-                time_interval_duration=[
-                    1 for i in range(number_steps)], 
-                unit_conversion_factor=ucf
+        
+        # with length and time interval duration
+        time_interval_duration = 10
+        lengths = [3+i for i in range(two_pipe_system.number_options())]
+        q = two_pipe_system.simplified_heat_transfer_surroundings(
+            specific_heat_transfer_coefficient=u, 
+            temperature_surroundings=ground_temperature, 
+            time_interval_duration=time_interval_duration,
+            length=lengths
+            )
+        assert type(q) == list or type(q) == tuple
+        assert len(q) == two_pipe_system.number_options()
+        for u_i, q_i, l in zip(u, q, lengths):
+            assert math.isclose(
+                q_i, 
+                u_i*(dh_flow_temperature/2+dh_return_temperature/2-ground_temperature)*l*time_interval_duration,
+                abs_tol=1e-3
                 )
-            assert type(sht) == list
-            assert len(sht) == number_steps
-            for _sht in sht:
-                assert isclose(_sht, sht_true*ucf, abs_tol=sht_tol*ucf)
-            
-            sht = srps.simplified_specific_heat_transfer_surroundings(
-                specific_heat_transfer_coefficient=u_coef, 
-                temperature_surroundings=[
-                    temperature_surroundings for i in range(number_steps)], 
-                unit_conversion_factor=ucf
+        # use unit conversion
+        q = two_pipe_system.simplified_heat_transfer_surroundings(
+            specific_heat_transfer_coefficient=u, 
+            temperature_surroundings=ground_temperature, 
+            time_interval_duration=time_interval_duration,
+            length=lengths,
+            unit_conversion_factor=u_factor
+            )
+        assert type(q) == list or type(q) == tuple
+        assert len(q) == two_pipe_system.number_options()
+        for u_i, q_i, l in zip(u, q, lengths):
+            assert math.isclose(
+                q_i, 
+                u_i*(dh_flow_temperature/2+dh_return_temperature/2-ground_temperature)*l*time_interval_duration*u_factor,
+                abs_tol=1e-3
                 )
-            assert type(sht) == list
-            assert len(sht) == number_steps
-            for _sht in sht:
-                assert isclose(_sht, sht_true*ucf, abs_tol=sht_tol*ucf)
-                
-    # *************************************************************************         
     
-    def test_simplified_heat_transfer_surroundings_vector(self):
-        
-        # load pipe data
-            
-        singlepipedata_files = ['tests/data/enerpipe_caldopex_single.csv']
-        singlepipedb = StandardisedPipeDatabase(source=singlepipedata_files)
-        
-        # water 
+    # *************************************************************************
+    # *************************************************************************
         
-        waterdata_file = 'tests/data/incropera2006_saturated_water.csv'
-        phase = FluidDatabase.fluid_LIQUID
-        fluid_db = FluidDatabase(
-            fluid='fluid',
-            phase=phase,
-            source=waterdata_file
-            )
+    def test_simplified_heat_transfer_surroundings_vector_multi_step(self):
         
-        supply_temperature = 80+273.15
-        return_temperature = 50+273.15
-        pressure = 1e5
-        max_specific_pressure_loss = 100 # Pa/m
-        temperature_surroundings = 10+273.15
-        
-        pipe_tuples = [
-            (20,2),     # 1
-            (25,2),     # 2
-            (32,2),     # 3
-            (40,2),     # 4 
-            (50,2),     # 5 
-            (65,2),     # 6 twin pipe too! the DN is all that matters here
-            (80,2),     # 7
-            (100,2),    # 8
-            (125,2),    # 9
-            (150,1)     # 11
-            ]
-    
-        pipes = {
-            pipe_tuple[0:2]: StandardisedPipe(
-                pipe_tuple=pipe_tuple,
-                db=singlepipedb)
-            for pipe_tuple in singlepipedb.pipe_tuples
-            }
+        # 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
 
-        u_coefficients = [
-            0.1392/2,   # 1
-            0.1571/2,   # 2
-            0.1741/2,   # 3
-            0.1662/2,   # 4
-            0.2075/2,   # 5
-            0.2802/2,   # 6
-            0.1747,     # 7
-            0.1992,     # 8
-            0.2771,     # 9 
-            0.3028      # 11
+        # pipe absolute/effective roughness
+        pipe_e_eff = 0.01/1000 # m
+
+        # pipe length
+        pipe_length = 1000
+        
+        # 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',
             ]
-                
-        # specific heat losses extracted graphically from the paper
-        specific_heat_losses = [
-            7.653,      # 1
-            8.6724,     # 2
-            9.5565,     # 3
-            9.1742,     # 4
-            11.3841,    # 5
-            15.4729,    # 6
-            19.2494,    # 7
-            21.9193,    # 8
-            30.4776,    # 9
-            33.3056     # 11
+        
+        # pipe data
+        pipedb = StandardisedPipeDatabase(source=pipedata_files)
+
+        # fluid data   
+        fluiddata_file = 'tests/data/incropera2006_saturated_water.csv'
+        fluiddb = FluidDatabase(fluid='water', phase='l', source=fluiddata_file)
+    
+        # pipe tuples
+        list_pipe_tuples = [
+            pipe_tuple
+            for pipe_tuple in pipedb.pipe_tuples
             ]
         
-        specific_heat_losses_tol = [
-            3.1e-3, # 1 : 0.0030000000000001137
-            0.0320,  # 2 : 0.03190000000000026
-            0.0191,  # 3 : 0.019000000000000128
-            0.0333,  # 4 : 0.033200000000000784
-            0.0284,  # 5 : 0.028399999999999537
-            0.0619, # 6 : 0.06189999999999962
-            0.0325, # 7 : 0.03240000000000265
-            0.0074, # 8 : 0.00730000000000075
-            0.0035, # 9 : 0.003400000000002734
-            0.0025, # 11: 0.002400000000001512
+        # 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
             ]
-        number_options = len(pipe_tuples)
-        srps = SupplyReturnPipeSystem(
-            fluid_db=fluid_db,
-            phase=phase, 
-            pressure=[pressure for i in range(number_options)], 
-            supply_temperature=[
-                supply_temperature for i in range(number_options)
-                ], 
-            return_temperature=[
-                return_temperature for i in range(number_options)
-                ], 
-            supply_pipe=[
-                pipes[pipe_tuple]
-                for pipe_tuple in pipe_tuples
-                ],
-            max_specific_pressure_loss=[
-                max_specific_pressure_loss for i in range(number_options)
-                ]
+        number_options = len(list_pipes)        
+        # two pipe system
+        two_pipe_system = SupplyReturnPipeSystem(
+            fluid_db=fluiddb, 
+            phase=fluiddb.fluid_LIQUID, 
+            pressure=[1e5 for i in range(number_options)],
+            supply_temperature=[dh_flow_temperature for i in range(number_options)], 
+            return_temperature=[dh_return_temperature for i in range(number_options)],  
+            max_specific_pressure_loss=[max_specific_pressure_loss for i in range(number_options)],
+            supply_pipe=list_pipes
             )
-        
-        # simplified heat transfer surroundings
-        sht = srps.simplified_heat_transfer_surroundings(
-            length=[1 for i in range(number_options)], 
-            specific_heat_transfer_coefficient=u_coefficients, 
-            temperature_surroundings=temperature_surroundings, 
-            time_interval_duration=1
+        assert two_pipe_system.vector_mode
+        assert two_pipe_system.number_options() == len(list_pipes)
+        
+        # vector mode, one u value per time interval and option, multiple temperatures
+        number_time_intervals = 3
+        u = [
+            [i+1+k for k in range(number_time_intervals)]
+            for i in range(len(list_pipe_tuples))
+            ]
+        t_surroundings = [ground_temperature+k for k in range(number_time_intervals)]
+        # q should be a list of lists
+        q = two_pipe_system.simplified_specific_heat_transfer_surroundings(
+            specific_heat_transfer_coefficient=u, 
+            temperature_surroundings=t_surroundings, 
             )
-            
-        assert type(sht) == list
-        assert len(sht) == number_options
-        for _sht, sht_true, sht_tol in zip(
-                sht, 
-                specific_heat_losses, 
-                specific_heat_losses_tol
-                ):
-            assert isclose(_sht, sht_true, abs_tol=sht_tol)       
-            
-        sht = srps.simplified_specific_heat_transfer_surroundings(
-            specific_heat_transfer_coefficient=u_coefficients, 
-            temperature_surroundings=temperature_surroundings,
+        assert type(q) == list or type(q) == tuple
+        assert len(q) == two_pipe_system.number_options()
+        for u_i, q_i in zip(u, q):
+            assert type(q_i) == list or type(q_i) == tuple
+            assert len(q_i) == number_time_intervals
+            for k in range(number_time_intervals):
+                assert math.isclose(
+                    q_i[k], 
+                    u_i[k]*(dh_flow_temperature/2+dh_return_temperature/2-t_surroundings[k]),
+                    abs_tol=1e-3
+                    )
+        # use unit conversion
+        u_factor = 2
+        q = two_pipe_system.simplified_specific_heat_transfer_surroundings(
+            specific_heat_transfer_coefficient=u, 
+            temperature_surroundings=t_surroundings, 
+            unit_conversion_factor=u_factor
             )
-            
-        assert type(sht) == list
-        assert len(sht) == number_options
-        for _sht, sht_true, sht_tol in zip(
-                sht, 
-                specific_heat_losses, 
-                specific_heat_losses_tol
-                ):
-            assert isclose(_sht, sht_true, abs_tol=sht_tol)
-            
-        # simplified heat transfer surroundings (temporal vector mode)
-        number_steps = 3
-        sht = srps.simplified_heat_transfer_surroundings(
-            length=[1 for i in range(number_options)], 
-            specific_heat_transfer_coefficient=u_coefficients, 
-            temperature_surroundings=[
-                temperature_surroundings for i in range(number_steps)], 
-            time_interval_duration=[1 for i in range(number_steps)]
+        assert type(q) == list or type(q) == tuple
+        assert len(q) == two_pipe_system.number_options()
+        for u_i, q_i in zip(u, q):
+            assert type(q_i) == list or type(q_i) == tuple
+            assert len(q_i) == number_time_intervals
+            for k in range(number_time_intervals):
+                assert math.isclose(
+                    q_i[k], 
+                    u_i[k]*(dh_flow_temperature/2+dh_return_temperature/2-t_surroundings[k])*u_factor,
+                    abs_tol=1e-3
+                    )
+        
+        # with length and time interval duration
+        time_interval_durations = [10+k for k in range(number_time_intervals)]
+        lengths = [3+i for i in range(two_pipe_system.number_options())]
+        q = two_pipe_system.simplified_heat_transfer_surroundings(
+            specific_heat_transfer_coefficient=u, 
+            temperature_surroundings=t_surroundings, 
+            time_interval_duration=time_interval_durations,
+            length=lengths
             )
-            
-        assert type(sht) == list
-        assert len(sht) == number_options
-        for _sht, sht_true, sht_tol in zip(
-                sht, 
-                specific_heat_losses, 
-                specific_heat_losses_tol
-                ):
-            assert type(_sht) == list
-            assert len(_sht) == number_steps
-            for __sht in _sht:
-                assert isclose(__sht, sht_true, abs_tol=sht_tol)
-                
-        number_steps = 3
-        sht = srps.simplified_specific_heat_transfer_surroundings(
-            specific_heat_transfer_coefficient=u_coefficients, 
-            temperature_surroundings=[
-                temperature_surroundings for i in range(number_steps)]
+        assert type(q) == list or type(q) == tuple
+        assert len(q) == two_pipe_system.number_options()
+        for u_i, q_i, l in zip(u, q, lengths):
+            assert type(q_i) == list or type(q_i) == tuple
+            assert len(q_i) == number_time_intervals
+            for k in range(number_time_intervals):
+                assert math.isclose(
+                    q_i[k], 
+                    u_i[k]*(dh_flow_temperature/2+dh_return_temperature/2-t_surroundings[k])*l*time_interval_durations[k],
+                    abs_tol=1e-3
+                    )
+        # use unit conversion
+        q = two_pipe_system.simplified_heat_transfer_surroundings(
+            specific_heat_transfer_coefficient=u, 
+            temperature_surroundings=t_surroundings, 
+            time_interval_duration=time_interval_durations,
+            length=lengths,
+            unit_conversion_factor=u_factor
             )
-            
-        assert type(sht) == list
-        assert len(sht) == number_options
-        for _sht, sht_true, sht_tol in zip(
-                sht, 
-                specific_heat_losses, 
-                specific_heat_losses_tol
-                ):
-            assert type(_sht) == list
-            assert len(_sht) == number_steps
-            for __sht in _sht:
-                assert isclose(__sht, sht_true, abs_tol=sht_tol)
+        assert type(q) == list or type(q) == tuple
+        assert len(q) == two_pipe_system.number_options()
+        for u_i, q_i, l in zip(u, q, lengths):
+            assert type(q_i) == list or type(q_i) == tuple
+            assert len(q_i) == number_time_intervals
+            for k in range(number_time_intervals):
+                assert math.isclose(
+                    q_i[k], 
+                    u_i[k]*(dh_flow_temperature/2+dh_return_temperature/2-t_surroundings[k])*l*time_interval_durations[k]*u_factor,
+                    abs_tol=1e-3
+                    )
+    
+    # *************************************************************************
+    # *************************************************************************
+    
+    # TODO: cover the main errors
+    
+    def test_simplified_heat_transfer_errors_single_option_single_step(self):
         
-        # try using incorrect option-based inputs
-            
-        error_raised = False
-        try:
-            srps.simplified_heat_transfer_surroundings(
-                length=1, 
-                specific_heat_transfer_coefficient=u_coefficients[0], 
-                temperature_surroundings=temperature_surroundings, 
-                time_interval_duration=1
-                )
-        except TypeError:
-            error_raised = True
-        assert error_raised
-            
-        error_raised = False
-        try:
-            srps.simplified_specific_heat_transfer_surroundings(
-                specific_heat_transfer_coefficient=u_coefficients[0], 
-                temperature_surroundings=temperature_surroundings
-                )
-        except TypeError:
-            error_raised = True
-        assert error_raised
+        pass
         
-        # *********************************************************************
+    def test_simplified_heat_transfer_errors_single_option_multi_step(self):
         
-        # try other units
-        ucf = 2
-        # simplified heat transfer surroundings
-        sht = srps.simplified_heat_transfer_surroundings(
-            length=[1 for i in range(number_options)], 
-            specific_heat_transfer_coefficient=u_coefficients, 
-            temperature_surroundings=temperature_surroundings, 
-            time_interval_duration=1,
-            unit_conversion_factor=ucf
-            )
-            
-        assert type(sht) == list
-        assert len(sht) == number_options
-        for _sht, sht_true, sht_tol in zip(
-                sht, 
-                specific_heat_losses, 
-                specific_heat_losses_tol
-                ):
-            assert isclose(_sht, sht_true*ucf, abs_tol=sht_tol*ucf)
-            
-        sht = srps.simplified_specific_heat_transfer_surroundings(
-            specific_heat_transfer_coefficient=u_coefficients, 
-            temperature_surroundings=temperature_surroundings,
-            unit_conversion_factor=ucf
-            )
-            
-        assert type(sht) == list
-        assert len(sht) == number_options
-        for _sht, sht_true, sht_tol in zip(
-                sht, 
-                specific_heat_losses, 
-                specific_heat_losses_tol
-                ):
-            assert isclose(_sht, sht_true*ucf, abs_tol=sht_tol*ucf)
-            
-        # simplified heat transfer surroundings (temporal vector mode)
-        number_steps = 3
-        sht = srps.simplified_heat_transfer_surroundings(
-            length=[1 for i in range(number_options)], 
-            specific_heat_transfer_coefficient=u_coefficients, 
-            temperature_surroundings=[
-                temperature_surroundings for i in range(number_steps)], 
-            time_interval_duration=[1 for i in range(number_steps)],
-            unit_conversion_factor=ucf
-            )
-            
-        assert type(sht) == list
-        assert len(sht) == number_options
-        for _sht, sht_true, sht_tol in zip(
-                sht, 
-                specific_heat_losses, 
-                specific_heat_losses_tol
-                ):
-            assert type(_sht) == list
-            assert len(_sht) == number_steps
-            for __sht in _sht:
-                assert isclose(__sht, sht_true*ucf, abs_tol=sht_tol*ucf)
+        pass
         
-        number_steps = 3
-        sht = srps.simplified_specific_heat_transfer_surroundings(
-            specific_heat_transfer_coefficient=u_coefficients, 
-            temperature_surroundings=[
-                temperature_surroundings for i in range(number_steps)],
-            unit_conversion_factor=ucf
-            )
-            
-        assert type(sht) == list
-        assert len(sht) == number_options
-        for _sht, sht_true, sht_tol in zip(
-                sht, 
-                specific_heat_losses, 
-                specific_heat_losses_tol
-                ):
-            assert type(_sht) == list
-            assert len(_sht) == number_steps
-            for __sht in _sht:
-                assert isclose(__sht, sht_true*ucf, abs_tol=sht_tol*ucf)
+    def test_simplified_heat_transfer_errors_multi_option_single_step(self):
+        
+        pass
+        
+    def test_simplified_heat_transfer_errors_multi_option_multi_step(self):
+        
+        pass
                 
     # *************************************************************************
     # *************************************************************************
diff --git a/tests/test_trenches.py b/tests/test_trenches.py
index 2ce14ed34657a8bdd580c614f1519ec105368a7f..308cf95fa5a64513b1138c7af6157ffa3a410b6a 100644
--- a/tests/test_trenches.py
+++ b/tests/test_trenches.py
@@ -2710,7 +2710,7 @@ class TestPipeTrench:
     # *************************************************************************
     # *************************************************************************
     
-    def test_simplified_heat_transfer_twin_pipes_single_option(self):
+    def test_simplified_heat_transfer_twin_pipes_single_option_single_step(self):
             
         # maximum specific pressure drop
         max_specific_pressure_loss = 80 # Pa/m
@@ -2808,6 +2808,76 @@ class TestPipeTrench:
             u*(dh_flow_temperature/2+dh_return_temperature/2-ground_temperature)*length*time_interval_duration,
             abs_tol=1e-3
             )
+              
+    # *************************************************************************
+    # *************************************************************************
+    
+    def test_simplified_heat_transfer_twin_pipes_single_option_multi_step(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_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    
+        twin_trench = trenches.TwinPipeTrench(
+            pipe_center_depth=pipe_depth+list_pipes[0].d_cas/2,  
+            fluid_db=fluiddb, 
+            phase=fluiddb.fluid_LIQUID, 
+            pressure=1e5, 
+            supply_temperature=dh_flow_temperature, 
+            return_temperature=dh_return_temperature, 
+            max_specific_pressure_loss=max_specific_pressure_loss, 
+            supply_pipe=list_pipes[0]
+            )
+        assert not twin_trench.vector_mode
+        assert twin_trench.number_options() == 1
         
         # non-vector mode, one u value per time interval (and option), multiple temperatures
         number_time_intervals = 3
@@ -2831,7 +2901,7 @@ class TestPipeTrench:
     # *************************************************************************
     # *************************************************************************
     
-    def test_simplified_heat_transfer_twin_pipes(self):
+    def test_simplified_heat_transfer_twin_pipes_multi_option_single_step(self):
             
         # maximum specific pressure drop
         max_specific_pressure_loss = 80 # Pa/m
@@ -2914,6 +2984,91 @@ class TestPipeTrench:
                 abs_tol=1e-3
                 )
         
+        # with length and time interval duration
+        time_interval_duration = 10
+        lengths = [3+i for i in range(twin_trench.number_options())]
+        q = twin_trench.simplified_heat_transfer_surroundings(
+            specific_heat_transfer_coefficient=u, 
+            temperature_surroundings=ground_temperature, 
+            time_interval_duration=time_interval_duration,
+            length=lengths
+            )
+        assert type(q) == list or type(q) == tuple
+        assert len(q) == twin_trench.number_options()
+        for u_i, q_i, l in zip(u, q, lengths):
+            assert math.isclose(
+                q_i, 
+                u_i*(dh_flow_temperature/2+dh_return_temperature/2-ground_temperature)*l*time_interval_duration,
+                abs_tol=1e-3
+                )
+                    
+    # *************************************************************************
+    # *************************************************************************
+    
+    def test_simplified_heat_transfer_twin_pipes_multi_option_multi_step(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_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    
+        twin_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 twin_trench.vector_mode
+        assert twin_trench.number_options() == len(list_pipes)
+        
         # vector mode, one u value per time interval and option, multiple temperatures
         number_time_intervals = 3
         u = [
@@ -2937,6 +3092,27 @@ class TestPipeTrench:
                     u_i[k]*(dh_flow_temperature/2+dh_return_temperature/2-t_surroundings[k]),
                     abs_tol=1e-3
                     )
+                
+        # with length and time interval duration
+        time_interval_durations = [10+k for k in range(number_time_intervals)]
+        lengths = [3+i for i in range(twin_trench.number_options())]
+        q = twin_trench.simplified_heat_transfer_surroundings(
+            specific_heat_transfer_coefficient=u, 
+            temperature_surroundings=t_surroundings, 
+            time_interval_duration=time_interval_durations,
+            length=lengths
+            )
+        assert type(q) == list or type(q) == tuple
+        assert len(q) == twin_trench.number_options()
+        for u_i, q_i, l in zip(u, q, lengths):
+            assert type(q_i) == list or type(q_i) == tuple
+            assert len(q_i) == number_time_intervals
+            for k in range(number_time_intervals):
+                assert math.isclose(
+                    q_i[k], 
+                    u_i[k]*(dh_flow_temperature/2+dh_return_temperature/2-t_surroundings[k])*l*time_interval_durations[k],
+                    abs_tol=1e-3
+                    )
             
     # *************************************************************************
     # *************************************************************************