diff --git a/src/topupheat/pipes/system.py b/src/topupheat/pipes/system.py
index e4a29756fdb864c6416c56db91071ef48b4b91ad..867ed6dfa9010882cb3c99ca41416d40286b9acf 100644
--- a/src/topupheat/pipes/system.py
+++ b/src/topupheat/pipes/system.py
@@ -558,7 +558,7 @@ class SupplyReturnPipeSystem:
         raise NotImplementedError
         
     # *************************************************************************
-    # TODO: test all 3 modes (single input, list, list of lists,)
+    
     def simplified_specific_heat_transfer_surroundings(
             self,
             specific_heat_transfer_coefficient: float or list,
@@ -620,12 +620,18 @@ class SupplyReturnPipeSystem:
         
         if not self.vector_mode:
             # non-vector mode
-            if type(specific_heat_transfer_coefficient) != type(temperature_surroundings):
-                # case 1 (mismatched input types)
-                raise TypeError('Inconsistent inputs.')
-            elif temporal_vector_mode and len(specific_heat_transfer_coefficient) != len(temperature_surroundings):
-                # case 2 (matching input types but wrong sizes)
-                raise TypeError('Inconsistent inputs.')
+            if temporal_vector_mode:
+                # temporal vector mode
+                if type(specific_heat_transfer_coefficient) != type(temperature_surroundings):
+                    # case 1 (mismatched input types)
+                    raise TypeError('Inconsistent inputs.')
+                if len(specific_heat_transfer_coefficient) != len(temperature_surroundings):
+                    # case 2 (matching input types but wrong sizes)
+                    raise TypeError('Inconsistent inputs.')
+            else:
+                # non-temporal vector mode                
+                if not isinstance(specific_heat_transfer_coefficient, Real):
+                    raise TypeError('Inconsistent inputs.')
         else: 
             # vector mode
             if (type(specific_heat_transfer_coefficient) != tuple and 
@@ -691,10 +697,10 @@ class SupplyReturnPipeSystem:
                 # multiple time steps
                 out = [
                     ((self.supply_temperature+self.return_temperature)*0.5-ts)*
-                    specific_heat_transfer_coefficient*
+                    specific_heat_transfer_coefficient[i]*
                     unit_conversion_factor*
                     (1 if self.losses_are_positive else -1)
-                    for ts in temperature_surroundings
+                    for i, ts in enumerate(temperature_surroundings)
                     ]
             else:
                 # one time step
@@ -708,7 +714,7 @@ class SupplyReturnPipeSystem:
         return out
     
     # *************************************************************************
-    # TODO: implement using simplified_specific_heat_transfer_surroundings method
+    
     def simplified_heat_transfer_surroundings(
             self,
             length: float or list,
@@ -717,140 +723,88 @@ class SupplyReturnPipeSystem:
             time_interval_duration: float or list,
             unit_conversion_factor: float = 1.0
             ):
-        raise NotImplementedError
-        # """
-        # Calculates the heat transfer with the surroundings using a constant
-        # specific heat transfer coefficient.
+        """
+        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.
+        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.
+        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.
+        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
-        # # temperature_surroundings >> temporal
-        # # determine if the inputs are for temporal vector/normal mode
-        # if (isinstance(temperature_surroundings, Real) and
-        #     isinstance(time_interval_duration, Real)):
-        #     temporal_vector_mode = False
-        # elif ((type(temperature_surroundings) == list or 
-        #        type(temperature_surroundings) == tuple) and
-        #       (type(time_interval_duration) == list or 
-        #        type(time_interval_duration) == tuple)):
-        #     temporal_vector_mode = True
-        # else:
-        #     raise TypeError('Incompatible inputs.')
-        
-        # # confirm inputs are consistent with vector/normal mode
-        # if ((not self.vector_mode and 
-        #      (not isinstance(length, Real) or 
-        #       not isinstance(specific_heat_transfer_coefficient, Real))) or
-        #     (self.vector_mode and 
-        #      (type(length) != list and 
-        #       type(length) != tuple or
-        #       len(length) != self.number_options() or
-        #       len(specific_heat_transfer_coefficient) != len(length))
-        #      )
-        #     ):
-        #     raise TypeError('Inconsistent inputs.')
+        """
         
-        # if self.vector_mode:
-        #     # vector mode: multiple options
-        #     if temporal_vector_mode:
-        #         out = [
-        #             [((st+rt)*0.5-ts)*
-        #              shtc*
-        #              pipe.length*
-        #              dt*
-        #              unit_conversion_factor*
-        #              (1 if self.losses_are_positive else -1)*
-        #              (1 if self.twin_pipes else 2)
-        #              for ts, dt in zip(
-        #                      temperature_surroundings, 
-        #                      time_interval_duration
-        #                      )
-        #              ]
-        #             for st, rt, pipe, shtc in zip(
-        #                     self.supply_temperature, 
-        #                     self.return_temperature, 
-        #                     self.supply_pipe,
-        #                     specific_heat_transfer_coefficient
-        #                     )
-        #             ]
-        #     else:
-        #         # one time step
-        #         out = [
-        #             ((st+rt)*0.5-temperature_surroundings)*
-        #             shtc*
-        #             pipe.length*
-        #             time_interval_duration*
-        #             unit_conversion_factor*
-        #             (1 if self.losses_are_positive else -1)*
-        #             (1 if self.twin_pipes else 2)
-        #             for st, rt, pipe, shtc in zip(
-        #                     self.supply_temperature, 
-        #                     self.return_temperature, 
-        #                     self.supply_pipe,
-        #                     specific_heat_transfer_coefficient
-        #                     )
-        #             ]
-        # else:
-        #     # normal mode: one option
-        #     if temporal_vector_mode:
-        #         # multiple time steps
-        #         out = [
-        #             ((self.supply_temperature+self.return_temperature)*0.5-ts)*
-        #             specific_heat_transfer_coefficient*
-        #             length*
-        #             dt*
-        #             unit_conversion_factor*
-        #             (1 if self.losses_are_positive else -1)*
-        #             (1 if self.twin_pipes else 2)
-        #             for ts, dt in zip(
-        #                     temperature_surroundings, 
-        #                     time_interval_duration
-        #                     )
-        #             ]
-        #     else:
-        #         # one time step
-        #         out = (
-        #             ((self.supply_temperature+self.return_temperature)*0.5-
-        #              temperature_surroundings)*
-        #             specific_heat_transfer_coefficient*
-        #             length*
-        #             time_interval_duration*
-        #             unit_conversion_factor*
-        #             (1 if self.losses_are_positive else -1)*
-        #             (1 if self.twin_pipes else 2)
-        #             )
-        # return out
+        # length has to match the number of options
+        if self.vector_mode:
+            if type(length) != tuple and type(length) != list:
+                # not an iterable
+                raise TypeError('Inconsistent inputs.')
+        else: # non-vector mode
+            if not isinstance(length, Real):
+                # not a number
+                raise TypeError('Inconsistent inputs.')
+        # time interval durations have to match the number of intervals
+        if (type(temperature_surroundings) != type (time_interval_duration) or
+            len(temperature_surroundings) != len(time_interval_duration)):
+            if not isinstance(temperature_surroundings, Real) and not isinstance(time_interval_duration, Real):
+                raise TypeError('Inconsistent inputs.')
+        # compute the specific heat transfer
+        qs = self.simplified_specific_heat_transfer_surroundings(
+            specific_heat_transfer_coefficient,
+            temperature_surroundings,
+            unit_conversion_factor=unit_conversion_factor
+            )
+        # multiply by time interval duration (interval) and length (option)
+        if self.vector_mode:
+            # multiple options
+            if isinstance(time_interval_duration, Real):
+                # single time interval
+                return [
+                    _qs*_l*time_interval_duration
+                    for _qs, _l in zip(qs, length)
+                    ]
+            else:
+                # multiple time intervals
+                return [
+                    [_qs_i*_l*dt 
+                     for _qs_i, dt in zip(_qs, time_interval_duration)]
+                    for _qs, _l in zip(qs, length)
+                    ]
+        else:
+            # single option
+            if isinstance(time_interval_duration, Real):
+                # single time interval
+                return qs*length*time_interval_duration
+            else:
+                # multiple time intervals
+                return [
+                    _qs_i*length*dt 
+                    for _qs_i, dt in enumerate(qs, time_interval_duration)
+                    ]
     
     # *************************************************************************
     
diff --git a/tests/test_system.py b/tests/test_system.py
index 284036fc477a714b723b6a151693c30a3db77a5c..138ef4c17cd4f28f7b03086abceb38fee8c2465c 100644
--- a/tests/test_system.py
+++ b/tests/test_system.py
@@ -12,6 +12,25 @@ from math import isclose, inf, ceil
 
 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):
+        
+        pass
+        
+    def test_simplified_heat_transfer_surroundings_vector_multi_step(self):
+        
+        pass
+    
+    
     def test_simplified_heat_transfer_surroundings(self):
         
         # load pipe data
diff --git a/tests/test_trenches.py b/tests/test_trenches.py
index 0fe0fd309fa2c2b4ebccbaaa3d6389d3c50f5ce9..2ce14ed34657a8bdd580c614f1519ec105368a7f 100644
--- a/tests/test_trenches.py
+++ b/tests/test_trenches.py
@@ -1,6 +1,6 @@
 import math
 import numpy as np
-# from numbers import Real
+from numbers import Real
 from src.topupheat.common.fluids import FluidDatabase, Fluid
 from src.topupheat.pipes import trenches, fic
 from src.topupheat.pipes.single import InsulatedPipe, StandardisedPipeDatabase
@@ -2705,6 +2705,128 @@ class TestPipeTrench:
             #     )
             
             # assert type(trench.printable_description()) == str
+            
+                    
+    # *************************************************************************
+    # *************************************************************************
+    
+    def test_simplified_heat_transfer_twin_pipes_single_option(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 option, and one temperature
+        u = 1
+        # q should be a number
+        q = twin_trench.simplified_specific_heat_transfer_surroundings(
+            specific_heat_transfer_coefficient=u, 
+            temperature_surroundings=ground_temperature, 
+            )
+        assert isinstance(q, Real)
+        assert 1 == twin_trench.number_options()
+        assert math.isclose(
+            q, 
+            u*(dh_flow_temperature/2+dh_return_temperature/2-ground_temperature),
+            abs_tol=1e-3
+            )
+        
+        # with length and time interval duration
+        length = 3
+        time_interval_duration = 10
+        q = twin_trench.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 == twin_trench.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
+            )
+        
+        # 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 = twin_trench.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
+                )
                     
     # *************************************************************************
     # *************************************************************************