Repository: MassimoCimmino/pygfunction (language: python)

Issue:
Negative mass flow rate in pipe and network objects to model reversed flow
This is to implement the use of negative mass flow rates when calling `Pipe` and `Network` objects to signify reversed flow direction. This will support #215 and enable a more intuitive modelling of borehole thermal energy storage systems.


Discussion / hints:
Thanks for your work on reversible flow! The BTES often uses different directions of fluid in winter and summer, so it would be crucial to support this.

Is there any prototype code that can be tested?
@metab0t This issue is only for the impact of the flow reversal on the interior of boreholes. The impact on the ground heat transfer will be tackled in #215.

The present issue and #215 are my priorities at the moment. I will update #215 as development progresses.

Fix patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index dba5f28..f74c912 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -2,6 +2,10 @@
 
 ## Version 2.3 (in development)
 
+### New features
+
+* [Issue 282](https://github.com/MassimoCimmino/pygfunction/issues/282) - Enabled the use of negative mass flow rates in `Pipe` and `Network` classes to model reversed flow direction.
+
 ## Version 2.2.3 (2024-07-01)
 
 ### New features
diff --git a/pygfunction/gfunction.py b/pygfunction/gfunction.py
index 2a43154..abb73d6 100644
--- a/pygfunction/gfunction.py
+++ b/pygfunction/gfunction.py
@@ -1677,8 +1677,10 @@ def solve(self, time, alpha):
                     # The gFunction is equal to the effective borehole wall
                     # temperature
                     # Outlet fluid temperature
-                    T_f_out = T_f_in - 2*pi*self.network.p[0].k_s*H_tot/(
-                        np.sum(self.network.m_flow_network*self.network.cp_f))
+                    T_f_out = T_f_in - 2*pi*self.network.p[0].k_s*H_tot / (
+                        np.sum(
+                            np.abs(self.network.m_flow_network)
+                            * self.network.cp_f))
                     # Average fluid temperature
                     T_f = 0.5*(T_f_in + T_f_out)
                     # Borefield thermal resistance
diff --git a/pygfunction/networks.py b/pygfunction/networks.py
index f4ba8c7..af37aec 100644
--- a/pygfunction/networks.py
+++ b/pygfunction/networks.py
@@ -12,7 +12,8 @@ class Network(object):
     characteristics of the pipes and the grout material in each boreholes, the
     topology of the connections between boreholes, as well as methods to
     evaluate fluid temperatures and heat extraction rates based on the work of
-    Cimmino (2018, 2019) [#Network-Cimmin2018]_, [#Network-Cimmin2019]_.
+    Cimmino (2018, 2019, 2024) [#Network-Cimmin2018]_, [#Network-Cimmin2019]_,
+    [#Network-Cimmino2024]_.
 
     Attributes
     ----------
@@ -65,6 +66,9 @@ class Network(object):
        g-function calculation of bore fields with series- and
        parallel-connected boreholes. Science and Technology for the Built
        Environment, 25 (8), 1007-1022.
+    .. [#Network-Cimmino2024] Cimmino, M. (2024). g-Functions for fields of
+       series- and parallel-connected boreholes with variable fluid mass flow
+       rate and reversible flow direction. Renewable Energy, 228, 120661.
 
     """
     def __init__(self, boreholes, pipes, bore_connectivity=None,
@@ -81,18 +85,8 @@ def __init__(self, boreholes, pipes, bore_connectivity=None,
         self.cp_f = cp_f
 
         # Verify that borehole connectivity is valid
-        _verify_bore_connectivity(bore_connectivity, self.nBoreholes)
-        iInlets, nInlets, iOutlets, nOutlets, iCircuit = _find_inlets_outlets(
-                bore_connectivity, self.nBoreholes)
-
-        # Number of inlets and outlets in network
-        self.nInlets = nInlets
-        self.nOutlets = nOutlets
-        # Indices of inlets and outlets in network
-        self.iInlets = iInlets
-        self.iOutlets = iOutlets
-        # Indices of circuit of each borehole in network
-        self.iCircuit = iCircuit
+        self._verify_bore_connectivity()
+        self._find_inlets_outlets()
 
         # Initialize stored_coefficients
         self._initialize_coefficients_connectivity()
@@ -140,7 +134,7 @@ def get_inlet_temperature(
                 m_flow_network, cp_f, nSegments, segment_ratios=segment_ratios)
         # Evaluate outlet temperatures
         if np.isscalar(T_b):
-            T_b = np.tile(T_b, sum(self.nSegments))
+            T_b = np.tile(T_b, np.size(a_b, axis=1))
         T_f_in_borehole = a_in @ np.atleast_1d(T_f_in) + a_b @ T_b
         return T_f_in_borehole
 
@@ -185,7 +179,7 @@ def get_outlet_temperature(
                 m_flow_network, cp_f, nSegments, segment_ratios=segment_ratios)
         # Evaluate outlet temperatures
         if np.isscalar(T_b):
-            T_b = np.tile(T_b, sum(self.nSegments))
+            T_b = np.tile(T_b, np.size(a_b, axis=1))
         T_f_out = a_in @ np.atleast_1d(T_f_in) + a_b @ T_b
         return T_f_out
 
@@ -228,7 +222,7 @@ def get_borehole_heat_extraction_rate(
         a_in, a_b = self.coefficients_borehole_heat_extraction_rate(
                 m_flow_network, cp_f, nSegments, segment_ratios=segment_ratios)
         if np.isscalar(T_b):
-            T_b = np.tile(T_b, sum(self.nSegments))
+            T_b = np.tile(T_b, np.size(a_b, axis=1))
         Q_b = a_in @ np.atleast_1d(T_f_in) + a_b @ T_b
 
         return Q_b
@@ -272,7 +266,7 @@ def get_fluid_heat_extraction_rate(
         a_in, a_b = self.coefficients_fluid_heat_extraction_rate(
                 m_flow_network, cp_f, nSegments, segment_ratios=segment_ratios)
         if np.isscalar(T_b):
-            T_b = np.tile(T_b, sum(self.nSegments))
+            T_b = np.tile(T_b, np.size(a_b, axis=1))
         Q_f = a_in @ np.atleast_1d(T_f_in) + a_b @ T_b
 
         return Q_f
@@ -319,7 +313,7 @@ def get_network_inlet_temperature(
                 m_flow_network, cp_f, nSegments, segment_ratios=segment_ratios)
         # Evaluate outlet temperatures
         if np.isscalar(T_b):
-            T_b = np.tile(T_b, sum(self.nSegments))
+            T_b = np.tile(T_b, np.size(a_b, axis=1))
         T_f_in = a_q @ np.atleast_1d(Q_t) + a_b @ T_b
         if np.isscalar(Q_t):
             T_f_in = T_f_in.item()
@@ -367,7 +361,7 @@ def get_network_outlet_temperature(
                 m_flow_network, cp_f, nSegments, segment_ratios=segment_ratios)
         # Evaluate outlet temperatures
         if np.isscalar(T_b):
-            T_b = np.tile(T_b, sum(self.nSegments))
+            T_b = np.tile(T_b, np.size(a_b, axis=1))
         T_f_out = a_in @ np.atleast_1d(T_f_in) + a_b @ T_b
         if np.isscalar(T_f_in):
             T_f_out = T_f_out.item()
@@ -413,7 +407,7 @@ def get_network_heat_extraction_rate(
         a_in, a_b = self.coefficients_network_heat_extraction_rate(
                 m_flow_network, cp_f, nSegments, segment_ratios=segment_ratios)
         if np.isscalar(T_b):
-            T_b = np.tile(T_b, sum(self.nSegments))
+            T_b = np.tile(T_b, np.size(a_b, axis=1))
         Q_t = a_in @ np.atleast_1d(T_f_in) + a_b @ T_b
         if np.isscalar(T_f_in):
             Q_t = Q_t.item()
@@ -960,7 +954,8 @@ def _coefficients_mixing(self, m_flow_network):
         """
         if not self._check_mixing_coefficients(m_flow_network):
             self._mix_out = np.zeros((1, self.nBoreholes))
-            self._mix_out[0, self.iOutlets] = self._m_flow_in/np.sum(self._m_flow_in)
+            self._mix_out[0, self.iOutlets] = \
+                np.abs(self._m_flow_in) / np.sum(np.abs(self._m_flow_in))
             self._mixing_m_flow = m_flow_network
         return self._mix_out
 
@@ -975,13 +970,22 @@ def _initialize_coefficients_connectivity(self):
                 + \\mathbf{c_{out}} \\mathbf{T_{f,borehole,out}}
 
         """
-        self._c_in = np.zeros((self.nBoreholes, 1))
-        self._c_out = np.zeros((self.nBoreholes, self.nBoreholes))
+        # Flow in positive direction
+        self._c_in_pos = np.zeros((self.nBoreholes, 1))
+        self._c_out_pos = np.zeros((self.nBoreholes, self.nBoreholes))
+        for i in range(self.nInlets):
+            self._c_in_pos[self._iInlets_pos[i], 0] = 1.
+        for i in range(self.nBoreholes):
+            if not self.c[i] == -1:
+                self._c_out_pos[i, self.c[i]] = 1.
+        # Flow in negative direction
+        self._c_in_neg = np.zeros((self.nBoreholes, 1))
+        self._c_out_neg = np.zeros((self.nBoreholes, self.nBoreholes))
         for i in range(self.nInlets):
-            self._c_in[self.iInlets[i], 0] = 1.
+            self._c_in_neg[self._iInlets_neg[i], 0] = 1.
         for i in range(self.nBoreholes):
             if not self.c[i] == -1:
-                self._c_out[i, self.c[i]] = 1.
+                self._c_out_neg[self.c[i], i] = 1.
 
         return
 
@@ -1078,6 +1082,22 @@ def _format_inputs(self, m_flow_network, cp_f, nSegments, segment_ratios):
             raise ValueError(
                 'Incorrect length of mass flow vector.')
         self._m_flow_in = m_flow_in
+        # Flow direction
+        if np.all(m_flow_network >= 0.):
+            self._is_reversed = False
+            self._c_in = self._c_in_pos
+            self._c_out = self._c_out_pos
+            self.iInlets = self._iInlets_pos
+            self.iOutlets = self._iOutlets_pos
+        elif np.all(m_flow_network <= 0.):
+            self._is_reversed = True
+            self._c_in = self._c_in_neg
+            self._c_out = self._c_out_neg
+            self.iInlets = self._iInlets_neg
+            self.iOutlets = self._iOutlets_neg
+        else:
+            raise ValueError(
+                'All elements of m_flow_network should be of the same sign.')
 
         # Format heat capacity inputs
         # Heat capacity in each fluid circuit
@@ -1121,6 +1141,115 @@ def _format_inputs(self, m_flow_network, cp_f, nSegments, segment_ratios):
                 'Incorrect format of the segment ratios list.')
 
 
+    def _find_inlets_outlets(self):
+        """
+        Finds the numbers of boreholes connected to the inlet and outlet of the
+        network and the indices of the boreholes.
+
+        This function raises an error if the supplied borehole connectivity is
+        invalid.
+
+        Parameters
+        ----------
+        bore_connectivity : list
+            Index of fluid inlet into each borehole. -1 corresponds to a borehole
+            connected to the bore field inlet.
+        nBoreholes : int
+            Number of boreholes in the bore field.
+
+        """
+        # Number and indices of inlets
+        nInlets = self.c.count(-1)
+        iInlets = [i for i in range(self.nBoreholes) if self.c[i] == -1]
+        # Number and indices of outlets
+        iOutlets = [i for i in range(self.nBoreholes) if i not in self.c]
+        nOutlets = len(iOutlets)
+        iCircuit = [iInlets.index(self._path_to_inlet(i)[-1])
+                    for i in range(self.nBoreholes)]
+        if not nInlets == nOutlets:
+            raise ValueError(
+                'The network should have as many inlets as outlets.')
+
+        # Number of inlets and outlets in network
+        self.nInlets = nInlets
+        self.nOutlets = nOutlets
+        # Indices of inlets and outlets in network
+        self._iInlets_pos = iInlets
+        self._iOutlets_pos = iOutlets
+        self._iInlets_neg = iOutlets
+        self._iOutlets_neg = iInlets
+        # Indices of circuit of each borehole in network
+        self.iCircuit = iCircuit
+
+        return
+
+
+    def _path_to_inlet(self, bore_index):
+        """
+        Returns the path from a borehole to the bore field inlet.
+
+        Parameters
+        ----------
+        bore_index : int
+            Index of borehole to evaluate path.
+
+        Returns
+        -------
+        path : list
+            List of boreholes leading to the bore field inlet, starting from
+            borehole bore_index
+
+        """
+        # Initialize path
+        path = [bore_index]
+        # Index of borehole feeding into borehole (bore_index)
+        index_in = self.c[bore_index]
+        # Stop when bore field inlet is reached (index_in == -1)
+        while not index_in == -1:
+            # Add index of upstream borehole to path
+            path.append(index_in)
+            # Get index of next upstream borehole
+            index_in = self.c[index_in]
+
+        return path
+
+
+    def _verify_bore_connectivity(self):
+        """
+        Verifies that borehole connectivity is valid.
+
+        This function raises an error if the supplied borehole connectivity is
+        invalid.
+
+        """
+        if not len(self.c) == self.nBoreholes:
+            raise ValueError(
+                'The length of the borehole connectivity list does not correspond '
+                'to the number of boreholes in the bore field.')
+        if max(self.c) >= self.nBoreholes:
+            raise ValueError(
+                'The borehole connectivity list contains borehole indices that '
+                'are not part of the network.')
+        # Cycle through each borehole and verify that connections lead to -1
+        # (-1 is the bore field inlet) and that no two boreholes have the same
+        # index of fluid inlet (except for -1).
+        for index_in in self.c:
+            n = 0 # Initialize step counter
+            if index_in != -1 and self.c.count(index_in) > 1:
+                raise ValueError(
+                    'Two boreholes cannot have the same inlet, except fort the '
+                    'network inlet (index of -1).')
+            # Stop when bore field inlet is reached (index_in == -1)
+            while not index_in == -1:
+                index_in = self.c[index_in]
+                n += 1 # Increment step counter
+                # Raise error if n exceeds the number of boreholes
+                if n > self.nBoreholes:
+                    raise ValueError(
+                        'The borehole connectivity list is invalid.')
+        return
+
+
 class _EquivalentNetwork(Network):
     """
     Class for networks of equivalent boreholes with parallel connections
@@ -1130,7 +1259,8 @@ class _EquivalentNetwork(Network):
     characteristics of the pipes and the grout material in each boreholes, the
     topology of the connections between boreholes, as well as methods to
     evaluate fluid temperatures and heat extraction rates based on the work of
-    Cimmino (2018, 2019) [#Network-Cimmin2018]_, [#Network-Cimmin2019]_.
+    Cimmino (2018, 2019, 2024) [#Network-Cimmin2018]_, [#Network-Cimmin2019]_,
+    [#Network-Cimmin2024]_.
 
     Attributes
     ----------
@@ -1178,6 +1308,9 @@ class _EquivalentNetwork(Network):
        g-function calculation of bore fields with series- and
        parallel-connected boreholes. Science and Technology for the Built
        Environment, 25 (8), 1007-1022.
+    .. [#Network-Cimmino2024] Cimmino, M. (2024). g-Functions for fields of
+       series- and parallel-connected boreholes with variable fluid mass flow
+       rate and reversible flow direction. Renewable Energy, 228, 120661.
 
     """
     def __init__(self, equivalentBoreholes, pipes, m_flow_network=None,
@@ -1194,17 +1327,7 @@ def __init__(self, equivalentBoreholes, pipes, m_flow_network=None,
         self.cp_f = cp_f
 
         # Verify that borehole connectivity is valid
-        iInlets, nInlets, iOutlets, nOutlets, iCircuit = _find_inlets_outlets(
-                self.c, self.nBoreholes)
-
-        # Number of inlets and outlets in network
-        self.nInlets = nInlets
-        self.nOutlets = nOutlets
-        # Indices of inlets and outlets in network
-        self.iInlets = iInlets
-        self.iOutlets = iOutlets
-        # Indices of circuit of each borehole in network
-        self.iCircuit = iCircuit
+        self._find_inlets_outlets()
 
         # Initialize stored_coefficients
         self._initialize_coefficients_connectivity()
@@ -1303,9 +1426,10 @@ def _coefficients_mixing(self, m_flow_network):
         """
         if not self._check_mixing_coefficients(m_flow_network):
             self._mix_out = np.zeros((1, self.nBoreholes))
+            m_flow_in = np.abs(self._m_flow_in)
+            wBoreholes = self.wBoreholes.flatten()
             self._mix_out[0, self.iOutlets] = \
-                self._m_flow_in * self.wBoreholes.flatten() \
-                / np.sum(self._m_flow_in * self.wBoreholes.flatten())
+                m_flow_in * wBoreholes / np.sum(m_flow_in * wBoreholes)
             self._mixing_m_flow = m_flow_network
         return self._mix_out
 
@@ -1323,6 +1447,18 @@ def _format_inputs(self, m_flow_network, cp_f, nSegments, segment_ratios):
             raise ValueError(
                 'Incorrect length of mass flow vector.')
         self._m_flow_in = m_flow_in
+        # Flow direction
+        self._is_reversed = m_flow_network < 0.
+        if self._is_reversed:
+            self._c_in = self._c_in_neg
+            self._c_out = self._c_out_neg
+            self.iInlets = self._iOutlets_neg
+            self.iOutlets = self._iInlets_neg
+        else:
+            self._c_in = self._c_in_pos
+            self._c_out = self._c_out_pos
+            self.iInlets = self._iOutlets_pos
+            self.iOutlets = self._iInlets_pos
 
         # Format heat capacity inputs
         # Heat capacity in each fluid circuit
@@ -1410,112 +1546,3 @@ def network_thermal_resistance(network, m_flow_network, cp_f):
         R_field = R_field.item()
 
     return R_field
-
-
-def _find_inlets_outlets(bore_connectivity, nBoreholes):
-    """
-    Finds the numbers of boreholes connected to the inlet and outlet of the
-    network and the indices of the boreholes.
-
-    This function raises an error if the supplied borehole connectivity is
-    invalid.
-
-    Parameters
-    ----------
-    bore_connectivity : list
-        Index of fluid inlet into each borehole. -1 corresponds to a borehole
-        connected to the bore field inlet.
-    nBoreholes : int
-        Number of boreholes in the bore field.
-
-    """
-    # Number and indices of inlets
-    nInlets = bore_connectivity.count(-1)
-    iInlets = [i for i in range(nBoreholes) if bore_connectivity[i]==-1]
-    # Number and indices of outlets
-    iOutlets = [i for i in range(nBoreholes) if i not in bore_connectivity]
-    nOutlets = len(iOutlets)
-    iCircuit = [iInlets.index(_path_to_inlet(bore_connectivity, i)[-1])
-                for i in range(nBoreholes)]
-    if not nInlets == nOutlets:
-        raise ValueError(
-            'The network should have as many inlets as outlets.')
-
-    return iInlets, nInlets, iOutlets, nOutlets, iCircuit
-
-
-def _path_to_inlet(bore_connectivity, bore_index):
-    """
-    Returns the path from a borehole to the bore field inlet.
-
-    Parameters
-    ----------
-    bore_connectivity : list
-        Index of fluid inlet into each borehole. -1 corresponds to a borehole
-        connected to the bore field inlet.
-    bore_index : int
-        Index of borehole to evaluate path.
-
-    Returns
-    -------
-    path : list
-        List of boreholes leading to the bore field inlet, starting from
-        borehole bore_index
-
-    """
-    # Initialize path
-    path = [bore_index]
-    # Index of borehole feeding into borehole (bore_index)
-    index_in = bore_connectivity[bore_index]
-    # Stop when bore field inlet is reached (index_in == -1)
-    while not index_in == -1:
-        # Add index of upstream borehole to path
-        path.append(index_in)
-        # Get index of next upstream borehole
-        index_in = bore_connectivity[index_in]
-
-    return path
-
-
-def _verify_bore_connectivity(bore_connectivity, nBoreholes):
-    """
-    Verifies that borehole connectivity is valid.
-
-    This function raises an error if the supplied borehole connectivity is
-    invalid.
-
-    Parameters
-    ----------
-    bore_connectivity : list
-        Index of fluid inlet into each borehole. -1 corresponds to a borehole
-        connected to the bore field inlet.
-    nBoreholes : int
-        Number of boreholes in the bore field.
-
-    """
-    if not len(bore_connectivity) == nBoreholes:
-        raise ValueError(
-            'The length of the borehole connectivity list does not correspond '
-            'to the number of boreholes in the bore field.')
-    if max(bore_connectivity) >= nBoreholes:
-        raise ValueError(
-            'The borehole connectivity list contains borehole indices that '
-            'are not part of the network.')
-    # Cycle through each borehole and verify that connections lead to -1
-    # (-1 is the bore field inlet) and that no two boreholes have the same
-    # index of fluid inlet (except for -1).
-    for index_in in bore_connectivity:
-        n = 0 # Initialize step counter
-        if index_in != -1 and bore_connectivity.count(index_in) > 1:
-            raise ValueError(
-                'Two boreholes cannot have the same inlet, except fort the '
-                'network inlet (index of -1).')
-        # Stop when bore field inlet is reached (index_in == -1)
-        while not index_in == -1:
-            index_in = bore_connectivity[index_in]
-            n += 1 # Increment step counter
-            # Raise error if n exceeds the number of boreholes
-            if n > nBoreholes:
-                raise ValueError(
-                    'The borehole connectivity list is invalid.')
-    return
diff --git a/pygfunction/pipes.py b/pygfunction/pipes.py
index f74f12f..0eb028e 100644
--- a/pygfunction/pipes.py
+++ b/pygfunction/pipes.py
@@ -24,6 +24,11 @@ class _BasePipe(object):
         Total number of pipe inlets, equals to 1.
     nOutlets : int
         Total number of pipe outlets, equals to 1.
+    reversible_flow : bool
+        True to treat a negative mass flow rate as the reversal of flow
+        direction within the borehole. If False, the direction of flow is not
+        reversed when the mass flow rate is negative, and the absolute value is
+        used for calculations. Equals to True.
 
     Notes
     -----
@@ -41,6 +46,7 @@ def __init__(self, borehole):
         self.nPipes = 1
         self.nInlets = 1
         self.nOutlets = 1
+        self.reversible_flow = True
 
     def get_temperature(
             self, z, T_f_in, T_b, m_flow_borehole, cp_f, segment_ratios=None):
@@ -405,19 +411,19 @@ def coefficients_outlet_temperature(
                 self.b.H, m_flow_borehole, cp_f, nSegments,
                 segment_ratios=segment_ratios)
 
+            # Coefficients for the fluid temperatures
+            # [T_f](z) = [f_fd]*[T_fd](z) + [f_fu][T_fu](z)
+            f_fd, f_fu = self._flow_direction(
+                m_flow_borehole, cp_f, nSegments,
+                segment_ratios=segment_ratios)
+
             # Final coefficient matrices for outlet temperatures:
             # [T_{f,out}] = [a_in]*[T_{f,in}] + [a_b]*[T_b]
             ImI = np.hstack(
                 (np.eye(self.nPipes),
                  -np.eye(self.nPipes)))
-            CfuoI = np.vstack(
-                (c_fu,
-                 np.eye(self.nPipes)))
-            CinoZ = np.vstack(
-                (c_in,
-                 np.zeros((self.nPipes, self.nInlets))))
-            A = d_fu @ np.linalg.solve(ImI @ e_f0 @ CfuoI, -ImI)
-            a_in = A @ e_f0 @ CinoZ
+            A = d_fu @ np.linalg.solve(ImI @ e_f0 @ (f_fu + f_fd @ c_fu), -ImI)
+            a_in = A @ e_f0 @ f_fd @ c_in
             a_b = A @ e_b
 
             # Store coefficients
@@ -501,19 +507,20 @@ def coefficients_temperature(
         f_f0, f_b = self._general_solution(
             z, m_flow_borehole, cp_f, nSegments, segment_ratios=segment_ratios)
 
+        # Coefficients for the fluid temperatures
+        # [T_f](z) = [g_fd]*[T_fd](z) + [g_fu][T_fu](z)
+        g_fd, g_fu = self._flow_direction(
+            m_flow_borehole, cp_f, nSegments,
+            segment_ratios=segment_ratios)
+
         # Final coefficient matrices for outlet temperatures:
         # [T_{f,out}] = [a_in]*[T_{f,in}] + [a_b]*[T_b]
         ImI = np.hstack(
             (np.eye(self.nPipes),
              -np.eye(self.nPipes)))
-        CfuoI = np.vstack(
-            (c_fu,
-             np.eye(self.nPipes)))
-        CinoZ = np.vstack(
-            (c_in,
-             np.zeros((self.nPipes, self.nInlets))))
-        A = CfuoI @ np.linalg.solve(ImI @ e_f0 @ CfuoI, -ImI)
-        a_in = f_f0 @ (np.eye(2*self.nPipes) + A @ e_f0) @ CinoZ
+        B = (g_fu + g_fd @ c_fu)
+        A = B @ np.linalg.solve(ImI @ e_f0 @ B, -ImI)
+        a_in = f_f0 @ (np.eye(2*self.nPipes) + A @ e_f0) @ g_fd @ c_in
         a_b = (f_f0 @ A @ e_b) + f_b
 
         return a_in, a_b
@@ -566,8 +573,7 @@ def coefficients_borehole_heat_extraction_rate(
                 m_flow_borehole, cp_f, nSegments, segment_ratios)
             m_flow_pipe = self._m_flow_pipe
             cp_pipe = self._cp_pipe
-            mcp = np.hstack((-m_flow_pipe[0:nPipes],
-                             m_flow_pipe[-nPipes:]))*cp_pipe
+            mcp = m_flow_pipe * cp_pipe
 
             # Initialize coefficient matrices
             a_in = np.zeros((nSegments, self.nInlets))
@@ -578,8 +584,8 @@ def coefficients_borehole_heat_extraction_rate(
             aTf, bTf = self.coefficients_temperature(
                 z, m_flow_borehole, cp_f, nSegments,
                 segment_ratios=segment_ratios)
-            a_in = mcp @ (aTf[:-1,:,:] - aTf[1:,:,:])
-            a_b = mcp @ (bTf[:-1,:,:] - bTf[1:,:,:])
+            a_in = mcp @ (aTf[1:,:,:] - aTf[:-1,:,:])
+            a_b = mcp @ (bTf[1:,:,:] - bTf[:-1,:,:])
 
             # Store coefficients
             self._set_stored_coefficients(
@@ -641,7 +647,7 @@ def coefficients_fluid_heat_extraction_rate(
 
             # Intermediate matrices for fluid heat extraction rates:
             # [Q_{f}] = [c_in]*[T_{f,in}] + [c_out]*[T_{f,out}]
-            MCP = self._m_flow_in * self._cp_in
+            MCP = np.abs(self._m_flow_in) * self._cp_in
             c_in = -np.diag(MCP)
             c_out = np.diag(MCP)
 
@@ -927,20 +933,20 @@ def _check_geometry(self):
 
         return True
 
-    def _pipe_connectivity(
+    def _continuity_condition_head(
             self, m_flow_borehole, cp_f, nSegments, segment_ratios=None):
         """ Returns coefficients for the relation
-            [T_{f,out}] = [c_fu]*[T_fu](z=0)
+            [T_fd](z=0) = [c_in]*[T_{f,in}] + [c_fu]*[T_fu](z=0)
         """
         raise NotImplementedError(
             '_continuity_condition_head class method not implemented, '
             'this method should return matrices for the relation: '
-            '[T_{f,out}] = [c_fu]*[T_fu](z=0)')
+            '[T_fd](z=0) = [c_in]*[T_{f,in}] + [c_fu]*[T_fu](z=0)')
 
-    def _continuity_condition_head(
+    def _flow_direction(
             self, m_flow_borehole, cp_f, nSegments, segment_ratios=None):
         """ Returns coefficients for the relation
-            [T_fd](z=0) = [c_in]*[T_{f,in}] + [c_fu]*[T_fu](z=0)
+            [T_f](z) = [c_fd]*[T_fd](z) + [c_fu]*[T_fu](z)
         """
         raise NotImplementedError(
             '_continuity_condition_head class method not implemented, '
@@ -957,6 +963,16 @@ def _general_solution(
             'this method should return matrices for the relation: '
             '[T_f](z) = [a_f0]*[T_f](0) + [a_b]*[T_b]')
 
+    def _pipe_connectivity(
+            self, m_flow_borehole, cp_f, nSegments, segment_ratios=None):
+        """ Returns coefficients for the relation
+            [T_{f,out}] = [c_fu]*[T_fu](z=0)
+        """
+        raise NotImplementedError(
+            '_continuity_condition_head class method not implemented, '
+            'this method should return matrices for the relation: '
+            '[T_{f,out}] = [c_fu]*[T_fu](z=0)')
+
     def _update_model_variables(
             self, m_flow_borehole, cp_f, nSegments, segment_ratios):
         """
@@ -1010,6 +1026,12 @@ class SingleUTube(_BasePipe):
     J : int, optional
         Number of multipoles per pipe to evaluate the thermal resistances.
         Default is 2.
+    reversible_flow : bool
+        True to treat a negative mass flow rate as the reversal of flow
+        direction within the borehole. If False, the direction of flow is not
+        reversed when the mass flow rate is negative, and the absolute value is
+        used for calculations.
+        Default is True.
     nPipes : int
         Number of U-Tubes, equals to 1.
     nInlets : int
@@ -1044,7 +1066,9 @@ class SingleUTube(_BasePipe):
        Environment, 25 (8), 1007-1022.
 
     """
-    def __init__(self, pos, r_in, r_out, borehole, k_s, k_g, R_fp, J=2):
+    def __init__(
+            self, pos, r_in, r_out, borehole, k_s, k_g, R_fp, J=2,
+            reversible_flow=True):
         self.pos = pos
         self.r_in = r_in
         self.r_out = r_out
@@ -1053,6 +1077,7 @@ def __init__(self, pos, r_in, r_out, borehole, k_s, k_g, R_fp, J=2):
         self.k_g = k_g
         self.R_fp = R_fp
         self.J = J
+        self.reversible_flow = reversible_flow
         self.nPipes = 1
         self.nInlets = 1
         self.nOutlets = 1
@@ -1206,14 +1231,14 @@ def coefficients_temperature(
         # Downward pipe
         a_in_00 = (1 - delta) * A + beta12 / gamma * B
         a_in_01 = (1 + delta) * B - beta12 / gamma * A
-        a_in[:, 0, 0] = 0.5 * C * (
+        a_in[:, self._iInlet, 0] = 0.5 * C * (
             a_in_00 * np.exp((beta + gamma) * z - 2 * gamma * H)
             + a_in_01 * np.exp((beta - gamma) * z)
             )
         # Upward pipe
         a_in_10 = (1 + delta) * B - beta12 / gamma * A
         a_in_11 = (1 - delta) * A + beta12 / gamma * B
-        a_in[:, 1, 0] = 0.5 * C * (
+        a_in[:, self._iOutlet, 0] = 0.5 * C * (
             a_in_10 * np.exp((beta + gamma) * z - 2 * gamma * H)
             + a_in_11 * np.exp((beta - gamma) * z)
             )
@@ -1242,7 +1267,7 @@ def coefficients_temperature(
         a_b_04 = beta12 / gamma * (beta1 + beta2) / gamma
         a_b_05 = -beta12 / gamma * (beta1 + beta2) / gamma
         # Coefficients for z < z_u
-        a_b[:, 0, :][index] = 0.5 * (
+        a_b[:, self._iInlet, :][index] = 0.5 * (
             a_b_01 * np.exp((beta + gamma) * dz_p - 2 * gamma * H)
             + a_b_02 * np.exp((beta - gamma) * dz_p)
             + a_b_03 * np.exp((beta - gamma) * dz_p - 2 * gamma * H)
@@ -1252,7 +1277,7 @@ def coefficients_temperature(
             + a_b_05 * np.exp((beta - gamma) * dz_p - 2 * gamma * H)
             )
         # Coefficients for z > z_u
-        a_b[:, 0, :][~index] = 0.5 * (
+        a_b[:, self._iInlet, :][~index] = 0.5 * (
             a_b_00 - a_b_04 * np.exp((beta + gamma) * dz_m)
             + a_b_01 * np.exp(-2 * gamma * H)
             + a_b_02
@@ -1270,7 +1295,7 @@ def coefficients_temperature(
         a_b_14 = (delta + 1) * (beta1 + beta2) / gamma
         a_b_15 = (delta - 1) * (beta1 + beta2) / gamma
         # Coefficients for z < z_u
-        a_b[:, 1, :][index] = 0.5 * (
+        a_b[:, self._iOutlet, :][index] = 0.5 * (
             a_b_11 * np.exp((beta + gamma) * dz_p - 2 * gamma * H)
             + a_b_12 * np.exp((beta - gamma) * dz_p)
             + a_b_13 * np.exp((beta - gamma) * dz_p - 2 * gamma * H)
@@ -1280,7 +1305,7 @@ def coefficients_temperature(
             - a_b_15 * np.exp((beta - gamma) * dz_p - 2 * gamma * H)
             )
         # Coefficients for z > z_u
-        a_b[:, 1, :][~index] = 0.5 * (
+        a_b[:, self._iOutlet, :][~index] = 0.5 * (
             a_b_10 - a_b_14 * np.exp((beta + gamma) * dz_m)
             + a_b_11 * np.exp(-2 * gamma * H)
             + a_b_12
@@ -1318,17 +1343,19 @@ def update_thermal_resistances(self, R_fp):
         self._initialize_stored_coefficients()
         return
 
-    def _pipe_connectivity(
+    def _continuity_condition_head(
             self, m_flow_borehole, cp_f, nSegments, segment_ratios=None):
         """
-        Equation that satisfies equal fluid temperatures in both legs of
-        each U-tube pipe at depth (z = H).
+        Build coefficient matrices to evaluate fluid temperatures at depth
+        (z = 0). These coefficients take into account connections between
+        U-tube pipes.
 
         Returns coefficients for the relation:
 
             .. math::
 
-                T_{f,out} = \\mathbf{c_{fu}} \\mathbf{T_{fu}}(z=0)
+                \\mathbf{T_fd}(z=0) = \\mathbf{c_{in}} \\mathbf{T_{f,in}}
+                + \\mathbf{c_{fu}} \\mathbf{T_{fu}}(z=0)
 
         Parameters
         ----------
@@ -1346,7 +1373,9 @@ def _pipe_connectivity(
 
         Returns
         -------
-        c_fu : (nOutlets, nPipes,) array
+        c_in : (nPipes, nInlets,) array
+            Array of coefficients for inlet fluid temperature.
+        c_fu : (nPipes, nPipes,) array
             Array of coefficients for upward fluid temperatures.
 
         """
@@ -1354,24 +1383,26 @@ def _pipe_connectivity(
         self._check_model_variables(
             m_flow_borehole, cp_f, nSegments, segment_ratios)
 
-        # The upward fluid temperature at z=0 is the outlet fluid temperature
-        c_fu = np.array([[1.]])
+        # The inlet is connected to pipe 0
+        c_in = np.array([[1.0]])
+        # The upward pipe is not connected to another pipe (only to the outlet)
+        c_fu = np.zeros((self.nPipes, self.nPipes))
 
-        return c_fu
+        return c_in, c_fu
 
-    def _continuity_condition_head(
+    def _flow_direction(
             self, m_flow_borehole, cp_f, nSegments, segment_ratios=None):
         """
         Build coefficient matrices to evaluate fluid temperatures at depth
-        (z = 0). These coefficients take into account connections between
-        U-tube pipes.
+        (z). These coefficients take into account connections between U-tube
+        pipes.
 
         Returns coefficients for the relation:
 
             .. math::
 
-                \\mathbf{T_fd}(z=0) = \\mathbf{c_{in}} \\mathbf{T_{f,in}}
-                + \\mathbf{c_{fu}} \\mathbf{T_{fu}}(z=0)
+                \\mathbf{T_f}(z) = \\mathbf{c_{fd}} \\mathbf{T_{fd}}(z)
+                + \\mathbf{c_{fu}} \\mathbf{T_{fu}}(z)
 
         Parameters
         ----------
@@ -1389,9 +1420,9 @@ def _continuity_condition_head(
 
         Returns
         -------
-        c_in : (nPipes, nInlets,) array
+        c_fd : (2*nPipes, nPipes,) array
             Array of coefficients for inlet fluid temperature.
-        c_fu : (nPipes, nPipes,) array
+        c_fu : (2*nPipes, nPipes,) array
             Array of coefficients for upward fluid temperatures.
 
         """
@@ -1399,12 +1430,26 @@ def _continuity_condition_head(
         self._check_model_variables(
             m_flow_borehole, cp_f, nSegments, segment_ratios)
 
-        # The inlet is connected to pipe 0
-        c_in = np.array([[1.0]])
-        # The upward pipe is not connected to another pipe (only to the outlet)
-        c_fu = np.zeros((self.nPipes, self.nPipes))
+        if self._is_reversed:
+            # The last nPipes pipes are downward flowing
+            c_fd = np.vstack(
+                (np.zeros((self.nPipes, self.nPipes)),
+                 np.eye(self.nPipes)))
+            # The first nPipes pipes are upward flowing
+            c_fu = np.vstack(
+                (np.eye(self.nPipes),
+                 np.zeros((self.nPipes, self.nPipes))))
+        else:
+            # The first nPipes pipes are downward flowing
+            c_fd = np.vstack(
+                (np.eye(self.nPipes),
+                 np.zeros((self.nPipes, self.nPipes))))
+            # The last nPipes pipes are upward flowing
+            c_fu = np.vstack(
+                (np.zeros((self.nPipes, self.nPipes)),
+                 np.eye(self.nPipes)))
 
-        return c_in, c_fu
+        return c_fd, c_fu
 
     def _general_solution(
             self, z, m_flow_borehole, cp_f, nSegments, segment_ratios=None):
@@ -1449,8 +1494,16 @@ def _general_solution(
         self._check_model_variables(
             m_flow_borehole, cp_f, nSegments, segment_ratios)
 
-        a_f0 = np.block([[[self._f1(z)], [self._f2(z)]],
-                         [[-self._f2(z)], [self._f3(z)]]]).transpose(2, 0, 1)
+        if self._is_reversed:
+            a_f0 = np.block(
+                [[[self._f3(z)], [-self._f2(z)]],
+                 [[self._f2(z)], [self._f1(z)]]]).transpose(
+                     2, 0, 1)
+        else:
+            a_f0 = np.block(
+                [[[self._f1(z)], [self._f2(z)]],
+                 [[-self._f2(z)], [self._f3(z)]]]).transpose(
+                     2, 0, 1)
 
         a_b = np.zeros((len(z_array), 2*self.nPipes, nSegments))
         z_edges = self.b._segment_edges(
@@ -1458,8 +1511,8 @@ def _general_solution(
         dz = np.maximum(np.subtract.outer(z_array, z_edges), 0.)
         dF4 = self._F4(dz)
         dF5 = self._F5(dz)
-        a_b[:, 0, :] = (dF4[:, :-1] - dF4[:, 1:])
-        a_b[:, 1, :] = -(dF5[:, :-1] - dF5[:, 1:])
+        a_b[:, self._iInlet, :] = (dF4[:, :-1] - dF4[:, 1:])
+        a_b[:, self._iOutlet, :] = -(dF5[:, :-1] - dF5[:, 1:])
         # Remove first dimension if z is a scalar
         if np.isscalar(z):
             a_f0 = a_f0[0, :, :]
@@ -1467,6 +1520,47 @@ def _general_solution(
 
         return a_f0, a_b
 
+    def _pipe_connectivity(
+            self, m_flow_borehole, cp_f, nSegments, segment_ratios=None):
+        """
+        Equation that satisfies equal fluid temperatures in both legs of
+        each U-tube pipe at depth (z = H).
+
+        Returns coefficients for the relation:
+
+            .. math::
+
+                T_{f,out} = \\mathbf{c_{fu}} \\mathbf{T_{fu}}(z=0)
+
+        Parameters
+        ----------
+        m_flow_borehole : float or (nInlets,) array
+            Inlet mass flow rate (in kg/s) into the borehole.
+        cp_f : float or (nInlets,) array
+            Fluid specific isobaric heat capacity (in J/kg.degC).
+        nSegments : int
+            Number of borehole segments.
+        segment_ratios : (nSegments,) array, optional
+            Ratio of the borehole length represented by each segment. The sum
+            of ratios must be equal to 1. If segment_ratios==None, segments of
+            equal lengths are considered.
+            Default is None.
+
+        Returns
+        -------
+        c_fu : (nOutlets, nPipes,) array
+            Array of coefficients for upward fluid temperatures.
+
+        """
+        # Check if model variables need to be updated
+        self._check_model_variables(
+            m_flow_borehole, cp_f, nSegments, segment_ratios)
+
+        # The upward fluid temperature at z=0 is the outlet fluid temperature
+        c_fu = np.array([[1.]])
+
+        return c_fu
+
     def _update_model_variables(
             self, m_flow_borehole, cp_f, nSegments, segment_ratios):
         """
@@ -1488,13 +1582,18 @@ def _update_model_variables(
 
         # Format mass flow rate and heat capacity inputs
         self._format_inputs(m_flow_borehole, cp_f, nSegments, segment_ratios)
-        m_flow_in = self._m_flow_in
+        m_flow_in = np.abs(self._m_flow_in)
         cp_in = self._cp_in
+        iInlet = self._iInlet
+        iOutlet = self._iOutlet
 
         # Dimensionless delta-circuit conductances
-        self._beta1 = 1./(self._Rd[0][0]*m_flow_in[0]*cp_in[0])
-        self._beta2 = 1./(self._Rd[1][1]*m_flow_in[0]*cp_in[0])
-        self._beta12 = 1./(self._Rd[0][1]*m_flow_in[0]*cp_in[0])
+        self._beta1 = 1. / (
+            self._Rd[iInlet, iInlet] * m_flow_in[0] * cp_in[0])
+        self._beta2 = 1. / (
+            self._Rd[iOutlet, iOutlet] * m_flow_in[0] * cp_in[0])
+        self._beta12 = 1. / (
+            self._Rd[iInlet, iOutlet] * m_flow_in[0] * cp_in[0])
         self._beta = 0.5*(self._beta2 - self._beta1)
         # Eigenvalues
         self._gamma = np.sqrt(0.25*(self._beta1+self._beta2)**2
@@ -1523,14 +1622,22 @@ def _format_inputs(self, m_flow_borehole, cp_f, nSegments, segment_ratios):
         # Format mass flow rate inputs
         if np.isscalar(m_flow_borehole):
             # Mass flow rate in each fluid circuit
-            m_flow_in = m_flow_borehole*np.ones(self.nInlets)
+            m_flow_in = m_flow_borehole * np.ones(self.nInlets)
         else:
             # Mass flow rate in each fluid circuit
             m_flow_in = m_flow_borehole
         self._m_flow_in = m_flow_in
         # Mass flow rate in pipes
-        m_flow_pipe = np.tile(m_flow_in, 2*self.nPipes)
+        m_flow_pipe = np.array([m_flow_borehole, -m_flow_borehole]).flatten()
         self._m_flow_pipe = m_flow_pipe
+        # Flow direction
+        self._is_reversed = m_flow_borehole < 0. and self.reversible_flow
+        if self._is_reversed:
+            self._iInlet = 1
+            self._iOutlet = 0
+        else:
+            self._iInlet = 0
+            self._iOutlet = 1
 
         # Format heat capacity inputs
         if np.isscalar(cp_f):
@@ -1683,6 +1790,12 @@ class MultipleUTube(_BasePipe):
     J : int, optional
         Number of multipoles per pipe to evaluate the thermal resistances.
         Default is 2.
+    reversible_flow : bool
+        True to treat a negative mass flow rate as the reversal of flow
+        direction within the borehole. If False, the direction of flow is not
+        reversed when the mass flow rate is negative, and the absolute value is
+        used for calculations.
+        Default is True.
     nPipes : int
         Number of U-Tubes.
     config : str, defaults to 'parallel'
@@ -1712,20 +1825,20 @@ class MultipleUTube(_BasePipe):
     .. [#Cimmino2016] Cimmino, M. (2016). Fluid and borehole wall temperature
        profiles in vertical geothermal boreholes with multiple U-tubes.
        Renewable Energy, 96, 137-147.
+    .. [#Multiple-Cimmin2019] Cimmino, M. (2019). Semi-analytical method for
+       g-function calculation of bore fields with series- and
+       parallel-connected boreholes. Science and Technology for the Built
+       Environment, 25 (8), 1007-1022.
     .. [#Cimmino2024] Cimmino, M. (2024). g-Functions for fields of series- and
        parallel-connected boreholes with variable fluid mass flow rate and
        reversible flow direction. Renewable Energy, 228, 120661.
     .. [#Multiple-Claesson2011b] Claesson, J., & Hellstrom, G. (2011).
        Multipole method to calculate borehole thermal resistances in a borehole
        heat exchanger. HVAC&R Research, 17(6), 895-911.
-    .. [#Multiple-Cimmin2019] Cimmino, M. (2019). Semi-analytical method for
-       g-function calculation of bore fields with series- and
-       parallel-connected boreholes. Science and Technology for the Built
-       Environment, 25 (8), 1007-1022.
 
     """
-    def __init__(self, pos, r_in, r_out, borehole, k_s,
-                 k_g, R_fp, nPipes, config='parallel', J=2):
+    def __init__(self, pos, r_in, r_out, borehole, k_s, k_g, R_fp, nPipes,
+                 config='parallel', J=2, reversible_flow=True):
         self.pos = pos
         self.r_in = r_in
         self.r_out = r_out
@@ -1734,6 +1847,7 @@ def __init__(self, pos, r_in, r_out, borehole, k_s,
         self.k_g = k_g
         self.R_fp = R_fp
         self.J = J
+        self.reversible_flow = reversible_flow
         self.nPipes = nPipes
         self.nInlets = 1
         self.nOutlets = 1
@@ -1767,17 +1881,19 @@ def update_thermal_resistances(self, R_fp):
         self._initialize_stored_coefficients()
         return
 
-    def _pipe_connectivity(
+    def _continuity_condition_head(
             self, m_flow_borehole, cp_f, nSegments, segment_ratios=None):
         """
-        Equation that satisfies equal fluid temperatures in both legs of
-        each U-tube pipe at depth (z = H).
+        Build coefficient matrices to evaluate fluid temperatures at depth
+        (z = 0). These coefficients take into account connections between
+        U-tube pipes.
 
         Returns coefficients for the relation:
 
             .. math::
 
-                T_{f,out} = \\mathbf{c_{fu}} \\mathbf{T_{fu}}(z=0)
+                \\mathbf{T_fd}(z=0) = \\mathbf{c_{in}} \\mathbf{T_{f,in}}
+                + \\mathbf{c_{fu}} \\mathbf{T_{fu}}(z=0)
 
         Parameters
         ----------
@@ -1795,42 +1911,53 @@ def _pipe_connectivity(
 
         Returns
         -------
-        c_fu : (nOutlets, nPipes,) array
+        c_in : (nPipes, nInlets,) array
+            Array of coefficients for inlet fluid temperature.
+        c_fu : (nPipes, nPipes,) array
             Array of coefficients for upward fluid temperatures.
 
         """
         # Check if model variables need to be updated
         self._check_model_variables(
             m_flow_borehole, cp_f, nSegments, segment_ratios)
-        m_flow_pipe = self._m_flow_pipe[-self.nPipes:]
-        cp_pipe = self._cp_pipe[-self.nPipes:]
 
         if self.config == 'parallel':
-            # The outlet temperature is a result of mixing from all the upward
-            # flowing pipes
-            c_fu = m_flow_pipe * cp_pipe / np.sum(m_flow_pipe * cp_pipe)
+            # The inlet is connected to all downward flowing pipes
+            c_in = np.ones(self.nPipes)
+            # None of the upward flowing pipes are connected to another pipe
+            c_fu = np.zeros((self.nPipes, self.nPipes))
         elif self.config == 'series':
-            # Only the last pipe is connected to the outlet
-            c_fu = np.concatenate((np.zeros(self.nPipes-1), np.ones(1)))
+            if self._is_reversed:
+                # The inlet is connected to the last downward flowing pipe
+                c_in = np.concatenate((np.zeros(self.nPipes-1), np.ones(1)))
+                # Each upward flowing pipe is connected to a downward flowing pipe
+                # (except for the last one connected to the outlet)
+                c_fu = np.eye(self.nPipes, k=+1)
+            else:
+                # The inlet is connected to the first downward flowing pipe
+                c_in = np.concatenate((np.ones(1), np.zeros(self.nPipes-1)))
+                # Each upward flowing pipe is connected to a downward flowing pipe
+                # (except for the last one connected to the outlet)
+                c_fu = np.eye(self.nPipes, k=-1)
         else:
             raise NotImplementedError(
                 f"Configuration '{self.config}' not implemented.")
 
-        return c_fu[np.newaxis, :]
+        return c_in[:, np.newaxis], c_fu
 
-    def _continuity_condition_head(
+    def _flow_direction(
             self, m_flow_borehole, cp_f, nSegments, segment_ratios=None):
         """
         Build coefficient matrices to evaluate fluid temperatures at depth
-        (z = 0). These coefficients take into account connections between
-        U-tube pipes.
+        (z). These coefficients take into account connections between U-tube
+        pipes.
 
         Returns coefficients for the relation:
 
             .. math::
 
-                \\mathbf{T_fd}(z=0) = \\mathbf{c_{in}} \\mathbf{T_{f,in}}
-                + \\mathbf{c_{fu}} \\mathbf{T_{fu}}(z=0)
+                \\mathbf{T_f}(z) = \\mathbf{c_{fd}} \\mathbf{T_{fd}}(z)
+                + \\mathbf{c_{fu}} \\mathbf{T_{fu}}(z)
 
         Parameters
         ----------
@@ -1848,9 +1975,9 @@ def _continuity_condition_head(
 
         Returns
         -------
-        c_in : (nPipes, nInlets,) array
+        c_fd : (2*nPipes, nPipes,) array
             Array of coefficients for inlet fluid temperature.
-        c_fu : (nPipes, nPipes,) array
+        c_fu : (2*nPipes, nPipes,) array
             Array of coefficients for upward fluid temperatures.
 
         """
@@ -1858,22 +1985,26 @@ def _continuity_condition_head(
         self._check_model_variables(
             m_flow_borehole, cp_f, nSegments, segment_ratios)
 
-        if self.config == 'parallel':
-            # The inlet is connected to all downward flowing pipes
-            c_in = np.ones(self.nPipes)
-            # None of the upward flowing pipes are connected to another pipe
-            c_fu = np.zeros((self.nPipes, self.nPipes))
-        elif self.config == 'series':
-            # The inlet is connected to the first downward flowing pipe
-            c_in = np.concatenate((np.ones(1), np.zeros(self.nPipes-1)))
-            # Each upward flowing pipe is connected to a downward flowing pipe
-            # (except for the last one connected to the outlet)
-            c_fu = np.eye(self.nPipes, k=-1)
+        if self._is_reversed:
+            # The last nPipes pipes are downward flowing
+            c_fd = np.vstack(
+                (np.zeros((self.nPipes, self.nPipes)),
+                 np.eye(self.nPipes)))
+            # The first nPipes pipes are upward flowing
+            c_fu = np.vstack(
+                (np.eye(self.nPipes),
+                 np.zeros((self.nPipes, self.nPipes))))
         else:
-            raise NotImplementedError(
-                f"Configuration '{self.config}' not implemented.")
+            # The first nPipes pipes are downward flowing
+            c_fd = np.vstack(
+                (np.eye(self.nPipes),
+                 np.zeros((self.nPipes, self.nPipes))))
+            # The last nPipes pipes are upward flowing
+            c_fu = np.vstack(
+                (np.zeros((self.nPipes, self.nPipes)),
+                 np.eye(self.nPipes)))
 
-        return c_in[:, np.newaxis], c_fu
+        return c_fd, c_fu
 
     def _general_solution(
             self, z, m_flow_borehole, cp_f, nSegments, segment_ratios=None):
@@ -1944,6 +2075,61 @@ def _general_solution(
 
         return a_f0, a_b
 
+    def _pipe_connectivity(
+            self, m_flow_borehole, cp_f, nSegments, segment_ratios=None):
+        """
+        Equation that satisfies equal fluid temperatures in both legs of
+        each U-tube pipe at depth (z = H).
+
+        Returns coefficients for the relation:
+
+            .. math::
+
+                T_{f,out} = \\mathbf{c_{fu}} \\mathbf{T_{fu}}(z=0)
+
+        Parameters
+        ----------
+        m_flow_borehole : float or (nInlets,) array
+            Inlet mass flow rate (in kg/s) into the borehole.
+        cp_f : float or (nInlets,) array
+            Fluid specific isobaric heat capacity (in J/kg.degC).
+        nSegments : int
+            Number of borehole segments.
+        segment_ratios : (nSegments,) array, optional
+            Ratio of the borehole length represented by each segment. The sum
+            of ratios must be equal to 1. If segment_ratios==None, segments of
+            equal lengths are considered.
+            Default is None.
+
+        Returns
+        -------
+        c_fu : (nOutlets, nPipes,) array
+            Array of coefficients for upward fluid temperatures.
+
+        """
+        # Check if model variables need to be updated
+        self._check_model_variables(
+            m_flow_borehole, cp_f, nSegments, segment_ratios)
+        m_flow_pipe = np.abs(self._m_flow_pipe[:self.nPipes])
+        cp_pipe = self._cp_pipe[:self.nPipes]
+
+        if self.config == 'parallel':
+            # The outlet temperature is a result of mixing from all the upward
+            # flowing pipes
+            c_fu = m_flow_pipe * cp_pipe / np.sum(m_flow_pipe * cp_pipe)
+        elif self.config == 'series':
+            if self._is_reversed:
+                # Only the first pipe is connected to the outlet
+                c_fu = np.concatenate((np.ones(1), np.zeros(self.nPipes-1)))
+            else:
+                # Only the last pipe is connected to the outlet
+                c_fu = np.concatenate((np.zeros(self.nPipes-1), np.ones(1)))
+        else:
+            raise NotImplementedError(
+                f"Configuration '{self.config}' not implemented.")
+
+        return c_fu[np.newaxis, :]
+
     def _update_model_variables(
             self, m_flow_borehole, cp_f, nSegments, segment_ratios):
         """
@@ -1975,8 +2161,6 @@ def _update_model_variables(
         for i in range(2*nPipes):
             self._A[i, i] = -self._A[i, i] - sum(
                 [self._A[i, j] for j in range(2*nPipes) if not i == j])
-        for i in range(nPipes, 2*nPipes):
-            self._A[i, :] = - self._A[i, :]
         self._sumA = np.sum(self._A, axis=1)
         # Eigenvalues and eigenvectors of A
         self._L, self._V = np.linalg.eig(self._A)
@@ -2008,13 +2192,20 @@ def _format_inputs(self, m_flow_borehole, cp_f, nSegments, segment_ratios):
         # Format mass flow rate inputs
         # Mass flow rate in pipes
         if self.config.lower() == 'parallel':
-            m_flow_pipe = np.tile(m_flow_borehole/nPipes, 2*self.nPipes)
+            m_flow_pipe = np.tile(
+                np.abs(m_flow_borehole) / nPipes,
+                self.nPipes)
         elif self.config.lower() == 'series':
-            m_flow_pipe = np.tile(m_flow_borehole, 2*self.nPipes)
-        self._m_flow_pipe = m_flow_pipe
+            m_flow_pipe = np.tile(
+                np.abs(m_flow_borehole),
+                self.nPipes)
+        self._m_flow_pipe = np.sign(m_flow_borehole) * np.concatenate(
+            (m_flow_pipe, -m_flow_pipe))
         # Mass flow rate in each fluid circuit
         m_flow_in = np.atleast_1d(m_flow_borehole)
         self._m_flow_in = m_flow_in
+        # Flow direction
+        self._is_reversed = m_flow_borehole < 0. and self.reversible_flow
 
         # Format heat capacity inputs
         # Heat capacity in each fluid circuit
@@ -2055,6 +2246,12 @@ class IndependentMultipleUTube(MultipleUTube):
     J : int, optional
         Number of multipoles per pipe to evaluate the thermal resistances.
         Default is 2.
+    reversible_flow : bool
+        True to treat a negative mass flow rate as the reversal of flow
+        direction within the borehole. If False, the direction of flow is not
+        reversed when the mass flow rate is negative, and the absolute value is
+        used for calculations.
+        Default is True.
     nPipes : int
         Number of U-Tubes.
     nInlets : int
@@ -2077,16 +2274,16 @@ class IndependentMultipleUTube(MultipleUTube):
     .. [#Cimmino2016b] Cimmino, M. (2016). Fluid and borehole wall temperature
        profiles in vertical geothermal boreholes with multiple U-tubes.
        Renewable Energy, 96, 137-147.
-    .. [#Cimmino2024b] Cimmino, M. (2024). g-Functions for fields of series- and
-       parallel-connected boreholes with variable fluid mass flow rate and
+    .. [#Cimmino2024b] Cimmino, M. (2024). g-Functions for fields of series-
+       and parallel-connected boreholes with variable fluid mass flow rate and
        reversible flow direction. Renewable Energy, 228, 120661.
     .. [#Independent-Claesson2011b] Claesson, J., & Hellstrom, G. (2011).
        Multipole method to calculate borehole thermal resistances in a borehole
        heat exchanger. HVAC&R Research, 17(6), 895-911.
 
     """
-    def __init__(self, pos, r_in, r_out, borehole, k_s,
-                 k_g, R_fp, nPipes, J=2):
+    def __init__(self, pos, r_in, r_out, borehole, k_s, k_g, R_fp, nPipes,
+                 J=2, reversible_flow=True):
         self.pos = pos
         self.r_in = r_in
         self.r_out = r_out
@@ -2095,6 +2292,7 @@ def __init__(self, pos, r_in, r_out, borehole, k_s,
         self.k_g = k_g
         self.R_fp = R_fp
         self.J = J
+        self.reversible_flow = reversible_flow
         self.nPipes = nPipes
         self.nInlets = nPipes
         self.nOutlets = nPipes
@@ -2127,6 +2325,58 @@ def update_thermal_resistances(self, R_fp):
         self._initialize_stored_coefficients()
         return
 
+    def _flow_direction(
+            self, m_flow_borehole, cp_f, nSegments, segment_ratios=None):
+        """
+        Build coefficient matrices to evaluate fluid temperatures at depth
+        (z). These coefficients take into account connections between U-tube
+        pipes.
+
+        Returns coefficients for the relation:
+
+            .. math::
+
+                \\mathbf{T_f}(z) = \\mathbf{c_{fd}} \\mathbf{T_{fd}}(z)
+                + \\mathbf{c_{fu}} \\mathbf{T_{fu}}(z)
+
+        Parameters
+        ----------
+        m_flow_borehole : float or (nInlets,) array
+            Inlet mass flow rate (in kg/s) into the borehole.
+        cp_f : float or (nInlets,) array
+            Fluid specific isobaric heat capacity (in J/kg.degC).
+        nSegments : int
+            Number of borehole segments.
+        segment_ratios : (nSegments,) array, optional
+            Ratio of the borehole length represented by each segment. The sum
+            of ratios must be equal to 1. If segment_ratios==None, segments of
+            equal lengths are considered.
+            Default is None.
+
+        Returns
+        -------
+        c_fd : (2*nPipes, nPipes,) array
+            Array of coefficients for inlet fluid temperature.
+        c_fu : (2*nPipes, nPipes,) array
+            Array of coefficients for upward fluid temperatures.
+
+        """
+        # Check if model variables need to be updated
+        self._check_model_variables(
+            m_flow_borehole, cp_f, nSegments, segment_ratios)
+
+        c_fd = np.zeros((2*self.nPipes, self.nPipes))
+        c_fu = np.zeros((2*self.nPipes, self.nPipes))
+        for i, is_reversed in enumerate(self._is_reversed):
+            if is_reversed:
+                c_fd[i+self.nPipes, i] = 1.
+                c_fu[i, i] = 1.
+            else:
+                c_fd[i, i] = 1.
+                c_fu[i+self.nPipes, i] = 1.
+
+        return c_fd, c_fu
+
     def _pipe_connectivity(
             self, m_flow_borehole, cp_f, nSegments, segment_ratios=None):
         """
@@ -2227,8 +2477,11 @@ def _format_inputs(self, m_flow_borehole, cp_f, nSegments, segment_ratios):
                 'Incorrect length of mass flow vector.')
         self._m_flow_in = m_flow_in
         # Mass flow rate in pipes
-        m_flow_pipe = np.tile(m_flow_in, 2)
-        self._m_flow_pipe = m_flow_pipe
+        self._m_flow_pipe = np.concatenate((m_flow_in, -m_flow_in))
+        # Flow direction
+        self._is_reversed = np.logical_and(
+            m_flow_borehole < 0.,
+            self.reversible_flow)
 
         # Format heat capacity inputs
         # Heat capacity in each fluid circuit
@@ -2280,6 +2533,12 @@ class Coaxial(SingleUTube):
     J : int, optional
         Number of multipoles per pipe to evaluate the thermal resistances.
         Default is 2.
+    reversible_flow : bool
+        True to treat a negative mass flow rate as the reversal of flow
+        direction within the borehole. If False, the direction of flow is not
+        reversed when the mass flow rate is negative, and the absolute value is
+        used for calculations.
+        Default is True.
     nPipes : int
         Number of U-Tubes, equals to 1.
     nInlets : int
@@ -2314,7 +2573,8 @@ class Coaxial(SingleUTube):
        Environment, 25 (8), 1007-1022.
 
     """
-    def __init__(self, pos, r_in, r_out, borehole, k_s, k_g, R_ff, R_fp, J=2):
+    def __init__(self, pos, r_in, r_out, borehole, k_s, k_g, R_ff, R_fp, J=2,
+                 reversible_flow=True):
         if isinstance(pos, tuple):
             pos = [pos]
         self.pos = pos
@@ -2326,6 +2586,7 @@ def __init__(self, pos, r_in, r_out, borehole, k_s, k_g, R_ff, R_fp, J=2):
         self.R_ff = R_ff
         self.R_fp = R_fp
         self.J = J
+        self.reversible_flow = reversible_flow
         self.nPipes = 1
         self.nInlets = 1
         self.nOutlets = 1
@@ -2804,7 +3065,7 @@ def fluid_friction_factor_circular_pipe(
     # Relative roughness
     E = epsilon / D
     # Fluid velocity
-    V_flow = m_flow_pipe / rho_f
+    V_flow = np.abs(m_flow_pipe) / rho_f
     A_cs = pi * r_in**2
     V = V_flow / A_cs
     # Reynolds number
@@ -2881,7 +3142,7 @@ def convective_heat_transfer_coefficient_circular_pipe(
     # Hydraulic diameter
     D = 2.*r_in
     # Fluid velocity
-    V_flow = m_flow_pipe / rho_f
+    V_flow = np.abs(m_flow_pipe) / rho_f
     A_cs = pi * r_in**2
     V = V_flow / A_cs
     # Reynolds number
@@ -2993,7 +3254,7 @@ def convective_heat_transfer_coefficient_concentric_annulus(
     # Cross-sectional area of the annulus region
     A_c = pi * ((r_a_out ** 2) - (r_a_in ** 2))
     # Volume flow rate
-    V_dot = m_flow_pipe / rho_f
+    V_dot = np.abs(m_flow_pipe) / rho_f
     # Average velocity
     V = V_dot / A_c
     # Reynolds number


Test patch:
diff --git a/tests/conftest.py b/tests/conftest.py
index 3fdd962..b14f4d4 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -154,7 +154,41 @@ def double_Utube_parallel(single_borehole):
     r_in = 0.015            # Pipe inner radius [m]
     epsilon = 1.0e-06       # Pipe surface roughness [m]
     m_flow_borehole = 0.05  # Nominal fluid mass flow rate [kg/s]
-    m_flow_pipe = m_flow_borehole
+    m_flow_pipe = m_flow_borehole / 2
+    # Fluid is propylene-glycol (20 %) at 20 degC
+    fluid = gt.media.Fluid('MPG', 20.)
+    # Pipe thermal resistance [m.K/W]
+    R_p = gt.pipes.conduction_thermal_resistance_circular_pipe(
+        r_in, r_out, k_p)
+    # Convection heat transfer coefficient [W/m2.K]
+    h_f = gt.pipes.convective_heat_transfer_coefficient_circular_pipe(
+        m_flow_pipe, r_in, fluid.mu, fluid.rho, fluid.k, fluid.cp,
+        epsilon)
+    # Film thermal resistance [m.K/W]
+    R_f = 1.0 / (h_f * 2 * np.pi * r_in)
+    # Initialize pipe
+    doubleUTube = gt.pipes.MultipleUTube(
+            pos_pipes, r_in, r_out, borehole, k_s, k_g, R_f + R_p, 2,
+            config='parallel')
+    return doubleUTube
+
+
+@pytest.fixture
+def double_Utube_parallel_asymmetrical(single_borehole):
+    # Extract borehole from fixture
+    borehole = single_borehole[0]
+    # Pipe positions [m]
+    D_s = 0.05
+    pos_pipes = [
+        (-D_s, 0.), (0.35*D_s, -0.6*D_s), (0.9*D_s, 0.), (0., 0.8*D_s)]
+    k_s = 2.0               # Ground thermal conductivity [W/m.K]
+    k_g = 1.0               # Grout thermal conductivity [W/m.K]
+    k_p = 0.4               # Pipe thermal conductivity [W/m.K]
+    r_out = 0.02            # Pipe outer radius [m]
+    r_in = 0.015            # Pipe inner radius [m]
+    epsilon = 1.0e-06       # Pipe surface roughness [m]
+    m_flow_borehole = 0.05  # Nominal fluid mass flow rate [kg/s]
+    m_flow_pipe = m_flow_borehole / 2
     # Fluid is propylene-glycol (20 %) at 20 degC
     fluid = gt.media.Fluid('MPG', 20.)
     # Pipe thermal resistance [m.K/W]
@@ -206,6 +240,40 @@ def double_Utube_series(single_borehole):
     return doubleUTube
 
 
+@pytest.fixture
+def double_Utube_series_asymmetrical(single_borehole):
+    # Extract borehole from fixture
+    borehole = single_borehole[0]
+    # Pipe positions [m]
+    D_s = 0.05
+    pos_pipes = [
+        (-D_s, 0.), (0.35*D_s, -0.6*D_s), (0.9*D_s, 0.), (0., 0.8*D_s)]
+    k_s = 2.0               # Ground thermal conductivity [W/m.K]
+    k_g = 1.0               # Grout thermal conductivity [W/m.K]
+    k_p = 0.4               # Pipe thermal conductivity [W/m.K]
+    r_out = 0.02            # Pipe outer radius [m]
+    r_in = 0.015            # Pipe inner radius [m]
+    epsilon = 1.0e-06       # Pipe surface roughness [m]
+    m_flow_borehole = 0.05  # Nominal fluid mass flow rate [kg/s]
+    m_flow_pipe = m_flow_borehole
+    # Fluid is propylene-glycol (20 %) at 20 degC
+    fluid = gt.media.Fluid('MPG', 20.)
+    # Pipe thermal resistance [m.K/W]
+    R_p = gt.pipes.conduction_thermal_resistance_circular_pipe(
+        r_in, r_out, k_p)
+    # Convection heat transfer coefficient [W/m2.K]
+    h_f = gt.pipes.convective_heat_transfer_coefficient_circular_pipe(
+        m_flow_pipe, r_in, fluid.mu, fluid.rho, fluid.k, fluid.cp,
+        epsilon)
+    # Film thermal resistance [m.K/W]
+    R_f = 1.0 / (h_f * 2 * np.pi * r_in)
+    # Initialize pipe
+    doubleUTube = gt.pipes.MultipleUTube(
+            pos_pipes, r_in, r_out, borehole, k_s, k_g, R_f + R_p, 2,
+            config='series')
+    return doubleUTube
+
+
 @pytest.fixture
 def double_Utube_independent(single_borehole):
     # Extract borehole from fixture
@@ -585,3 +653,45 @@ def ten_boreholes_network_rectangular_series(ten_boreholes_rectangular):
         boreField, UTubes, bore_connectivity=bore_connectivity,
         m_flow_network=m_flow_network, cp_f=fluid.cp)
     return network
+
+
+@pytest.fixture
+def three_boreholes_network_series_unequal(three_boreholes_unequal):
+    # Extract bore field from fixture
+    boreField = three_boreholes_unequal
+    # Pipe positions [m]
+    D_s = 0.04
+    pos_pipes = [(-D_s, 0.), (D_s, 0.)]
+    k_s = 2.0               # Ground thermal conductivity [W/m.K]
+    k_g = 1.0               # Grout thermal conductivity [W/m.K]
+    k_p = 0.4               # Pipe thermal conductivity [W/m.K]
+    r_out = 0.02            # Pipe outer radius [m]
+    r_in = 0.015            # Pipe inner radius [m]
+    epsilon = 1.0e-06       # Pipe surface roughness [m]
+    m_flow_borehole = 0.05  # Nominal fluid mass flow rate [kg/s]
+    m_flow_pipe = m_flow_borehole
+    m_flow_network = m_flow_borehole
+    # Fluid is propylene-glycol (20 %) at 20 degC
+    fluid = gt.media.Fluid('MPG', 20.)
+    # Pipe thermal resistance [m.K/W]
+    R_p = gt.pipes.conduction_thermal_resistance_circular_pipe(
+        r_in, r_out, k_p)
+    # Convection heat transfer coefficient [W/m2.K]
+    h_f = gt.pipes.convective_heat_transfer_coefficient_circular_pipe(
+        m_flow_pipe, r_in, fluid.mu, fluid.rho, fluid.k, fluid.cp,
+        epsilon)
+    # Film thermal resistance [m.K/W]
+    R_f = 1.0 / (h_f * 2 * np.pi * r_in)
+
+    # Build network of series-connected boreholes
+    bore_connectivity = list(range(-1, len(boreField)-1))
+    UTubes = []
+    for borehole in boreField:
+        UTube = gt.pipes.SingleUTube(
+            pos_pipes, r_in, r_out, borehole, k_s, k_g, R_f + R_p)
+        UTubes.append(UTube)
+    # Initialize network
+    network = gt.networks.Network(
+        boreField, UTubes, bore_connectivity=bore_connectivity,
+        m_flow_network=m_flow_network, cp_f=fluid.cp)
+    return network
diff --git a/tests/networks_test.py b/tests/networks_test.py
new file mode 100644
index 0000000..03a447a
--- /dev/null
+++ b/tests/networks_test.py
@@ -0,0 +1,269 @@
+# -*- coding: utf-8 -*-
+""" Test suite for networks module.
+"""
+import pytest
+
+import numpy as np
+
+import pygfunction as gt
+
+
+# =============================================================================
+# Test network classes
+# =============================================================================
+# Test get_inlet_temperature
+@pytest.mark.parametrize("network_fixture, m_flow_network, segment_ratios, T_b, nSegments, expected", [
+    # Single borehole
+    ('single_borehole_network', 0.2, None, 1., 1, np.array([5.])),
+    ('single_borehole_network', 0.2, None, np.array([1., 2., 3., 1.]), 4, np.array([5.])),
+    ('single_borehole_network', 0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 4, np.array([5.])),
+    ('single_borehole_network', -0.2, None, 1., 1, np.array([5.])),
+    ('single_borehole_network', -0.2, None, np.array([1., 2., 3., 1.]), 4, np.array([5.])),
+    ('single_borehole_network', -0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 4, np.array([5.])),
+    # Ten boreholes (Parallel-connected)
+    ('ten_boreholes_network_rectangular', 2., None, 1., 1, np.full(10, 5.)),
+    ('ten_boreholes_network_rectangular', 2., None, np.tile([1., 2., 3., 1.], 10), 4, np.full(10, 5.)),
+    ('ten_boreholes_network_rectangular', 2., np.array([0.1, 0.35, 0.40, 0.15]), np.tile([1., 2., 3., 1.], 10), 4, np.full(10, 5.)),
+    ('ten_boreholes_network_rectangular', -2., None, 1., 1, np.full(10, 5.)),
+    ('ten_boreholes_network_rectangular', -2., None, np.tile([1., 2., 3., 1.], 10), 4, np.full(10, 5.)),
+    ('ten_boreholes_network_rectangular', -2., np.array([0.1, 0.35, 0.40, 0.15]), np.tile([1., 2., 3., 1.], 10), 4, np.full(10, 5.)),
+    # Three boreholes (Series-connected, unequal lengths)
+    ('three_boreholes_network_series_unequal', 0.2, None, 1., 1, np.array([5., 2.85537985, 2.12376334])),
+    ('three_boreholes_network_series_unequal', 0.2, None, np.tile([1., 2., 3., 1.], 3), 4, np.array([5., 3.25437509, 2.66026346])),
+    ('three_boreholes_network_series_unequal', 0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.tile([1., 2., 3., 1.], 3), 4, np.array([5., 3.46848051, 2.94756734])),
+    ('three_boreholes_network_series_unequal', -0.2, None, 1., 1, np.array([2.053936, 2.74009201, 5.])),
+    ('three_boreholes_network_series_unequal', -0.2, None, np.tile([1., 2., 3., 1.], 3), 4, np.array([2.60290661, 3.15967656, 5.])),
+    ('three_boreholes_network_series_unequal', -0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.tile([1., 2., 3., 1.], 3), 4, np.array([2.89712276, 3.38519442, 5.])),
+    ])
+def test_inlet_temperature(network_fixture, m_flow_network, segment_ratios, T_b, nSegments, expected, request):
+    # Extract pipe from fixture
+    network = request.getfixturevalue(network_fixture)
+    # Fluid is propylene-glycol 20%
+    fluid = gt.media.Fluid('MPG', 20.)
+    T_f_in = 5.0           # Total heat transfer rate [W]
+    # Inlet fluid temperature [degC]
+    T_f_in_borehole = network.get_inlet_temperature(
+        T_f_in, T_b, m_flow_network, fluid.cp, nSegments, segment_ratios=segment_ratios)
+    assert np.allclose(T_f_in_borehole, expected)
+
+
+# Test get_outlet_temperature
+@pytest.mark.parametrize("network_fixture, m_flow_network, segment_ratios, T_b, nSegments, expected", [
+    # Single borehole
+    ('single_borehole_network', 0.2, None, 1., 1, 2.712371852688313),
+    ('single_borehole_network', 0.2, None, np.array([1., 2., 3., 1.]), 4, 3.1377635748663573),
+    ('single_borehole_network', 0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 4, 3.3661222931515784),
+    ('single_borehole_network', -0.2, None, 1., 1, 2.712371852688313),
+    ('single_borehole_network', -0.2, None, np.array([1., 2., 3., 1.]), 4, 3.1377635748663573),
+    ('single_borehole_network', -0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 4, 3.3661222931515784),
+    # Ten boreholes (Parallel-connected)
+    ('ten_boreholes_network_rectangular', 2., None, 1., 1, np.full(10, 2.712371852688313)),
+    ('ten_boreholes_network_rectangular', 2., None, np.tile([1., 2., 3., 1.], 10), 4, np.full(10, 3.1377635748663573)),
+    ('ten_boreholes_network_rectangular', 2., np.array([0.1, 0.35, 0.40, 0.15]), np.tile([1., 2., 3., 1.], 10), 4, np.full(10, 3.3661222931515784)),
+    ('ten_boreholes_network_rectangular', -2., None, 1., 1, np.full(10, 2.712371852688313)),
+    ('ten_boreholes_network_rectangular', -2., None, np.tile([1., 2., 3., 1.], 10), 4, np.full(10, 3.1377635748663573)),
+    ('ten_boreholes_network_rectangular', -2., np.array([0.1, 0.35, 0.40, 0.15]), np.tile([1., 2., 3., 1.], 10), 4, np.full(10, 3.3661222931515784)),
+    # Three boreholes (Series-connected, unequal lengths)
+    ('three_boreholes_network_series_unequal', 0.2, None, 1., 1, np.array([2.85537985, 2.12376334, 1.4888629])),
+    ('three_boreholes_network_series_unequal', 0.2, None, np.tile([1., 2., 3., 1.], 3), 4, np.array([3.25437509, 2.66026346, 2.14183735])),
+    ('three_boreholes_network_series_unequal', 0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.tile([1., 2., 3., 1.], 3), 4, np.array([3.46848051, 2.94756734, 2.492339])),
+    ('three_boreholes_network_series_unequal', -0.2, None, 1., 1, np.array([1.4888629, 2.053936, 2.74009201])),
+    ('three_boreholes_network_series_unequal', -0.2, None, np.tile([1., 2., 3., 1.], 3), 4, np.array([2.1424954, 2.60290661, 3.15967656])),
+    ('three_boreholes_network_series_unequal', -0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.tile([1., 2., 3., 1.], 3), 4, np.array([2.49307149, 2.89712276, 3.38519442])),
+    ])
+def test_outlet_temperature(
+        network_fixture, m_flow_network, segment_ratios, T_b, nSegments, expected, request):
+    # Extract pipe from fixture
+    network = request.getfixturevalue(network_fixture)
+    # Fluid is propylene-glycol 20%
+    fluid = gt.media.Fluid('MPG', 20.)
+    T_f_in = 5.0            # Inlet fluid temperature [degC]
+    # Outlet fluid temperature [degC]
+    T_f_out_borehole = network.get_outlet_temperature(
+        T_f_in, T_b, m_flow_network, fluid.cp, nSegments, segment_ratios=segment_ratios)
+    assert np.allclose(T_f_out_borehole, expected)
+
+
+# Test get_borehole_heat_extraction_rate
+@pytest.mark.parametrize("network_fixture, m_flow_network, segment_ratios, T_b, nSegments, expected", [
+    # Single borehole
+    ('single_borehole_network', 0.2, None, 1., 1, np.array([-1819.4736348927008])),
+    ('single_borehole_network', 0.2, None, np.array([1., 2., 3., 1.]), 4, np.array([-507.98022943, -330.29924271, -155.92399643, -486.93326314])),
+    ('single_borehole_network', 0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 4, np.array([-212.48822505, -496.72428991, -284.6464931, -305.65175979])),
+    ('single_borehole_network', -0.2, None, 1., 1, np.array([-1819.4736348927008])),
+    ('single_borehole_network', -0.2, None, np.array([1., 2., 3., 1.]), 4, np.array([-507.98022943, -330.29924271, -155.92399643, -486.93326314])),
+    ('single_borehole_network', -0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 4, np.array([-212.48822505, -496.72428991, -284.6464931, -305.65175979])),
+    # Ten boreholes (Parallel-connected)
+    ('ten_boreholes_network_rectangular', 2., None, 1., 1, np.full(10, -1819.4736348927008)),
+    ('ten_boreholes_network_rectangular', 2., None, np.tile([1., 2., 3., 1.], 10), 4, np.tile([-507.98022943, -330.29924271, -155.92399643, -486.93326314], 10)),
+    ('ten_boreholes_network_rectangular', 2., np.array([0.1, 0.35, 0.40, 0.15]), np.tile([1., 2., 3., 1.], 10), 4, np.tile([-212.48822505, -496.72428991, -284.6464931, -305.65175979], 10)),
+    ('ten_boreholes_network_rectangular', -2., None, 1., 1, np.full(10, -1819.4736348927008)),
+    ('ten_boreholes_network_rectangular', -2., None, np.tile([1., 2., 3., 1.], 10), 4, np.tile([-507.98022943, -330.29924271, -155.92399643, -486.93326314], 10)),
+    ('ten_boreholes_network_rectangular', -2., np.array([0.1, 0.35, 0.40, 0.15]), np.tile([1., 2., 3., 1.], 10), 4, np.tile([-212.48822505, -496.72428991, -284.6464931, -305.65175979], 10)),
+    # Three boreholes (Series-connected, unequal lengths)
+    ('three_boreholes_network_series_unequal', 0.2, None, 1., 1, np.array([-1705.73168677, -581.89393685, -504.97044758])),
+    ('three_boreholes_network_series_unequal', 0.2, None, np.tile([1., 2., 3., 1.], 3), 4, np.array([-472.32371371, -310.57683056, -151.75246159, -453.73631591, -194.1424052, -93.48480165, 6.82853736, -191.73025237, -231.10379372, -62.2454343, 105.87215406, -224.85507901])),
+    ('three_boreholes_network_series_unequal', 0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.tile([1., 2., 3., 1.], 3), 4, np.array([-196.92259261, -464.3542477, -272.99391166, -283.82902484, -87.82894928, -166.65785636, -29.55839836, -130.26503886, -114.15250293, -163.41153844, 82.98731428, -167.49084336])),
+    ('three_boreholes_network_series_unequal', -0.2, None, 1., 1, np.array([-449.43300626, -545.73676086, -1797.42630407])),
+    ('three_boreholes_network_series_unequal', -0.2, None, np.tile([1., 2., 3., 1.], 3), 4, np.array([-208.47432383, -53.86348191, 100.24050704, -204.09248464, -186.62226798, -86.04634283, 14.21271046, -184.37317225, -504.70429869, -326.58978133, -152.36194671, -480.05212833])),
+    ('three_boreholes_network_series_unequal', -0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.tile([1., 2., 3., 1.], 3), 4, np.array([-103.44758806, -145.80224228, 80.44510753, -152.55896958, -85.17323325, -157.46150512, -19.17059397, -126.38429075, -211.30877858, -492.28316225, -279.11269094, -301.63705046])),
+    ])
+def test_borehole_heat_extraction_rate(
+        network_fixture, m_flow_network, segment_ratios, T_b, nSegments, expected, request):
+    # Extract pipe from fixture
+    network = request.getfixturevalue(network_fixture)
+    # Fluid is propylene-glycol 20%
+    fluid = gt.media.Fluid('MPG', 20.)
+    T_f_in = 5.0            # Inlet fluid temperature [degC]
+    # Borehole heat extraction rates [W]
+    Q_b = network.get_borehole_heat_extraction_rate(
+        T_f_in, T_b, m_flow_network, fluid.cp, nSegments, segment_ratios=segment_ratios)
+    assert np.allclose(Q_b, expected)
+
+
+# Test get_fluid_heat_extraction_rate
+@pytest.mark.parametrize("network_fixture, m_flow_network, segment_ratios, T_b, nSegments, expected", [
+    # Single borehole
+    ('single_borehole_network', 0.2, None, 1., 1, np.array([-1819.4736348927008])),
+    ('single_borehole_network', 0.2, None, np.array([1., 2., 3., 1.]), 4, np.array([-1481.1367317058312])),
+    ('single_borehole_network', 0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 4, np.array([-1299.5107678418537])),
+    ('single_borehole_network', -0.2, None, 1., 1, np.array([-1819.4736348927008])),
+    ('single_borehole_network', -0.2, None, np.array([1., 2., 3., 1.]), 4, np.array([-1481.1367317058312])),
+    ('single_borehole_network', -0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 4, np.array([-1299.5107678418537])),
+    # Ten boreholes (Parallel-connected)
+    ('ten_boreholes_network_rectangular', 2., None, 1., 1, np.full(10, -1819.4736348927008)),
+    ('ten_boreholes_network_rectangular', 2., None, np.tile([1., 2., 3., 1.], 10), 4, np.full(10, -1481.1367317058312)),
+    ('ten_boreholes_network_rectangular', 2., np.array([0.1, 0.35, 0.40, 0.15]), np.tile([1., 2., 3., 1.], 10), 4, np.full(10, -1299.5107678418537)),
+    ('ten_boreholes_network_rectangular', -2., None, 1., 1, np.full(10, -1819.4736348927008)),
+    ('ten_boreholes_network_rectangular', -2., None, np.tile([1., 2., 3., 1.], 10), 4, np.full(10, -1481.1367317058312)),
+    ('ten_boreholes_network_rectangular', -2., np.array([0.1, 0.35, 0.40, 0.15]), np.tile([1., 2., 3., 1.], 10), 4, np.full(10, -1299.5107678418537)),
+    # Three boreholes (Series-connected, unequal lengths)
+    ('three_boreholes_network_series_unequal', 0.2, None, 1., 1, np.array([-1705.73168677, -581.89393685, -504.97044758])),
+    ('three_boreholes_network_series_unequal', 0.2, None, np.tile([1., 2., 3., 1.], 3), 4, np.array([-1388.38932177, -472.52892185, -412.33215297])),
+    ('three_boreholes_network_series_unequal', 0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.tile([1., 2., 3., 1.], 3), 4, np.array([-1218.09977682, -414.31024285, -362.06757046])),
+    ('three_boreholes_network_series_unequal', -0.2, None, 1., 1, np.array([-449.43300626, -545.73676086, -1797.42630407])),
+    ('three_boreholes_network_series_unequal', -0.2, None, np.tile([1., 2., 3., 1.], 3), 4, np.array([-366.18978334, -442.82907261, -1463.70815506])),
+    ('three_boreholes_network_series_unequal', -0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.tile([1., 2., 3., 1.], 3), 4, np.array([-321.3636924, -388.18962309, -1284.34168224])),
+    ])
+def test_fluid_heat_extraction_rate(
+        network_fixture, m_flow_network, segment_ratios, T_b, nSegments, expected, request):
+    # Extract pipe from fixture
+    network = request.getfixturevalue(network_fixture)
+    # Fluid is propylene-glycol 20%
+    fluid = gt.media.Fluid('MPG', 20.)
+    T_f_in = 5.0            # Inlet fluid temperature [degC]
+    # Fluid heat extraction rate [W]
+    Q_f = network.get_fluid_heat_extraction_rate(
+        T_f_in, T_b, m_flow_network, fluid.cp, nSegments, segment_ratios=segment_ratios)
+    assert np.allclose(Q_f, expected)
+
+
+# Test get_network_inlet_temperature
+@pytest.mark.parametrize("network_fixture, m_flow_network, segment_ratios, T_b, nSegments, expected", [
+    # Single borehole
+    ('single_borehole_network', 0.2, None, 1., 1, 7.595314034714041),
+    ('single_borehole_network', 0.2, None, np.array([1., 2., 3., 1.]), 4, 8.33912674339739),
+    ('single_borehole_network', 0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 4, 8.73842016624424),
+    ('single_borehole_network', -0.2, None, 1., 1, 7.595314034714041),
+    ('single_borehole_network', -0.2, None, np.array([1., 2., 3., 1.]), 4, 8.33912674339739),
+    ('single_borehole_network', -0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 4, 8.73842016624424),
+    # Ten boreholes (Parallel-connected)
+    ('ten_boreholes_network_rectangular', 2., None, 1., 1, 7.595314034714041),
+    ('ten_boreholes_network_rectangular', 2., None, np.tile([1., 2., 3., 1.], 10), 4, 8.33912674339739),
+    ('ten_boreholes_network_rectangular', 2., np.array([0.1, 0.35, 0.40, 0.15]), np.tile([1., 2., 3., 1.], 10), 4, 8.73842016624424),
+    ('ten_boreholes_network_rectangular', -2., None, 1., 1, 7.595314034714041),
+    ('ten_boreholes_network_rectangular', -2., None, np.tile([1., 2., 3., 1.], 10), 4, 8.33912674339739),
+    ('ten_boreholes_network_rectangular', -2., np.array([0.1, 0.35, 0.40, 0.15]), np.tile([1., 2., 3., 1.], 10), 4, 8.73842016624424),
+    # Three boreholes (Series-connected, unequal lengths)
+    ('three_boreholes_network_series_unequal', 0.2, None, 1., 1, 13.891230626356375),
+    ('three_boreholes_network_series_unequal', 0.2, None, np.tile([1., 2., 3., 1.], 3), 4, 14.635120055915039),
+    ('three_boreholes_network_series_unequal', 0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.tile([1., 2., 3., 1.], 3), 4, 15.034422782654072),
+    ('three_boreholes_network_series_unequal', -0.2, None, 1., 1, 13.891230626356379),
+    ('three_boreholes_network_series_unequal', -0.2, None, np.tile([1., 2., 3., 1.], 3), 4, 14.635869731958419),
+    ('three_boreholes_network_series_unequal', -0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.tile([1., 2., 3., 1.], 3), 4, 15.03525726407867),
+    ])
+def test_network_inlet_temperature(network_fixture, m_flow_network, segment_ratios, T_b, nSegments, expected, request):
+    # Extract pipe from fixture
+    network = request.getfixturevalue(network_fixture)
+    nBoreholes = len(network.b)
+    # Fluid is propylene-glycol 20%
+    fluid = gt.media.Fluid('MPG', 20.)
+    Q_f = -3000.0 * nBoreholes  # Total heat transfer rate [W]
+    # Inlet fluid temperature [degC]
+    T_f_in = network.get_network_inlet_temperature(
+        Q_f, T_b, m_flow_network, fluid.cp, nSegments, segment_ratios=segment_ratios)
+    assert np.isclose(T_f_in, expected)
+
+
+# Test get_network_outlet_temperature
+@pytest.mark.parametrize("network_fixture, m_flow_network, segment_ratios, T_b, nSegments, expected", [
+    # Single borehole
+    ('single_borehole_network', 0.2, None, 1., 1, 2.712371852688313),
+    ('single_borehole_network', 0.2, None, np.array([1., 2., 3., 1.]), 4, 3.1377635748663573),
+    ('single_borehole_network', 0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 4, 3.3661222931515784),
+    ('single_borehole_network', -0.2, None, 1., 1, 2.712371852688313),
+    ('single_borehole_network', -0.2, None, np.array([1., 2., 3., 1.]), 4, 3.1377635748663573),
+    ('single_borehole_network', -0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 4, 3.3661222931515784),
+    # Ten boreholes (Parallel-connected)
+    ('ten_boreholes_network_rectangular', 2., None, 1., 1, 2.712371852688313),
+    ('ten_boreholes_network_rectangular', 2., None, np.tile([1., 2., 3., 1.], 10), 4, 3.1377635748663573),
+    ('ten_boreholes_network_rectangular', 2., np.array([0.1, 0.35, 0.40, 0.15]), np.tile([1., 2., 3., 1.], 10), 4, 3.3661222931515784),
+    ('ten_boreholes_network_rectangular', -2., None, 1., 1, 2.712371852688313),
+    ('ten_boreholes_network_rectangular', -2., None, np.tile([1., 2., 3., 1.], 10), 4, 3.1377635748663573),
+    ('ten_boreholes_network_rectangular', -2., np.array([0.1, 0.35, 0.40, 0.15]), np.tile([1., 2., 3., 1.], 10), 4, 3.3661222931515784),
+    # Three boreholes (Series-connected, unequal lengths)
+    ('three_boreholes_network_series_unequal', 0.2, None, 1., 1, 1.4888629029742049),
+    ('three_boreholes_network_series_unequal', 0.2, None, np.tile([1., 2., 3., 1.], 3), 4, 2.141837346026401),
+    ('three_boreholes_network_series_unequal', 0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.tile([1., 2., 3., 1.], 3), 4, 2.492339000225644),
+    ('three_boreholes_network_series_unequal', -0.2, None, 1., 1, 1.4888629029742049),
+    ('three_boreholes_network_series_unequal', -0.2, None, np.tile([1., 2., 3., 1.], 3), 4, 2.142495399868063),
+    ('three_boreholes_network_series_unequal', -0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.tile([1., 2., 3., 1.], 3), 4, 2.4930714948973103),
+    ])
+def test_network_outlet_temperature(
+        network_fixture, m_flow_network, segment_ratios, T_b, nSegments, expected, request):
+    # Extract pipe from fixture
+    network = request.getfixturevalue(network_fixture)
+    # Fluid is propylene-glycol 20%
+    fluid = gt.media.Fluid('MPG', 20.)
+    T_f_in = 5.0            # Inlet fluid temperature [degC]
+    # Outlet fluid temperature [degC]
+    T_f_out = network.get_network_outlet_temperature(
+        T_f_in, T_b, m_flow_network, fluid.cp, nSegments, segment_ratios=segment_ratios)
+    assert np.isclose(T_f_out, expected)
+
+
+# Test get_network_heat_extraction_rate
+@pytest.mark.parametrize("network_fixture, m_flow_network, segment_ratios, T_b, nSegments, expected", [
+    # Single borehole
+    ('single_borehole_network', 0.2, None, 1., 1, -1819.4736348927008),
+    ('single_borehole_network', 0.2, None, np.array([1., 2., 3., 1.]), 4, -1481.1367317058312),
+    ('single_borehole_network', 0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 4, -1299.5107678418537),
+    ('single_borehole_network', -0.2, None, 1., 1, -1819.4736348927008),
+    ('single_borehole_network', -0.2, None, np.array([1., 2., 3., 1.]), 4, -1481.1367317058312),
+    ('single_borehole_network', -0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 4, -1299.5107678418537),
+    # Ten boreholes (Parallel-connected)
+    ('ten_boreholes_network_rectangular', 2., None, 1., 1, -18194.736348927008),
+    ('ten_boreholes_network_rectangular', 2., None, np.tile([1., 2., 3., 1.], 10), 4, -14811.367317058312),
+    ('ten_boreholes_network_rectangular', 2., np.array([0.1, 0.35, 0.40, 0.15]), np.tile([1., 2., 3., 1.], 10), 4, -12995.107678418537),
+    ('ten_boreholes_network_rectangular', -2., None, 1., 1, -18194.736348927008),
+    ('ten_boreholes_network_rectangular', -2., None, np.tile([1., 2., 3., 1.], 10), 4, -14811.367317058312),
+    ('ten_boreholes_network_rectangular', -2., np.array([0.1, 0.35, 0.40, 0.15]), np.tile([1., 2., 3., 1.], 10), 4, -12995.107678418537),
+    # Three boreholes (Series-connected, unequal lengths)
+    ('three_boreholes_network_series_unequal', 0.2, None, 1., 1, -2792.5960711925586),
+    ('three_boreholes_network_series_unequal', 0.2, None, np.tile([1., 2., 3., 1.], 3), 4, -2273.2503965957585),
+    ('three_boreholes_network_series_unequal', 0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.tile([1., 2., 3., 1.], 3), 4, -1994.477590118784),
+    ('three_boreholes_network_series_unequal', -0.2, None, 1., 1, -2792.596071192558),
+    ('three_boreholes_network_series_unequal', -0.2, None, np.tile([1., 2., 3., 1.], 3), 4, -2272.7270110024083),
+    ('three_boreholes_network_series_unequal', -0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.tile([1., 2., 3., 1.], 3), 4, -1993.8949977318334),
+    ])
+def test_total_heat_extraction_rate(
+        network_fixture, m_flow_network, segment_ratios, T_b, nSegments, expected, request):
+    # Extract pipe from fixture
+    network = request.getfixturevalue(network_fixture)
+    # Fluid is propylene-glycol 20%
+    fluid = gt.media.Fluid('MPG', 20.)
+    T_f_in = 5.0            # Inlet fluid temperature [degC]
+    # Total heat extraction rate [W]
+    Q_t = network.get_network_heat_extraction_rate(
+        T_f_in, T_b, m_flow_network, fluid.cp, nSegments, segment_ratios=segment_ratios)
+    assert np.isclose(Q_t, expected)
diff --git a/tests/pipes_test.py b/tests/pipes_test.py
index 67328b1..9d35f87 100644
--- a/tests/pipes_test.py
+++ b/tests/pipes_test.py
@@ -13,9 +13,12 @@
 # =============================================================================
 # Test convective_heat_transfer_coefficient_circular_pipe
 @pytest.mark.parametrize("m_flow, expected", [
-    (0.05, 90.07260000000001),  # Laminar flow
-    (0.10, 572.256944167273),   # Transition flow
-    (0.50, 4036.196217814895)   # Turbulent flow
+    (0.05, 90.07260000000001),      # Laminar flow
+    (0.10, 572.256944167273),       # Transition flow
+    (0.50, 4036.196217814895),      # Turbulent flow
+    (-0.05, 90.07260000000001),     # Laminar flow
+    (-0.10, 572.256944167273),      # Transition flow
+    (-0.50, 4036.196217814895),     # Turbulent flow
     ])
 def test_convective_heat_transfer_coefficient_circular_pipe(m_flow, expected):
     r_in = 0.01         # Inner radius [m]
@@ -34,7 +37,10 @@ def test_convective_heat_transfer_coefficient_circular_pipe(m_flow, expected):
 @pytest.mark.parametrize("m_flow, expected", [
     (0.05, (141.4907984705223, 110.95487746200112)),    # Laminar flow
     (0.40, (904.4869811625874, 904.4869811625874)),     # Transition flow
-    (0.60, (1411.2063074288633, 1411.2063074288633))    # Turbulent flow
+    (0.60, (1411.2063074288633, 1411.2063074288633)),    # Turbulent flow
+    (-0.05, (141.4907984705223, 110.95487746200112)),    # Laminar flow
+    (-0.40, (904.4869811625874, 904.4869811625874)),     # Transition flow
+    (-0.60, (1411.2063074288633, 1411.2063074288633)),    # Turbulent flow
     ])
 def test_convective_heat_transfer_coefficient_concentric_annulus(
         m_flow, expected):
@@ -65,7 +71,9 @@ def test_conduction_thermal_resistance_circular_pipe():
 # Test fluid_friction_factor_circular_pipe
 @pytest.mark.parametrize("m_flow, expected", [
     (0.05, 0.04081718025087723),    # Laminar flow
-    (0.50, 0.027641340780182006)    # Turbulent flow
+    (0.50, 0.027641340780182006),   # Turbulent flow
+    (-0.05, 0.04081718025087723),   # Laminar flow
+    (-0.50, 0.027641340780182006),  # Turbulent flow
     ])
 def test_fluid_friction_factor_circular_pipe(m_flow, expected):
     r_in = 0.01         # Inner radius [m]
@@ -148,57 +156,118 @@ def test_multipole(J, expected):
 # Test pipe classes
 # =============================================================================
 # Test get_temperature
-@pytest.mark.parametrize("pipe_fixture, segment_ratios, T_b, z, expected", [
+@pytest.mark.parametrize("pipe_fixture, m_flow_borehole, segment_ratios, T_b, z, expected", [
     # Single U-tube
-    ('single_Utube', None, 1., 65., np.array([4.34676755, 3.07354134])),
-    ('single_Utube', None, np.array([1., 2., 3., 1.]), 65., np.array([4.41754093, 3.49949295])),
-    ('single_Utube', np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 65., np.array([4.47282344, 3.70132653])),
-    ('single_Utube', None, 1., np.array([65., 75., 150.]), np.array([[4.34676755, 3.07354134], [4.25566624, 3.13435325], [3.64199359, 3.64199359]])),
-    ('single_Utube', None, np.array([1., 2., 3., 1.]), np.array([65., 75., 150.]), np.array([[4.41754093, 3.49949295], [4.35173147, 3.54346564], [3.89942492, 3.89942492]])),
-    ('single_Utube', np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), np.array([65., 75., 150.]), np.array([[4.47282344, 3.70132653], [4.4261555, 3.72955733], [4.03655472, 4.03655472]])),
+    ('single_Utube', 0.2, None, 1., 65., np.array([4.34676755, 3.07354134])),
+    ('single_Utube', 0.2, None, np.array([1., 2., 3., 1.]), 65., np.array([4.41754093, 3.49949295])),
+    ('single_Utube', 0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 65., np.array([4.47282344, 3.70132653])),
+    ('single_Utube', 0.2, None, 1., np.array([65., 75., 150.]), np.array([[4.34676755, 3.07354134], [4.25566624, 3.13435325], [3.64199359, 3.64199359]])),
+    ('single_Utube', 0.2, None, np.array([1., 2., 3., 1.]), np.array([65., 75., 150.]), np.array([[4.41754093, 3.49949295], [4.35173147, 3.54346564], [3.89942492, 3.89942492]])),
+    ('single_Utube', 0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), np.array([65., 75., 150.]), np.array([[4.47282344, 3.70132653], [4.4261555, 3.72955733], [4.03655472, 4.03655472]])),
+    ('single_Utube', -0.2, None, 1., 65., np.array([3.07354134, 4.34676755])),
+    ('single_Utube', -0.2, None, np.array([1., 2., 3., 1.]), 65., np.array([3.49949295, 4.41754093])),
+    ('single_Utube', -0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 65., np.array([3.70132653, 4.47282344])),
+    ('single_Utube', -0.2, None, 1., np.array([65., 75., 150.]), np.array([[3.07354134, 4.34676755], [3.13435325, 4.25566624], [3.64199359, 3.64199359]])),
+    ('single_Utube', -0.2, None, np.array([1., 2., 3., 1.]), np.array([65., 75., 150.]), np.array([[3.49949295, 4.41754093], [3.54346564, 4.35173147], [3.89942492, 3.89942492]])),
+    ('single_Utube', -0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), np.array([65., 75., 150.]), np.array([[3.70132653, 4.47282344], [3.72955733, 4.4261555], [4.03655472, 4.03655472]])),
     # Double U-tube (Parallel)
-    ('double_Utube_parallel', None, 1., 65., np.array([3.87525104, 3.87525104, 2.20313908, 2.20313908])),
-    ('double_Utube_parallel', None, np.array([1., 2., 3., 1.]), 65., np.array([4.00464852, 4.00464852, 2.84788608, 2.84788608])),
-    ('double_Utube_parallel', np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 65., np.array([4.09844256, 4.09844256, 3.14231047, 3.14231047])),
-    ('double_Utube_parallel', None, 1., np.array([65., 75., 150.]), np.array([[3.87525104, 3.87525104, 2.20313908, 2.20313908], [3.73265141, 3.73265141, 2.26719823, 2.26719823], [2.86396102, 2.86396102, 2.86396102, 2.86396102]])),
-    ('double_Utube_parallel', None, np.array([1., 2., 3., 1.]), np.array([65., 75., 150.]), np.array([[4.00464852, 4.00464852, 2.84788608, 2.84788608], [3.90522192, 3.90522192, 2.89301847, 2.89301847], [3.26970513, 3.26970513, 3.26970513, 3.26970513]])),
-    ('double_Utube_parallel', np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), np.array([65., 75., 150.]), np.array([[4.09844256, 4.09844256, 3.14231047, 3.14231047], [4.03189362, 4.03189362, 3.16367201, 3.16367201], [3.48739254, 3.48739254, 3.48739254, 3.48739254]])),
+    ('double_Utube_parallel', 0.2, None, 1., 65., np.array([3.87525104, 3.87525104, 2.20313908, 2.20313908])),
+    ('double_Utube_parallel', 0.2, None, np.array([1., 2., 3., 1.]), 65., np.array([4.00464852, 4.00464852, 2.84788608, 2.84788608])),
+    ('double_Utube_parallel', 0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 65., np.array([4.09844256, 4.09844256, 3.14231047, 3.14231047])),
+    ('double_Utube_parallel', 0.2, None, 1., np.array([65., 75., 150.]), np.array([[3.87525104, 3.87525104, 2.20313908, 2.20313908], [3.73265141, 3.73265141, 2.26719823, 2.26719823], [2.86396102, 2.86396102, 2.86396102, 2.86396102]])),
+    ('double_Utube_parallel', 0.2, None, np.array([1., 2., 3., 1.]), np.array([65., 75., 150.]), np.array([[4.00464852, 4.00464852, 2.84788608, 2.84788608], [3.90522192, 3.90522192, 2.89301847, 2.89301847], [3.26970513, 3.26970513, 3.26970513, 3.26970513]])),
+    ('double_Utube_parallel', 0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), np.array([65., 75., 150.]), np.array([[4.09844256, 4.09844256, 3.14231047, 3.14231047], [4.03189362, 4.03189362, 3.16367201, 3.16367201], [3.48739254, 3.48739254, 3.48739254, 3.48739254]])),
+    ('double_Utube_parallel', -0.2, None, 1., 65., np.array([2.20313908, 2.20313908, 3.87525104, 3.87525104])),
+    ('double_Utube_parallel', -0.2, None, np.array([1., 2., 3., 1.]), 65., np.array([2.84788608, 2.84788608, 4.00464852, 4.00464852])),
+    ('double_Utube_parallel', -0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 65., np.array([3.14231047, 3.14231047, 4.09844256, 4.09844256])),
+    ('double_Utube_parallel', -0.2, None, 1., np.array([65., 75., 150.]), np.array([[2.20313908, 2.20313908, 3.87525104, 3.87525104], [2.26719823, 2.26719823, 3.73265141, 3.73265141], [2.86396102, 2.86396102, 2.86396102, 2.86396102]])),
+    ('double_Utube_parallel', -0.2, None, np.array([1., 2., 3., 1.]), np.array([65., 75., 150.]), np.array([[2.84788608, 2.84788608, 4.00464852, 4.00464852], [2.89301847, 2.89301847, 3.90522192, 3.90522192], [3.26970513, 3.26970513, 3.26970513, 3.26970513]])),
+    ('double_Utube_parallel', -0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), np.array([65., 75., 150.]), np.array([[3.14231047, 3.14231047, 4.09844256, 4.09844256], [3.16367201, 3.16367201, 4.03189362, 4.03189362], [3.48739254, 3.48739254, 3.48739254, 3.48739254]])),
+    # Double U-tube (Parallel, Asymmetrical)
+    ('double_Utube_parallel_asymmetrical', 0.2, None, 1., 65., np.array([3.88377937, 3.98193395, 2.38708007, 2.47953756])),
+    ('double_Utube_parallel_asymmetrical', 0.2, None, np.array([1., 2., 3., 1.]), 65., np.array([4.01466532, 4.11668145, 2.97676364, 3.05353369])),
+    ('double_Utube_parallel_asymmetrical', 0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 65., np.array([4.10803611, 4.20442544, 3.25024389, 3.31798041])),
+    ('double_Utube_parallel_asymmetrical', 0.2, None, 1., np.array([65., 75., 150.]), np.array([[3.88377937, 3.98193395, 2.38708007, 2.47953756], [3.74282975, 3.85383443, 2.43329782, 2.53690382], [2.88927312, 3.08795748, 2.88927312, 3.08795748]])),
+    ('double_Utube_parallel_asymmetrical', 0.2, None, np.array([1., 2., 3., 1.]), np.array([65., 75., 150.]), np.array([[4.01466532, 4.11668145, 2.97676364, 3.05353369], [3.91631953, 4.02560538, 3.00903257, 3.0945893], [3.29041317, 3.45243657, 3.29041317, 3.45243657]])),
+    ('double_Utube_parallel_asymmetrical', 0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), np.array([65., 75., 150.]), np.array([[4.10803611, 4.20442544, 3.25024389, 3.31798041], [4.04195626, 4.14006795, 3.26395583, 3.33909292], [3.50562748, 3.64753152, 3.50562748, 3.64753152]])),
+    ('double_Utube_parallel_asymmetrical', -0.2, None, 1., 65., np.array([2.35721262, 2.55271346, 3.97798837, 3.98287957])),
+    ('double_Utube_parallel_asymmetrical', -0.2, None, np.array([1., 2., 3., 1.]), 65., np.array([2.97085504, 3.1050909, 4.10590079, 4.10778519])),
+    ('double_Utube_parallel_asymmetrical', -0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 65., np.array([3.24913725, 3.36045033, 4.19291661, 4.1940752])),
+    ('double_Utube_parallel_asymmetrical', -0.2, None, 1., np.array([65., 75., 150.]), np.array([[2.35721262, 2.55271346, 3.97798837, 3.98287957], [2.42477992, 2.59763795, 3.84829923, 3.85304893], [3.06304785, 3.06132252, 3.06304785, 3.06132252]])),
+    ('double_Utube_parallel_asymmetrical', -0.2, None, np.array([1., 2., 3., 1.]), np.array([65., 75., 150.]), np.array([[2.97085504, 3.1050909, 4.10590079, 4.10778519], [3.01922892, 3.13720999, 4.01422481, 4.01610335], [3.43184683, 3.43038357, 3.43184683, 3.43038357]])),
+    ('double_Utube_parallel_asymmetrical', -0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), np.array([65., 75., 150.]), np.array([[3.24913725, 3.36045033, 4.19291661, 4.1940752], [3.27383723, 3.37530358, 4.12926985, 4.13067701], [3.62926994, 3.62789036, 3.62926994, 3.62789036]])),
     # Double U-tube (Series)
-    ('double_Utube_series', None, 1., 65., np.array([4.36908096, 2.53231146, 3.13441957, 2.03763963])),
-    ('double_Utube_series', None, np.array([1., 2., 3., 1.]), 65., np.array([4.44022419, 2.94528677, 3.54323578, 2.65057213])),
-    ('double_Utube_series', np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 65., np.array([4.49417249, 3.18261918, 3.73818275, 2.95502223])),
-    ('double_Utube_series', None, 1., np.array([65., 75., 150.]), np.array([[4.36908096, 2.53231146, 3.13441957, 2.03763963], [4.28094228, 2.49706752, 3.19348974, 2.0612353], [3.68632591, 2.25874399, 3.68632591, 2.25874399]])),
-    ('double_Utube_series', None, np.array([1., 2., 3., 1.]), np.array([65., 75., 150.]), np.array([[4.44022419, 2.94528677, 3.54323578, 2.65057213], [4.37608674, 2.92420008, 3.58625745, 2.66472128], [3.93523525, 2.77333354, 3.93523525, 2.77333354]])),
-    ('double_Utube_series', np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), np.array([65., 75., 150.]), np.array([[4.49417249, 3.18261918, 3.73818275, 2.95502223], [4.44797829, 3.17419469, 3.76650271, 2.95801146], [4.06780952, 3.04843839, 4.06780952, 3.04843839]])),
+    ('double_Utube_series', 0.2, None, 1., 65., np.array([4.36908096, 2.53231146, 3.13441957, 2.03763963])),
+    ('double_Utube_series', 0.2, None, np.array([1., 2., 3., 1.]), 65., np.array([4.44022419, 2.94528677, 3.54323578, 2.65057213])),
+    ('double_Utube_series', 0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 65., np.array([4.49417249, 3.18261918, 3.73818275, 2.95502223])),
+    ('double_Utube_series', 0.2, None, 1., np.array([65., 75., 150.]), np.array([[4.36908096, 2.53231146, 3.13441957, 2.03763963], [4.28094228, 2.49706752, 3.19348974, 2.0612353], [3.68632591, 2.25874399, 3.68632591, 2.25874399]])),
+    ('double_Utube_series', 0.2, None, np.array([1., 2., 3., 1.]), np.array([65., 75., 150.]), np.array([[4.44022419, 2.94528677, 3.54323578, 2.65057213], [4.37608674, 2.92420008, 3.58625745, 2.66472128], [3.93523525, 2.77333354, 3.93523525, 2.77333354]])),
+    ('double_Utube_series', 0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), np.array([65., 75., 150.]), np.array([[4.49417249, 3.18261918, 3.73818275, 2.95502223], [4.44797829, 3.17419469, 3.76650271, 2.95801146], [4.06780952, 3.04843839, 4.06780952, 3.04843839]])),
+    ('double_Utube_series', -0.2, None, 1., 65., np.array([2.03763963, 3.13441957, 2.53231146, 4.36908096])),
+    ('double_Utube_series', -0.2, None, np.array([1., 2., 3., 1.]), 65., np.array([2.65057213, 3.54323578, 2.94528677, 4.44022419])),
+    ('double_Utube_series', -0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 65., np.array([2.95502223, 3.73818275, 3.18261918, 4.49417249])),
+    ('double_Utube_series', -0.2, None, 1., np.array([65., 75., 150.]), np.array([[2.03763963, 3.13441957, 2.53231146, 4.36908096], [2.0612353, 3.19348974, 2.49706752, 4.28094228], [2.25874399, 3.68632591, 2.25874399, 3.68632591]])),
+    ('double_Utube_series', -0.2, None, np.array([1., 2., 3., 1.]), np.array([65., 75., 150.]), np.array([[2.65057213, 3.54323578, 2.94528677, 4.44022419], [2.66472128, 3.58625745, 2.92420008, 4.37608674], [2.77333354, 3.93523525, 2.77333354, 3.93523525]])),
+    ('double_Utube_series', -0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), np.array([65., 75., 150.]), np.array([[2.95502223, 3.73818275, 3.18261918, 4.49417249], [2.95801146, 3.76650271, 3.17419469, 4.44797829], [3.04843839, 4.06780952, 3.04843839, 4.06780952]])),
+    # Double U-tube (Series, Asymmetrical)
+    ('double_Utube_series_asymmetrical', 0.2, None, 1., 65., np.array([4.37509052, 2.69495427, 3.20649518, 2.28963469])),
+    ('double_Utube_series_asymmetrical', 0.2, None, np.array([1., 2., 3., 1.]), 65., np.array([4.44658902, 3.09141515, 3.59118747, 2.84576981])),
+    ('double_Utube_series_asymmetrical', 0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 65., np.array([4.50016074, 3.31495225, 3.77742245, 3.12368927])),
+    ('double_Utube_series_asymmetrical', 0.2, None, 1., np.array([65., 75., 150.]), np.array([[4.37509052, 2.69495427, 3.20649518, 2.28963469], [4.2879064, 2.6675356, 3.25915038, 2.31078176], [3.70074184, 2.48937779, 3.70074184, 2.48937779]])),
+    ('double_Utube_series_asymmetrical', 0.2, None, np.array([1., 2., 3., 1.]), np.array([65., 75., 150.]), np.array([[4.44658902, 3.09141515, 3.59118747, 2.84576981], [4.38308285, 3.07494832, 3.63005149, 2.85887109], [3.94700606, 2.96119851, 3.94700606, 2.96119851]])),
+    ('double_Utube_series_asymmetrical', 0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), np.array([65., 75., 150.]), np.array([[4.50016074, 3.31495225, 3.77742245, 3.12368927], [4.45427774, 3.3080315, 3.80394222, 3.12746125], [4.07815937, 3.21327688, 4.07815937, 3.21327688]])),
+    ('double_Utube_series_asymmetrical', -0.2, None, 1., 65., np.array([2.33841354, 3.35868363, 2.85292145, 4.43249565])),
+    ('double_Utube_series_asymmetrical', -0.2, None, np.array([1., 2., 3., 1.]), 65., np.array([2.89383516, 3.71125236, 3.21591736, 4.50047147])),
+    ('double_Utube_series_asymmetrical', -0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 65., np.array([3.16788193, 3.88151739, 3.42296055, 4.54969912])),
+    ('double_Utube_series_asymmetrical', -0.2, None, 1., np.array([65., 75., 150.]), np.array([[2.33841354, 3.35868363, 2.85292145, 4.43249565], [2.36782915, 3.40736985, 2.82202214, 4.35330199], [2.61655429, 3.82066305, 2.61655429, 3.82066305]])),
+    ('double_Utube_series_asymmetrical', -0.2, None, np.array([1., 2., 3., 1.]), np.array([65., 75., 150.]), np.array([[2.89383516, 3.71125236, 3.21591736, 4.50047147], [2.91289702, 3.74732351, 3.19700954, 4.44232958], [3.06431085, 4.04418788, 3.06431085, 4.04418788]])),
+    ('double_Utube_series_asymmetrical', -0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), np.array([65., 75., 150.]), np.array([[3.16788193, 3.88151739, 3.42296055, 4.54969912], [3.17542873, 3.90622357, 3.41452794, 4.5069777], [3.30353387, 4.16317605, 3.30353387, 4.16317605]])),
     # Coaxial (Annular pipe is inlet pipe)
-    ('coaxial_annular_in', None, 1., 65., np.array([3.15203088, 2.18408362])),
-    ('coaxial_annular_in', None, np.array([1., 2., 3., 1.]), 65., np.array([3.4176666 , 2.73205968])),
-    ('coaxial_annular_in', np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 65., np.array([3.57440691, 3.0183749])),
-    ('coaxial_annular_in', None, 1., np.array([65., 75., 150.]), np.array([[3.15203088, 2.18408362], [2.96401382, 2.15051705], [2.0474566, 2.0474566]])),
-    ('coaxial_annular_in', None, np.array([1., 2., 3., 1.]), np.array([65., 75., 150.]), np.array([[3.4176666, 2.73205968], [3.2920645, 2.7081367], [2.59520774, 2.59520774]])),
-    ('coaxial_annular_in', np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), np.array([65., 75., 150.]), np.array([[3.57440691, 3.0183749], [3.49639377, 2.99877788], [2.90174463, 2.90174463]])),
+    ('coaxial_annular_in', 0.2, None, 1., 65., np.array([3.15203088, 2.18408362])),
+    ('coaxial_annular_in', 0.2, None, np.array([1., 2., 3., 1.]), 65., np.array([3.4176666, 2.73205968])),
+    ('coaxial_annular_in', 0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 65., np.array([3.57440691, 3.0183749])),
+    ('coaxial_annular_in', 0.2, None, 1., np.array([65., 75., 150.]), np.array([[3.15203088, 2.18408362], [2.96401382, 2.15051705], [2.0474566, 2.0474566]])),
+    ('coaxial_annular_in', 0.2, None, np.array([1., 2., 3., 1.]), np.array([65., 75., 150.]), np.array([[3.4176666, 2.73205968], [3.2920645, 2.7081367], [2.59520774, 2.59520774]])),
+    ('coaxial_annular_in', 0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), np.array([65., 75., 150.]), np.array([[3.57440691, 3.0183749], [3.49639377, 2.99877788], [2.90174463, 2.90174463]])),
+    ('coaxial_annular_in', -0.2, None, 1., 65., np.array([2.92933532, 4.50649998])),
+    ('coaxial_annular_in', -0.2, None, np.array([1., 2., 3., 1.]), 65., np.array([3.50307539, 4.62416027])),
+    ('coaxial_annular_in', -0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 65., np.array([3.72110166, 4.67083279])),
+    ('coaxial_annular_in', -0.2, None, 1., np.array([65., 75., 150.]), np.array([[2.92933532, 4.50649998], [3.02086677, 4.44976224], [4.23165721, 4.23165721]])),
+    ('coaxial_annular_in', -0.2, None, np.array([1., 2., 3., 1.]), np.array([65., 75., 150.]), np.array([[3.50307539, 4.62416027], [3.57860389, 4.58402116], [4.39160319, 4.39160319]])),
+    ('coaxial_annular_in', -0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), np.array([65., 75., 150.]), np.array([[3.72110166, 4.67083279], [3.76190344, 4.63666726], [4.47159269, 4.47159269]])),
     # Coaxial (Annular pipe is outlet pipe)
-    ('coaxial_annular_out', None, 1., 65., np.array([4.50649998, 2.92933532])),
-    ('coaxial_annular_out', None, np.array([1., 2., 3., 1.]), 65., np.array([4.62416027, 3.50307539])),
-    ('coaxial_annular_out', np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 65., np.array([4.67083279, 3.72110166])),
-    ('coaxial_annular_out', None, 1., np.array([65., 75., 150.]), np.array([[4.50649998, 2.92933532], [4.44976224, 3.02086677], [4.23165721, 4.23165721]])),
-    ('coaxial_annular_out', None, np.array([1., 2., 3., 1.]), np.array([65., 75., 150.]), np.array([[4.62416027, 3.50307539], [4.58402116, 3.57860389], [4.39160319, 4.39160319]])),
-    ('coaxial_annular_out', np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), np.array([65., 75., 150.]), np.array([[4.67083279, 3.72110166], [4.63666726, 3.76190344], [4.47159269, 4.47159269]])),
+    ('coaxial_annular_out', 0.2, None, 1., 65., np.array([4.50649998, 2.92933532])),
+    ('coaxial_annular_out', 0.2, None, np.array([1., 2., 3., 1.]), 65., np.array([4.62416027, 3.50307539])),
+    ('coaxial_annular_out', 0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 65., np.array([4.67083279, 3.72110166])),
+    ('coaxial_annular_out', 0.2, None, 1., np.array([65., 75., 150.]), np.array([[4.50649998, 2.92933532], [4.44976224, 3.02086677], [4.23165721, 4.23165721]])),
+    ('coaxial_annular_out', 0.2, None, np.array([1., 2., 3., 1.]), np.array([65., 75., 150.]), np.array([[4.62416027, 3.50307539], [4.58402116, 3.57860389], [4.39160319, 4.39160319]])),
+    ('coaxial_annular_out', 0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), np.array([65., 75., 150.]), np.array([[4.67083279, 3.72110166], [4.63666726, 3.76190344], [4.47159269, 4.47159269]])),
+    ('coaxial_annular_out', -0.2, None, 1., 65., np.array([2.18408362, 3.15203088])),
+    ('coaxial_annular_out', -0.2, None, np.array([1., 2., 3., 1.]), 65., np.array([2.73205968, 3.4176666])),
+    ('coaxial_annular_out', -0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 65., np.array([3.0183749, 3.57440691])),
+    ('coaxial_annular_out', -0.2, None, 1., np.array([65., 75., 150.]), np.array([[2.18408362, 3.15203088], [2.15051705, 2.96401382], [2.0474566, 2.0474566]])),
+    ('coaxial_annular_out', -0.2, None, np.array([1., 2., 3., 1.]), np.array([65., 75., 150.]), np.array([[2.73205968, 3.4176666], [2.7081367, 3.2920645], [2.59520774, 2.59520774]])),
+    ('coaxial_annular_out', -0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), np.array([65., 75., 150.]), np.array([[3.0183749, 3.57440691], [2.99877788, 3.49639377], [2.90174463, 2.90174463]])),
     # Coaxial (Annular pipe is inlet pipe, no grout)
-    ('coaxial_no_grout', None, 1., 140., np.array([2.8355954, 1.48415788])),
-    ('coaxial_no_grout', None, np.array([1., 2., 3., 1.]), 140., np.array([3.2799902, 1.60137175])),
-    ('coaxial_no_grout', np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 140., np.array([3.51716507, 1.66392916])),
-    ('coaxial_no_grout', None, 1., np.array([65., 140., 150.]), np.array([[3.77463898, 1.10936191], [2.8355954, 1.48415788], [2.78756714, 2.78756714]])),
-    ('coaxial_no_grout', None, np.array([1., 2., 3., 1.]), np.array([65., 140., 150.]), np.array([[3.95937381, 2.30543206], [3.2799902, 1.60137175], [3.22033437, 3.22033437]])),
-    ('coaxial_no_grout', np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), np.array([65., 140., 150.]), np.array([[4.07193759, 2.7536371], [3.51716507, 1.66392916], [3.45130357, 3.45130357]])),
+    ('coaxial_no_grout', 0.2, None, 1., 140., np.array([2.8355954, 1.48415788])),
+    ('coaxial_no_grout', 0.2, None, np.array([1., 2., 3., 1.]), 140., np.array([3.2799902, 1.60137175])),
+    ('coaxial_no_grout', 0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 140., np.array([3.51716507, 1.66392916])),
+    ('coaxial_no_grout', 0.2, None, 1., np.array([65., 140., 150.]), np.array([[3.77463898, 1.10936191], [2.8355954, 1.48415788], [2.78756714, 2.78756714]])),
+    ('coaxial_no_grout', 0.2, None, np.array([1., 2., 3., 1.]), np.array([65., 140., 150.]), np.array([[3.95937381, 2.30543206], [3.2799902, 1.60137175], [3.22033437, 3.22033437]])),
+    ('coaxial_no_grout', 0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), np.array([65., 140., 150.]), np.array([[4.07193759, 2.7536371], [3.51716507, 1.66392916], [3.45130357, 3.45130357]])),
+    ('coaxial_no_grout', -0.2, None, 1., 140., np.array([1., 1.00000001])),
+    ('coaxial_no_grout', -0.2, None, np.array([1., 2., 3., 1.]), 140., np.array([1.01006961, 1.03817708])),
+    ('coaxial_no_grout', -0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 140., np.array([1.0859313, 1.3257927])),
+    ('coaxial_no_grout', -0.2, None, 1., np.array([65., 140., 150.]), np.array([[1.00001469, 1.00037275], [1., 1.00000001], [1., 1.]])),
+    ('coaxial_no_grout', -0.2, None, np.array([1., 2., 3., 1.]), np.array([65., 140., 150.]), np.array([[1.44148538, 1.95942376], [1.01006961, 1.03817708], [1.00943254, 1.00943254]])),
+    ('coaxial_no_grout', -0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), np.array([65., 140., 150.]), np.array([[1.6430022 , 1.98553786], [1.0859313 , 1.3257927], [1.08049473, 1.08049473]])),
     ])
 def test_temperature(
-        pipe_fixture, segment_ratios, T_b, z, expected, request):
+        pipe_fixture, m_flow_borehole, segment_ratios, T_b, z, expected, request):
     # Extract pipe from fixture
     pipe = request.getfixturevalue(pipe_fixture)
     # Fluid is propylene-glycol 20%
     fluid = gt.media.Fluid('MPG', 20.)
-    m_flow_borehole = 0.2   # Fluid mass flow rate [kg/s]
     T_f_in = 5.0            # Inlet fluid temperature [degC]
     # Fluid temperatures [degC]
     T_f = pipe.get_temperature(
@@ -208,39 +277,70 @@ def test_temperature(
 
 
 # Test get_outlet_temperature
-@pytest.mark.parametrize("pipe_fixture, segment_ratios, T_b, expected", [
+@pytest.mark.parametrize("pipe_fixture, m_flow_borehole, segment_ratios, T_b, expected", [
     # Single U-tube
-    ('single_Utube', None, 1., 2.712371852688313),
-    ('single_Utube', None, np.array([1., 2., 3., 1.]), 3.1377635748663573),
-    ('single_Utube', np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 3.3661222931515784),
+    ('single_Utube', 0.2, None, 1., 2.712371852688313),
+    ('single_Utube', 0.2, None, np.array([1., 2., 3., 1.]), 3.1377635748663573),
+    ('single_Utube', 0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 3.3661222931515784),
+    ('single_Utube', -0.2, None, 1., 2.712371852688313),
+    ('single_Utube', -0.2, None, np.array([1., 2., 3., 1.]), 3.1377635748663573),
+    ('single_Utube', -0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 3.3661222931515784),
     # Double U-tube (Parallel)
-    ('double_Utube_parallel', None, 1., 1.8553031331306218),
-    ('double_Utube_parallel', None, np.array([1., 2., 3., 1.]), 2.4278457017624655),
-    ('double_Utube_parallel', np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 2.7404064201268588),
+    ('double_Utube_parallel', 0.2, None, 1., 1.8553031331306218),
+    ('double_Utube_parallel', 0.2, None, np.array([1., 2., 3., 1.]), 2.4278457017624655),
+    ('double_Utube_parallel', 0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 2.7404064201268588),
+    ('double_Utube_parallel', -0.2, None, 1., 1.8553031331306218),
+    ('double_Utube_parallel', -0.2, None, np.array([1., 2., 3., 1.]), 2.4278457017624655),
+    ('double_Utube_parallel', -0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 2.7404064201268588),
+    # Double U-tube (Parallel, Asymmetrical)
+    ('double_Utube_parallel_asymmetrical', 0.2, None, 1., 2.1704888954025483),
+    ('double_Utube_parallel_asymmetrical', 0.2, None, np.array([1., 2., 3., 1.]), 2.6878415147516206),
+    ('double_Utube_parallel_asymmetrical', 0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 2.969467601101167),
+    ('double_Utube_parallel_asymmetrical', -0.2, None, 1., 2.1704889162427445),
+    ('double_Utube_parallel_asymmetrical', -0.2, None, np.array([1., 2., 3., 1.]), 2.6859171311224035),
+    ('double_Utube_parallel_asymmetrical', -0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 2.9670350185625898),
     # Double U-tube (Series)
-    ('double_Utube_series', None, 1., 1.8983711735742066),
-    ('double_Utube_series', None, np.array([1., 2., 3., 1.]), 2.4755999700741578),
-    ('double_Utube_series', np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 2.7852683378456216),
+    ('double_Utube_series', 0.2, None, 1., 1.8983711735742066),
+    ('double_Utube_series', 0.2, None, np.array([1., 2., 3., 1.]), 2.4755999700741578),
+    ('double_Utube_series', 0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 2.7852683378456216),
+    ('double_Utube_series', -0.2, None, 1., 1.8983711735742066),
+    ('double_Utube_series', -0.2, None, np.array([1., 2., 3., 1.]), 2.4755999700741578),
+    ('double_Utube_series', -0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 2.7852683378456216),
+    # Double U-tube (Series, Asymmetrical)
+    ('double_Utube_series_asymmetrical', 0.2, None, 1., 2.1665182792127133),
+    ('double_Utube_series_asymmetrical', 0.2, None, np.array([1., 2., 3., 1.]), 2.694612953355227),
+    ('double_Utube_series_asymmetrical', 0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 2.9776742881734553),
+    ('double_Utube_series_asymmetrical', -0.2, None, 1., 2.1665183016853726),
+    ('double_Utube_series_asymmetrical', -0.2, None, np.array([1., 2., 3., 1.]), 2.6935999535110744),
+    ('double_Utube_series_asymmetrical', -0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 2.976372588013424),
     # Coaxial (Annular pipe is inlet pipe)
-    ('coaxial_annular_in', None, 1., 2.581130521333567),
-    ('coaxial_annular_in', None, np.array([1., 2., 3., 1.]), 3.0276625795763357),
-    ('coaxial_annular_in', np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 3.2740053393585953),
+    ('coaxial_annular_in', 0.2, None, 1., 2.581130521333567),
+    ('coaxial_annular_in', 0.2, None, np.array([1., 2., 3., 1.]), 3.0276625795763357),
+    ('coaxial_annular_in', 0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 3.2740053393585953),
+    ('coaxial_annular_in', -0.2, None, 1., 2.5811305213335674),
+    ('coaxial_annular_in', -0.2, None, np.array([1., 2., 3., 1.]), 2.981638747649938),
+    ('coaxial_annular_in', -0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 3.216388652724723),
     # Coaxial (Annular pipe is outlet pipe)
-    ('coaxial_annular_out', None, 1., 2.5811305213335674),
-    ('coaxial_annular_out', None, np.array([1., 2., 3., 1.]), 2.981638747649938),
-    ('coaxial_annular_out', np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 3.216388652724723),
+    ('coaxial_annular_out', 0.2, None, 1., 2.5811305213335674),
+    ('coaxial_annular_out', 0.2, None, np.array([1., 2., 3., 1.]), 2.981638747649938),
+    ('coaxial_annular_out', 0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 3.216388652724723),
+    ('coaxial_annular_out', -0.2, None, 1., 2.581130521333567),
+    ('coaxial_annular_out', -0.2, None, np.array([1., 2., 3., 1.]), 3.0276625795763357),
+    ('coaxial_annular_out', -0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 3.2740053393585953),
     # Coaxial (Annular pipe is inlet pipe, no grout)
-    ('coaxial_no_grout', None, 1., 1.157646033266169),
-    ('coaxial_no_grout', None, np.array([1., 2., 3., 1.]), 1.1622084785087514),
-    ('coaxial_no_grout', np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 1.2705296410278433),
+    ('coaxial_no_grout', 0.2, None, 1., 1.157646033266169),
+    ('coaxial_no_grout', 0.2, None, np.array([1., 2., 3., 1.]), 1.1622084785087514),
+    ('coaxial_no_grout', 0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 1.2705296410278433),
+    ('coaxial_no_grout', -0.2, None, 1., 1.1576460332661673),
+    ('coaxial_no_grout', -0.2, None, np.array([1., 2., 3., 1.]), 1.549177719025166),
+    ('coaxial_no_grout', -0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 1.7932400538024473),
     ])
 def test_outlet_temperature(
-        pipe_fixture, segment_ratios, T_b, expected, request):
+        pipe_fixture, m_flow_borehole, segment_ratios, T_b, expected, request):
     # Extract pipe from fixture
     pipe = request.getfixturevalue(pipe_fixture)
     # Fluid is propylene-glycol 20%
     fluid = gt.media.Fluid('MPG', 20.)
-    m_flow_borehole = 0.2   # Fluid mass flow rate [kg/s]
     T_f_in = 5.0            # Inlet fluid temperature [degC]
     # Outlet fluid temperature [degC]
     T_f_out = pipe.get_outlet_temperature(
@@ -249,38 +349,69 @@ def test_outlet_temperature(
 
 
 # Test get_inlet_temperature
-@pytest.mark.parametrize("pipe_fixture, segment_ratios, T_b, expected", [
+@pytest.mark.parametrize("pipe_fixture, m_flow_borehole, segment_ratios, T_b, expected", [
     # Single U-tube
-    ('single_Utube', None, 1., 7.595314034714041),
-    ('single_Utube', None, np.array([1., 2., 3., 1.]), 8.33912674339739),
-    ('single_Utube', np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 8.73842016624424),
+    ('single_Utube', 0.2, None, 1., 7.595314034714041),
+    ('single_Utube', 0.2, None, np.array([1., 2., 3., 1.]), 8.33912674339739),
+    ('single_Utube', 0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 8.73842016624424),
+    ('single_Utube', -0.2, None, 1., 7.595314034714041),
+    ('single_Utube', -0.2, None, np.array([1., 2., 3., 1.]), 8.33912674339739),
+    ('single_Utube', -0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 8.73842016624424),
     # Double U-tube (Parallel)
-    ('double_Utube_parallel', None, 1., 5.7977998086638305),
-    ('double_Utube_parallel', None, np.array([1., 2., 3., 1.]), 6.526064048901171),
-    ('double_Utube_parallel', np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 6.923635874226969),
+    ('double_Utube_parallel', 0.2, None, 1., 5.7977998086638305),
+    ('double_Utube_parallel', 0.2, None, np.array([1., 2., 3., 1.]), 6.526064048901171),
+    ('double_Utube_parallel', 0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 6.923635874226969),
+    ('double_Utube_parallel', -0.2, None, 1., 5.7977998086638305),
+    ('double_Utube_parallel', -0.2, None, np.array([1., 2., 3., 1.]), 6.526064048901171),
+    ('double_Utube_parallel', -0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 6.923635874226969),
+    # Double U-tube (Parallel, Asymmetrical)
+    ('double_Utube_parallel_asymmetrical', 0.2, None, 1., 6.332237785409975),
+    ('double_Utube_parallel_asymmetrical', 0.2, None, np.array([1., 2., 3., 1.]), 7.063604583735618),
+    ('double_Utube_parallel_asymmetrical', 0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 7.4617314345431),
+    ('double_Utube_parallel_asymmetrical', -0.2, None, 1., 6.332237824683503),
+    ('double_Utube_parallel_asymmetrical', -0.2, None, np.array([1., 2., 3., 1.]), 7.060884152083968),
+    ('double_Utube_parallel_asymmetrical', -0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 7.458292579361653),
     # Double U-tube (Series)
-    ('double_Utube_series', None, 1., 5.864420235466437),
-    ('double_Utube_series', None, np.array([1., 2., 3., 1.]), 6.608840446656091),
-    ('double_Utube_series', np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 7.008202698042327),
+    ('double_Utube_series', 0.2, None, 1., 5.864420235466437),
+    ('double_Utube_series', 0.2, None, np.array([1., 2., 3., 1.]), 6.608840446656091),
+    ('double_Utube_series', 0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 7.008202698042327),
+    ('double_Utube_series', -0.2, None, 1., 5.864420235466437),
+    ('double_Utube_series', -0.2, None, np.array([1., 2., 3., 1.]), 6.608840446656091),
+    ('double_Utube_series', -0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 7.008202698042327),
+    # Double U-tube (Series, Asymmetrical)
+    ('double_Utube_series_asymmetrical', 0.2, None, 1., 6.324765610974025),
+    ('double_Utube_series_asymmetrical', 0.2, None, np.array([1., 2., 3., 1.]), 7.070271989601069),
+    ('double_Utube_series_asymmetrical', 0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 7.469867064087138),
+    ('double_Utube_series_asymmetrical', -0.2, None, 1., 6.324765653205336),
+    ('double_Utube_series_asymmetrical', -0.2, None, np.array([1., 2., 3., 1.]), 7.068841963476496),
+    ('double_Utube_series_asymmetrical', -0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 7.46802948555655),
     # Coaxial (Annular pipe is inlet pipe)
-    ('coaxial_annular_in', None, 1., 7.237470090568812),
-    ('coaxial_annular_in', None, np.array([1., 2., 3., 1.]), 7.97588456424095),
-    ('coaxial_annular_in', np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 8.383252984827369),
+    ('coaxial_annular_in', 0.2, None, 1., 7.237470090568812),
+    ('coaxial_annular_in', 0.2, None, np.array([1., 2., 3., 1.]), 7.97588456424095),
+    ('coaxial_annular_in', 0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 8.383252984827369),
+    ('coaxial_annular_in', -0.2, None, 1., 7.237470090568813),
+    ('coaxial_annular_in', -0.2, None, np.array([1., 2., 3., 1.]), 7.899776560345228),
+    ('coaxial_annular_in', -0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 8.287974281876208),
     # Coaxial (Annular pipe is outlet pipe)
-    ('coaxial_annular_out', None, 1., 7.237470090568813),
-    ('coaxial_annular_out', None, np.array([1., 2., 3., 1.]), 7.899776560345228),
-    ('coaxial_annular_out', np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 8.287974281876208),
+    ('coaxial_annular_out', 0.2, None, 1., 7.237470090568813),
+    ('coaxial_annular_out', 0.2, None, np.array([1., 2., 3., 1.]), 7.899776560345228),
+    ('coaxial_annular_out', 0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 8.287974281876208),
+    ('coaxial_annular_out', -0.2, None, 1., 7.237470090568812),
+    ('coaxial_annular_out', -0.2, None, np.array([1., 2., 3., 1.]), 7.97588456424095),
+    ('coaxial_annular_out', -0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 8.383252984827369),
     # Coaxial (Annular pipe is inlet pipe, no grout)
-    ('coaxial_no_grout', None, 1., 4.926662185940351),
-    ('coaxial_no_grout', None, np.array([1., 2., 3., 1.]), 4.93141182148365),
-    ('coaxial_no_grout', np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 5.044177239200924),
+    ('coaxial_no_grout', 0.2, None, 1., 4.926662185940351),
+    ('coaxial_no_grout', 0.2, None, np.array([1., 2., 3., 1.]), 4.93141182148365),
+    ('coaxial_no_grout', 0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 5.044177239200924),
+    ('coaxial_no_grout', -0.2, None, 1., 4.9266621859403505),
+    ('coaxial_no_grout', -0.2, None, np.array([1., 2., 3., 1.]), 5.33425783084843),
+    ('coaxial_no_grout', -0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 5.588333677988297),
     ])
-def test_inlet_temperature(pipe_fixture, segment_ratios, T_b, expected, request):
+def test_inlet_temperature(pipe_fixture, m_flow_borehole, segment_ratios, T_b, expected, request):
     # Extract pipe from fixture
     pipe = request.getfixturevalue(pipe_fixture)
     # Fluid is propylene-glycol 20%
     fluid = gt.media.Fluid('MPG', 20.)
-    m_flow_borehole = 0.2   # Fluid mass flow rate [kg/s]
     Q_f = -3000.0           # Total heat transfer rate [W]
     # Inlet fluid temperature [degC]
     T_f_in = pipe.get_inlet_temperature(
@@ -289,39 +420,70 @@ def test_inlet_temperature(pipe_fixture, segment_ratios, T_b, expected, request)
 
 
 # Test get_borehole_heat_extraction_rate
-@pytest.mark.parametrize("pipe_fixture, segment_ratios, T_b, expected", [
+@pytest.mark.parametrize("pipe_fixture,  m_flow_borehole, segment_ratios, T_b, expected", [
     # Single U-tube
-    ('single_Utube', None, 1., -1819.4736348927008),
-    ('single_Utube', None, np.array([1., 2., 3., 1.]), np.array([-507.98022943, -330.29924271, -155.92399643, -486.93326314])),
-    ('single_Utube', np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), np.array([-212.48822505, -496.72428991, -284.6464931, -305.65175979])),
+    ('single_Utube', 0.2, None, 1., -1819.4736348927008),
+    ('single_Utube', 0.2, None, np.array([1., 2., 3., 1.]), np.array([-507.98022943, -330.29924271, -155.92399643, -486.93326314])),
+    ('single_Utube', 0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), np.array([-212.48822505, -496.72428991, -284.6464931, -305.65175979])),
+    ('single_Utube', -0.2, None, 1., -1819.4736348927008),
+    ('single_Utube', -0.2, None, np.array([1., 2., 3., 1.]), np.array([-507.98022943, -330.29924271, -155.92399643, -486.93326314])),
+    ('single_Utube', -0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), np.array([-212.48822505, -496.72428991, -284.6464931, -305.65175979])),
     # Double U-tube (Parallel)
-    ('double_Utube_parallel', None, 1., -2501.14645849),
-    ('double_Utube_parallel', None, np.array([1., 2., 3., 1.]), np.array([-796.48662356, -444.22614316, -108.02227066, -697.03753979])),
-    ('double_Utube_parallel', np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), np.array([-345.6258447, -724.00374747, -270.96795891, -456.57868535])),
+    ('double_Utube_parallel', 0.2, None, 1., -2501.14645849),
+    ('double_Utube_parallel', 0.2, None, np.array([1., 2., 3., 1.]), np.array([-796.48662356, -444.22614316, -108.02227066, -697.03753979])),
+    ('double_Utube_parallel', 0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), np.array([-345.6258447, -724.00374747, -270.96795891, -456.57868535])),
+    ('double_Utube_parallel', -0.2, None, 1., -2501.14645849),
+    ('double_Utube_parallel', -0.2, None, np.array([1., 2., 3., 1.]), np.array([-796.48662356, -444.22614316, -108.02227066, -697.03753979])),
+    ('double_Utube_parallel', -0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), np.array([-345.6258447, -724.00374747, -270.96795891, -456.57868535])),
+    # Double U-tube (Parallel, Asymmetrical)
+    ('double_Utube_parallel_asymmetrical', 0.2, None, 1., -2250.46227924),
+    ('double_Utube_parallel_asymmetrical', 0.2, None, np.array([1., 2., 3., 1.]), np.array([-703.79676211, -404.13655923, -120.76862091, -610.28198901])),
+    ('double_Utube_parallel_asymmetrical', 0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), np.array([-302.88310864, -647.8184961, -268.449683, -395.8402785])),
+    ('double_Utube_parallel_asymmetrical', -0.2, None, 1., -2250.46226266),
+    ('double_Utube_parallel_asymmetrical', -0.2, None, np.array([1., 2., 3., 1.]), np.array([-692.25229944, -403.05977335, -124.3681675, -620.83425671])),
+    ('double_Utube_parallel_asymmetrical', -0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), np.array([-297.55241658, -641.99757362, -275.22968865, -402.14665104])),
     # Double U-tube (Series)
-    ('double_Utube_series', None, 1., -2466.89213085),
-    ('double_Utube_series', None, np.array([1., 2., 3., 1.]), np.array([-745.16518357, -428.05472293, -114.8035859, -719.7675482])),
-    ('double_Utube_series', np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), np.array([-323.00538071, -688.58422816, -281.10403131, -468.80150159])),
+    ('double_Utube_series', 0.2, None, 1., -2466.89213085),
+    ('double_Utube_series', 0.2, None, np.array([1., 2., 3., 1.]), np.array([-745.16518357, -428.05472293, -114.8035859, -719.7675482])),
+    ('double_Utube_series', 0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), np.array([-323.00538071, -688.58422816, -281.10403131, -468.80150159])),
+    ('double_Utube_series', -0.2, None, 1., -2466.89213085),
+    ('double_Utube_series', -0.2, None, np.array([1., 2., 3., 1.]), np.array([-745.16518357, -428.05472293, -114.8035859, -719.7675482])),
+    ('double_Utube_series', -0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), np.array([-323.00538071, -688.58422816, -281.10403131, -468.80150159])),
+    # Double U-tube (Series, Asymmetrical)
+    ('double_Utube_series_asymmetrical', 0.2, None, 1., -2253.62032373),
+    ('double_Utube_series_asymmetrical', 0.2, None, np.array([1., 2., 3., 1.]), np.array([-665.58263158, -397.23145264, -133.56003739, -637.22412054])),
+    ('double_Utube_series_asymmetrical', 0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), np.array([-285.37371409, -625.46956592, -287.25059674, -410.37047011])),
+    ('double_Utube_series_asymmetrical', -0.2, None, 1., -2253.62030586),
+    ('double_Utube_series_asymmetrical', -0.2, None, np.array([1., 2., 3., 1.]), np.array([-659.22194084, -396.43635342, -135.56899714, -643.17664398])),
+    ('double_Utube_series_asymmetrical', -0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), np.array([-282.48526241, -622.0735276 , -291.01925835, -413.92161061])),
     # Coaxial (Annular pipe is inlet pipe)
-    ('coaxial_annular_in', None, 1., -1923.85692048),
-    ('coaxial_annular_in', None, np.array([1., 2., 3., 1.]), np.array([-757.51176437, -346.76503548, -48.92829119, -415.50088061])),
-    ('coaxial_annular_in', np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), np.array([-335.1152758, -618.6742642, -139.8377907, -279.14900216])),
+    ('coaxial_annular_in', 0.2, None, 1., -1923.85692048),
+    ('coaxial_annular_in', 0.2, None, np.array([1., 2., 3., 1.]), np.array([-757.51176437, -346.76503548, -48.92829119, -415.50088061])),
+    ('coaxial_annular_in', 0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), np.array([-335.1152758, -618.6742642, -139.8377907, -279.14900216])),
+    ('coaxial_annular_in', -0.2, None, 1., -1923.85692048),
+    ('coaxial_annular_in', -0.2, None, np.array([1., 2., 3., 1.]), np.array([-480.81667849, -324.83211948, -133.10520419, -666.55719699])),
+    ('coaxial_annular_in', -0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), np.array([-205.74675004, -483.36307319, -299.14468223, -430.34747541])),
     # Coaxial (Annular pipe is outlet pipe)
-    ('coaxial_annular_out', None, 1., -1923.85692048),
-    ('coaxial_annular_out', None, np.array([1., 2., 3., 1.]), np.array([-480.81667849, -324.83211948, -133.10520419, -666.55719699])),
-    ('coaxial_annular_out', np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), np.array([-205.74675004, -483.36307319, -299.14468223, -430.34747541])),
+    ('coaxial_annular_out', 0.2, None, 1., -1923.85692048),
+    ('coaxial_annular_out', 0.2, None, np.array([1., 2., 3., 1.]), np.array([-480.81667849, -324.83211948, -133.10520419, -666.55719699])),
+    ('coaxial_annular_out', 0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), np.array([-205.74675004, -483.36307319, -299.14468223, -430.34747541])),
+    ('coaxial_annular_out', -0.2, None, 1., -1923.85692048),
+    ('coaxial_annular_out', -0.2, None, np.array([1., 2., 3., 1.]), np.array([-757.51176437, -346.76503548, -48.92829119, -415.50088061])),
+    ('coaxial_annular_out', -0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), np.array([-335.1152758, -618.6742642, -139.8377907, -279.14900216])),
     # Coaxial (Annular pipe is inlet pipe, no grout)
-    ('coaxial_no_grout', None, 1., -3056.03065193),
-    ('coaxial_no_grout', None, np.array([1., 2., 3., 1.]), np.array([-1316.35767661, -1055.61151182, 1327.31501919, -2007.74772473])),
-    ('coaxial_no_grout', np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), np.array([-896.16349915, -1263.97717837, 1166.37595544, -1972.48352864])),
+    ('coaxial_no_grout', 0.2, None, 1., -3056.03065193),
+    ('coaxial_no_grout', 0.2, None, np.array([1., 2., 3., 1.]), np.array([-1316.35767661, -1055.61151182, 1327.31501919, -2007.74772473])),
+    ('coaxial_no_grout', 0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), np.array([-896.16349915, -1263.97717837, 1166.37595544, -1972.48352864])),
+    ('coaxial_no_grout', -0.2, None, 1., -3056.03065193),
+    ('coaxial_no_grout', -0.2, None, np.array([1., 2., 3., 1.]), np.array([-3100.21951788, 803.00604377, 1013.56111771, -1460.97209942])),
+    ('coaxial_no_grout', -0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), np.array([-2723.0188446, 448.90146775, 1141.41601057, -1417.80719535])),
     ])
 def test_borehole_heat_extraction_rate(
-        pipe_fixture, segment_ratios, T_b, expected, request):
+        pipe_fixture, m_flow_borehole, segment_ratios, T_b, expected, request):
     # Extract pipe from fixture
     pipe = request.getfixturevalue(pipe_fixture)
     # Fluid is propylene-glycol 20%
     fluid = gt.media.Fluid('MPG', 20.)
-    m_flow_borehole = 0.2   # Fluid mass flow rate [kg/s]
     T_f_in = 5.0            # Inlet fluid temperature [degC]
     # Borehole heat extraction rates [W]
     Q_b = pipe.get_borehole_heat_extraction_rate(
@@ -330,39 +492,70 @@ def test_borehole_heat_extraction_rate(
 
 
 # Test get_fluid_heat_extraction_rate
-@pytest.mark.parametrize("pipe_fixture, segment_ratios, T_b, expected", [
+@pytest.mark.parametrize("pipe_fixture, m_flow_borehole, segment_ratios, T_b, expected", [
     # Single U-tube
-    ('single_Utube', None, 1., -1819.4736348927008),
-    ('single_Utube', None, np.array([1., 2., 3., 1.]), -1481.1367317058312),
-    ('single_Utube', np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), -1299.5107678418537),
+    ('single_Utube', 0.2, None, 1., -1819.4736348927008),
+    ('single_Utube', 0.2, None, np.array([1., 2., 3., 1.]), -1481.1367317058312),
+    ('single_Utube', 0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), -1299.5107678418537),
+    ('single_Utube', -0.2, None, 1., -1819.4736348927008),
+    ('single_Utube', -0.2, None, np.array([1., 2., 3., 1.]), -1481.1367317058312),
+    ('single_Utube', -0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), -1299.5107678418537),
     # Double U-tube (Parallel)
-    ('double_Utube_parallel', None, 1., -2501.146458493431),
-    ('double_Utube_parallel', None, np.array([1., 2., 3., 1.]), -2045.7725771641726),
-    ('double_Utube_parallel', np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), -1797.176236436576),
+    ('double_Utube_parallel', 0.2, None, 1., -2501.146458493431),
+    ('double_Utube_parallel', 0.2, None, np.array([1., 2., 3., 1.]), -2045.7725771641726),
+    ('double_Utube_parallel', 0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), -1797.176236436576),
+    ('double_Utube_parallel', -0.2, None, 1., -2501.146458493431),
+    ('double_Utube_parallel', -0.2, None, np.array([1., 2., 3., 1.]), -2045.7725771641726),
+    ('double_Utube_parallel', -0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), -1797.176236436576),
+    # Double U-tube (Parallel, Asymmetrical)
+    ('double_Utube_parallel_asymmetrical', 0.2, None, 1., -2250.4622792393657),
+    ('double_Utube_parallel_asymmetrical', 0.2, None, np.array([1., 2., 3., 1.]), -1838.9839312593854),
+    ('double_Utube_parallel_asymmetrical', 0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), -1614.9915662357362),
+    ('double_Utube_parallel_asymmetrical', -0.2, None, 1., -2250.4622626640376),
+    ('double_Utube_parallel_asymmetrical', -0.2, None, np.array([1., 2., 3., 1.]), -1840.514497003164),
+    ('double_Utube_parallel_asymmetrical', -0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), -1616.9263298898907),
     # Double U-tube (Series)
-    ('double_Utube_series', None, 1., -2466.892130845958),
-    ('double_Utube_series', None, np.array([1., 2., 3., 1.]), -2007.7910405893476),
-    ('double_Utube_series', np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), -1761.49514176394),
+    ('double_Utube_series', 0.2, None, 1., -2466.892130845958),
+    ('double_Utube_series', 0.2, None, np.array([1., 2., 3., 1.]), -2007.7910405893476),
+    ('double_Utube_series', 0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), -1761.49514176394),
+    ('double_Utube_series', -0.2, None, 1., -2466.892130845958),
+    ('double_Utube_series', -0.2, None, np.array([1., 2., 3., 1.]), -2007.7910405893476),
+    ('double_Utube_series', -0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), -1761.49514176394),
+    # Double U-tube (Series, Asymmetrical)
+    ('double_Utube_series_asymmetrical', 0.2, None, 1., -2253.620323731943),
+    ('double_Utube_series_asymmetrical', 0.2, None, np.array([1., 2., 3., 1.]), -1833.5982421455162),
+    ('double_Utube_series_asymmetrical', 0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), -1608.4643468642694),
+    ('double_Utube_series_asymmetrical', -0.2, None, 1., -2253.6203058582287),
+    ('double_Utube_series_asymmetrical', -0.2, None, np.array([1., 2., 3., 1.]), -1834.4039353744404),
+    ('double_Utube_series_asymmetrical', -0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), -1609.4996589732305),
     # Coaxial (Annular pipe is inlet pipe)
-    ('coaxial_annular_in', None, 1., -1923.85692048),
-    ('coaxial_annular_in', None, np.array([1., 2., 3., 1.]), -1568.705971637178),
-    ('coaxial_annular_in', np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), -1372.776332855085),
+    ('coaxial_annular_in', 0.2, None, 1., -1923.85692048),
+    ('coaxial_annular_in', 0.2, None, np.array([1., 2., 3., 1.]), -1568.705971637178),
+    ('coaxial_annular_in', 0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), -1372.776332855085),
+    ('coaxial_annular_in', -0.2, None, 1., -1923.85692048),
+    ('coaxial_annular_in', -0.2, None, np.array([1., 2., 3., 1.]), -1605.3111991367698),
+    ('coaxial_annular_in', -0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), -1418.6019808667165),
     # Coaxial (Annular pipe is outlet pipe)
-    ('coaxial_annular_out', None, 1., -1923.85692048),
-    ('coaxial_annular_out', None, np.array([1., 2., 3., 1.]), -1605.3111991367698),
-    ('coaxial_annular_out', np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), -1418.6019808667165),
+    ('coaxial_annular_out', 0.2, None, 1., -1923.85692048),
+    ('coaxial_annular_out', 0.2, None, np.array([1., 2., 3., 1.]), -1605.3111991367698),
+    ('coaxial_annular_out', 0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), -1418.6019808667165),
+    ('coaxial_annular_out', -0.2, None, 1., -1923.85692048),
+    ('coaxial_annular_out', -0.2, None, np.array([1., 2., 3., 1.]), -1568.705971637178),
+    ('coaxial_annular_out', -0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), -1372.776332855085),
     # Coaxial (Annular pipe is inlet pipe, no grout)
-    ('coaxial_no_grout', None, 1., -3056.0306519279184),
-    ('coaxial_no_grout', None, np.array([1., 2., 3., 1.]), -3052.4018939764674),
-    ('coaxial_no_grout', np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), -2966.2482507211057),
+    ('coaxial_no_grout', 0.2, None, 1., -3056.0306519279184),
+    ('coaxial_no_grout', 0.2, None, np.array([1., 2., 3., 1.]), -3052.4018939764674),
+    ('coaxial_no_grout', 0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), -2966.2482507211057),
+    ('coaxial_no_grout', -0.2, None, 1., -3056.0306519279193),
+    ('coaxial_no_grout', -0.2, None, np.array([1., 2., 3., 1.]), -2744.6244558200647),
+    ('coaxial_no_grout', -0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), -2550.5085616265683),
     ])
 def test_fluid_heat_extraction_rate(
-        pipe_fixture, segment_ratios, T_b, expected, request):
+        pipe_fixture, m_flow_borehole, segment_ratios, T_b, expected, request):
     # Extract pipe from fixture
     pipe = request.getfixturevalue(pipe_fixture)
     # Fluid is propylene-glycol 20%
     fluid = gt.media.Fluid('MPG', 20.)
-    m_flow_borehole = 0.2   # Fluid mass flow rate [kg/s]
     T_f_in = 5.0            # Inlet fluid temperature [degC]
     # Fluid heat extraction rate [W]
     Q_f = pipe.get_fluid_heat_extraction_rate(
@@ -371,39 +564,70 @@ def test_fluid_heat_extraction_rate(
 
 
 # Test get_total_heat_extraction_rate
-@pytest.mark.parametrize("pipe_fixture, segment_ratios, T_b, expected", [
+@pytest.mark.parametrize("pipe_fixture, m_flow_borehole, segment_ratios, T_b, expected", [
     # Single U-tube
-    ('single_Utube', None, 1., -1819.4736348927008),
-    ('single_Utube', None, np.array([1., 2., 3., 1.]), -1481.1367317058312),
-    ('single_Utube', np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), -1299.5107678418537),
+    ('single_Utube', 0.2, None, 1., -1819.4736348927008),
+    ('single_Utube', 0.2, None, np.array([1., 2., 3., 1.]), -1481.1367317058312),
+    ('single_Utube', 0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), -1299.5107678418537),
+    ('single_Utube', -0.2, None, 1., -1819.4736348927008),
+    ('single_Utube', -0.2, None, np.array([1., 2., 3., 1.]), -1481.1367317058312),
+    ('single_Utube', -0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), -1299.5107678418537),
     # Double U-tube (Parallel)
-    ('double_Utube_parallel', None, 1., -2501.146458493431),
-    ('double_Utube_parallel', None, np.array([1., 2., 3., 1.]), -2045.7725771641726),
-    ('double_Utube_parallel', np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), -1797.176236436576),
+    ('double_Utube_parallel', 0.2, None, 1., -2501.146458493431),
+    ('double_Utube_parallel', 0.2, None, np.array([1., 2., 3., 1.]), -2045.7725771641726),
+    ('double_Utube_parallel', 0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), -1797.176236436576),
+    ('double_Utube_parallel', -0.2, None, 1., -2501.146458493431),
+    ('double_Utube_parallel', -0.2, None, np.array([1., 2., 3., 1.]), -2045.7725771641726),
+    ('double_Utube_parallel', -0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), -1797.176236436576),
+    # Double U-tube (Parallel, Asymmetrical)
+    ('double_Utube_parallel_asymmetrical', 0.2, None, 1., -2250.4622792393657),
+    ('double_Utube_parallel_asymmetrical', 0.2, None, np.array([1., 2., 3., 1.]), -1838.9839312593854),
+    ('double_Utube_parallel_asymmetrical', 0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), -1614.9915662357362),
+    ('double_Utube_parallel_asymmetrical', -0.2, None, 1., -2250.4622626640376),
+    ('double_Utube_parallel_asymmetrical', -0.2, None, np.array([1., 2., 3., 1.]), -1840.514497003164),
+    ('double_Utube_parallel_asymmetrical', -0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), -1616.9263298898907),
     # Double U-tube (Series)
-    ('double_Utube_series', None, 1., -2466.892130845958),
-    ('double_Utube_series', None, np.array([1., 2., 3., 1.]), -2007.7910405893476),
-    ('double_Utube_series', np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), -1761.49514176394),
+    ('double_Utube_series', 0.2, None, 1., -2466.892130845958),
+    ('double_Utube_series', 0.2, None, np.array([1., 2., 3., 1.]), -2007.7910405893476),
+    ('double_Utube_series', 0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), -1761.49514176394),
+    ('double_Utube_series', -0.2, None, 1., -2466.892130845958),
+    ('double_Utube_series', -0.2, None, np.array([1., 2., 3., 1.]), -2007.7910405893476),
+    ('double_Utube_series', -0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), -1761.49514176394),
+    # Double U-tube (Series, Asymmetrical)
+    ('double_Utube_series_asymmetrical', 0.2, None, 1., -2253.620323731943),
+    ('double_Utube_series_asymmetrical', 0.2, None, np.array([1., 2., 3., 1.]), -1833.5982421455162),
+    ('double_Utube_series_asymmetrical', 0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), -1608.4643468642694),
+    ('double_Utube_series_asymmetrical', -0.2, None, 1., -2253.6203058582287),
+    ('double_Utube_series_asymmetrical', -0.2, None, np.array([1., 2., 3., 1.]), -1834.4039353744404),
+    ('double_Utube_series_asymmetrical', -0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), -1609.4996589732305),
     # Coaxial (Annular pipe is inlet pipe)
-    ('coaxial_annular_in', None, 1., -1923.85692048),
-    ('coaxial_annular_in', None, np.array([1., 2., 3., 1.]), -1568.705971637178),
-    ('coaxial_annular_in', np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), -1372.776332855085),
+    ('coaxial_annular_in', 0.2, None, 1., -1923.85692048),
+    ('coaxial_annular_in', 0.2, None, np.array([1., 2., 3., 1.]), -1568.705971637178),
+    ('coaxial_annular_in', 0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), -1372.776332855085),
+    ('coaxial_annular_in', -0.2, None, 1., -1923.85692048),
+    ('coaxial_annular_in', -0.2, None, np.array([1., 2., 3., 1.]), -1605.3111991367698),
+    ('coaxial_annular_in', -0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), -1418.6019808667165),
     # Coaxial (Annular pipe is outlet pipe)
-    ('coaxial_annular_out', None, 1., -1923.85692048),
-    ('coaxial_annular_out', None, np.array([1., 2., 3., 1.]), -1605.3111991367698),
-    ('coaxial_annular_out', np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), -1418.6019808667165),
+    ('coaxial_annular_out', 0.2, None, 1., -1923.85692048),
+    ('coaxial_annular_out', 0.2, None, np.array([1., 2., 3., 1.]), -1605.3111991367698),
+    ('coaxial_annular_out', 0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), -1418.6019808667165),
+    ('coaxial_annular_out', -0.2, None, 1., -1923.85692048),
+    ('coaxial_annular_out', -0.2, None, np.array([1., 2., 3., 1.]), -1568.705971637178),
+    ('coaxial_annular_out', -0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), -1372.776332855085),
     # Coaxial (Annular pipe is inlet pipe, no grout)
-    ('coaxial_no_grout', None, 1., -3056.0306519279184),
-    ('coaxial_no_grout', None, np.array([1., 2., 3., 1.]), -3052.4018939764674),
-    ('coaxial_no_grout', np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), -2966.2482507211057),
+    ('coaxial_no_grout', 0.2, None, 1., -3056.0306519279184),
+    ('coaxial_no_grout', 0.2, None, np.array([1., 2., 3., 1.]), -3052.4018939764674),
+    ('coaxial_no_grout', 0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), -2966.2482507211057),
+    ('coaxial_no_grout', -0.2, None, 1., -3056.0306519279193),
+    ('coaxial_no_grout', -0.2, None, np.array([1., 2., 3., 1.]), -2744.6244558200647),
+    ('coaxial_no_grout', -0.2, np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), -2550.5085616265683),
     ])
 def test_total_heat_extraction_rate(
-        pipe_fixture, segment_ratios, T_b, expected, request):
+        pipe_fixture, m_flow_borehole, segment_ratios, T_b, expected, request):
     # Extract pipe from fixture
     pipe = request.getfixturevalue(pipe_fixture)
     # Fluid is propylene-glycol 20%
     fluid = gt.media.Fluid('MPG', 20.)
-    m_flow_borehole = 0.2   # Fluid mass flow rate [kg/s]
     T_f_in = 5.0            # Inlet fluid temperature [degC]
     # Total heat extraction rate [W]
     Q_t = pipe.get_total_heat_extraction_rate(