diff --git a/CHANGELOG.md b/CHANGELOG.md
index 993648b..d4238a7 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -17,6 +17,7 @@
 * [Issue 274](https://github.com/MassimoCimmino/pygfunction/issues/274) - Fixed scalar assignment from ndim-1 array. It is deprecated as of `numpy` version `1.25`. Only ndim-0 arrays can be treated as scalars.
 * [Issue 285](https://github.com/MassimoCimmino/pygfunction/issues/285) - Use `numpy.complex128` instead of `numpy.cfloat`. This is to comply with backward-incompatible changes introduced in `numpy` version `2.0`.
 * [Issue 286](https://github.com/MassimoCimmino/pygfunction/issues/286) - Fixed incorrect coefficients in `pipes.SingleUTube._continuity_condition_base` which caused errors in all dependent class methods when `segment_ratios` were not symmetric around the borehole mid-length.
+* [Issue 298](https://github.com/MassimoCimmino/pygfunction/issues/298) - Fixed incorrect coefficients in `pipes._basePipe`, `pipes.MultipleUTube` and `pipes.IndependentMultipleUTube` which caused errors in fluid temperature profiles and outlet fluid temperatures.
 
 ## Version 2.2.2 (2023-01-09)
 
diff --git a/pygfunction/pipes.py b/pygfunction/pipes.py
index bfc0115..f74f12f 100644
--- a/pygfunction/pipes.py
+++ b/pygfunction/pipes.py
@@ -376,34 +376,54 @@ def coefficients_outlet_temperature(
                 m_flow_borehole, cp_f, nSegments, segment_ratios, method_id):
             a_in, a_b = self._get_stored_coefficients(method_id)
         else:
-            # Check if _continuity_condition_base need to be called
-            # method_id for _continuity_condition_base is 0
+            # Check if _continuity_condition_head need to be called
+            # method_id for _continuity_condition_head is 1
             if self._check_coefficients(
-                    m_flow_borehole, cp_f, nSegments, segment_ratios, 0):
-                b_in, b_out, b_b = self._get_stored_coefficients(0)
+                    m_flow_borehole, cp_f, nSegments, segment_ratios, 1):
+                c_in, c_fu = self._get_stored_coefficients(1)
             else:
-                # Coefficient matrices from continuity condition:
-                # [b_out]*[T_{f,out}] = [b_in]*[T_{f,in}] + [b_b]*[T_b]
-                b_in, b_out, b_b = self._continuity_condition_base(
+                # Coefficient matrices for temperatures at depth (z = 0):
+                # [T_fd](z=0) = [c_in]*[T_{f,in}] + [c_fu]*[T_fu](z=0)
+                c_in, c_fu = self._continuity_condition_head(
                     m_flow_borehole, cp_f, nSegments,
                     segment_ratios=segment_ratios)
 
                 # Store coefficients
                 self._set_stored_coefficients(
                     m_flow_borehole, cp_f, nSegments, segment_ratios,
-                    (b_in, b_out, b_b), 0)
+                    (c_in, c_fu), 1)
+
+            # Coefficient matrix for connectivity at borehole outlet:
+            # [T_{f,out}] = [d_fu]*[T_fu](z=0)
+            d_fu = self._pipe_connectivity(
+                m_flow_borehole, cp_f, nSegments,
+                segment_ratios=segment_ratios)
+
+            # Coefficient matrices from general solution:
+            # [T_f](z=H) = [e_f0]*[T_f](0) + [e_b]*[T_b]
+            e_f0, e_b = self._general_solution(
+                self.b.H, 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]
-            b_out_m1 = np.linalg.inv(b_out)
-            a_in = b_out_m1 @ b_in
-            a_b = b_out_m1 @ b_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_b = A @ e_b
 
             # Store coefficients
             self._set_stored_coefficients(
                 m_flow_borehole, cp_f, nSegments, segment_ratios, (a_in, a_b),
                 method_id)
-
         return a_in, a_b
 
     def coefficients_temperature(
@@ -457,28 +477,44 @@ def coefficients_temperature(
         # method_id for _continuity_condition_head is 1
         if self._check_coefficients(
                 m_flow_borehole, cp_f, nSegments, segment_ratios, 1):
-            c_in, c_out, c_b = self._get_stored_coefficients(1)
+            c_in, c_fu = self._get_stored_coefficients(1)
         else:
             # Coefficient matrices for temperatures at depth (z = 0):
-            # [T_f](0) = [c_in]*[T_{f,in}] + [c_out]*[T_{f,out}]
-            #                              + [c_b]*[T_b]
-            c_in, c_out, c_b = self._continuity_condition_head(
-                m_flow_borehole, cp_f, nSegments, segment_ratios=segment_ratios)
+            # [T_fd](z=0) = [c_in]*[T_{f,in}] + [c_fu]*[T_fu](z=0)
+            c_in, c_fu = self._continuity_condition_head(
+                m_flow_borehole, cp_f, nSegments,
+                segment_ratios=segment_ratios)
 
             # Store coefficients
             self._set_stored_coefficients(
                 m_flow_borehole, cp_f, nSegments, segment_ratios,
-                (c_in, c_out, c_b), 1)
+                (c_in, c_fu), 1)
+
+        # Coefficient matrices from general solution:
+        # [T_f](z=H) = [e_f0]*[T_f](0) + [e_b]*[T_b]
+        e_f0, e_b = self._general_solution(
+            self.b.H, m_flow_borehole, cp_f, nSegments,
+            segment_ratios=segment_ratios)
 
         # Coefficient matrices from general solution:
-        # [T_f](z) = [d_f0]*[T_f](0) + [d_b]*[T_b]
-        d_f0, d_b = self._general_solution(
+        # [T_f](z) = [e_f0]*[T_f](0) + [e_b]*[T_b]
+        f_f0, f_b = self._general_solution(
             z, m_flow_borehole, cp_f, nSegments, segment_ratios=segment_ratios)
 
-        # Final coefficient matrices for temperatures at depth (z):
-        # [T_f](z) = [a_in]*[T_{f,in}] + [a_b]*[T_b]
-        a_in = d_f0 @ (c_in + c_out @ b_in)
-        a_b = d_f0 @ (c_b + c_out @ b_b) + d_b
+        # 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
+        a_b = (f_f0 @ A @ e_b) + f_b
 
         return a_in, a_b
 
@@ -891,26 +927,25 @@ def _check_geometry(self):
 
         return True
 
-    def _continuity_condition_base(
+    def _pipe_connectivity(
             self, m_flow_borehole, cp_f, nSegments, segment_ratios=None):
         """ Returns coefficients for the relation
-            [a_out]*[T_{f,out}] = [a_in]*[T_{f,in}] + [a_b]*[T_b]
+            [T_{f,out}] = [c_fu]*[T_fu](z=0)
         """
         raise NotImplementedError(
-            '_continuity_condition_base class method not implemented, '
+            '_continuity_condition_head class method not implemented, '
             'this method should return matrices for the relation: '
-            '[a_out]*[T_{f,out}] = [a_in]*[T_{f,in}] + [a_b]*[T_b]')
+            '[T_{f,out}] = [c_fu]*[T_fu](z=0)')
 
     def _continuity_condition_head(
             self, m_flow_borehole, cp_f, nSegments, segment_ratios=None):
         """ Returns coefficients for the relation
-            [T_f](z=0) = [a_in]*[T_{f,in}] + [a_out]*[T_{f,out}] + [a_b]*[T_b]
+            [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](z=0) = [a_in]*[T_{f,in}] + [a_out]*[T_{f,out}] '
-            '+ [a_b]*[T_b]')
+            '[T_fd](z=0) = [c_in]*[T_{f,in}] + [c_fu]*[T_fu](z=0)')
 
     def _general_solution(
             self, z, m_flow_borehole, cp_f, nSegments, segment_ratios=None):
@@ -928,7 +963,7 @@ def _update_model_variables(
         Evaluate common coefficients needed in other class methods.
         """
         raise NotImplementedError(
-            '_update_coefficients class method not implemented, '
+            '_update_model_variables class method not implemented, '
             'this method should evaluate common coefficients needed in other '
             'class methods.')
 
@@ -1283,7 +1318,7 @@ def update_thermal_resistances(self, R_fp):
         self._initialize_stored_coefficients()
         return
 
-    def _continuity_condition_base(
+    def _pipe_connectivity(
             self, m_flow_borehole, cp_f, nSegments, segment_ratios=None):
         """
         Equation that satisfies equal fluid temperatures in both legs of
@@ -1293,9 +1328,7 @@ def _continuity_condition_base(
 
             .. math::
 
-                \\mathbf{a_{out}} T_{f,out} =
-                \\mathbf{a_{in}} \\mathbf{T_{f,in}}
-                + \\mathbf{a_{b}} \\mathbf{T_b}
+                T_{f,out} = \\mathbf{c_{fu}} \\mathbf{T_{fu}}(z=0)
 
         Parameters
         ----------
@@ -1313,35 +1346,18 @@ def _continuity_condition_base(
 
         Returns
         -------
-        a_in : (nOutlets, nInlets,) array
-            Array of coefficients for inlet fluid temperature.
-        a_out : (nOutlets, nOutlets,) array
-            Array of coefficients for outlet fluid temperature.
-        a_b : (nOutlets, nSegments,) array
-            Array of coefficients for borehole wall temperatures.
+        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)
 
-        # Evaluate coefficient matrices from Hellstrom (1991):
-        a_in = ((self._f1(self.b.H) + self._f2(self.b.H))
-                / (self._f3(self.b.H) - self._f2(self.b.H)))
-        a_in = np.array([[a_in]])
+        # The upward fluid temperature at z=0 is the outlet fluid temperature
+        c_fu = np.array([[1.]])
 
-        a_out = np.array([[1.0]])
-
-        a_b = np.zeros((self.nOutlets, nSegments))
-
-        z = self.b._segment_edges(nSegments, segment_ratios=segment_ratios)
-        F4 = self._F4(self.b.H - z)
-        dF4 = F4[:-1] - F4[1:]
-        F5 = self._F5(self.b.H - z)
-        dF5 = F5[:-1] - F5[1:]
-        a_b[0, :] = (dF4 + dF5) / (self._f3(self.b.H) - self._f2(self.b.H))
-
-        return a_in, a_out, a_b
+        return c_fu
 
     def _continuity_condition_head(
             self, m_flow_borehole, cp_f, nSegments, segment_ratios=None):
@@ -1354,9 +1370,8 @@ def _continuity_condition_head(
 
             .. math::
 
-                \\mathbf{T_f}(z=0) = \\mathbf{a_{in}} \\mathbf{T_{f,in}}
-                + \\mathbf{a_{out}} \\mathbf{T_{f,out}}
-                + \\mathbf{a_{b}} \\mathbf{T_{b}}
+                \\mathbf{T_fd}(z=0) = \\mathbf{c_{in}} \\mathbf{T_{f,in}}
+                + \\mathbf{c_{fu}} \\mathbf{T_{fu}}(z=0)
 
         Parameters
         ----------
@@ -1374,24 +1389,22 @@ def _continuity_condition_head(
 
         Returns
         -------
-        a_in : (2*nPipes, nInlets,) array
+        c_in : (nPipes, nInlets,) array
             Array of coefficients for inlet fluid temperature.
-        a_out : (2*nPipes, nOutlets,) array
-            Array of coefficients for outlet fluid temperature.
-        a_b : (2*nPipes, nSegments,) array
-            Array of coefficients for borehole wall 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)
 
-        # There is only one pipe
-        a_in = np.array([[1.0], [0.0]])
-        a_out = np.array([[0.0], [1.0]])
-        a_b = np.zeros((2, nSegments))
+        # 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 a_in, a_out, a_b
+        return c_in, c_fu
 
     def _general_solution(
             self, z, m_flow_borehole, cp_f, nSegments, segment_ratios=None):
@@ -1647,9 +1660,9 @@ class MultipleUTube(_BasePipe):
     Contains information regarding the physical dimensions and thermal
     characteristics of the pipes and the grout material, as well as methods to
     evaluate fluid temperatures and heat extraction rates based on the work of
-    Cimmino [#Cimmino2016]_ for boreholes with any number of U-tubes. Internal
-    borehole thermal resistances are evaluated using the multipole method of
-    Claesson and Hellstrom [#Multiple-Claesson2011b]_.
+    Cimmino [#Cimmino2016]_, [#Cimmino2024]_ for boreholes with any number of
+    U-tubes. Internal borehole thermal resistances are evaluated using the
+    multipole method of Claesson and Hellstrom [#Multiple-Claesson2011b]_.
 
     Attributes
     ----------
@@ -1699,6 +1712,9 @@ 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.
+    .. [#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.
@@ -1751,7 +1767,7 @@ def update_thermal_resistances(self, R_fp):
         self._initialize_stored_coefficients()
         return
 
-    def _continuity_condition_base(
+    def _pipe_connectivity(
             self, m_flow_borehole, cp_f, nSegments, segment_ratios=None):
         """
         Equation that satisfies equal fluid temperatures in both legs of
@@ -1761,8 +1777,7 @@ def _continuity_condition_base(
 
             .. math::
 
-                \\mathbf{a_{out}} T_{f,out} = \\mathbf{a_{in}} T_{f,in}
-                + \\mathbf{a_{b}} \\mathbf{T_b}
+                T_{f,out} = \\mathbf{c_{fu}} \\mathbf{T_{fu}}(z=0)
 
         Parameters
         ----------
@@ -1780,67 +1795,28 @@ def _continuity_condition_base(
 
         Returns
         -------
-        a_in : (nOutlets, nInlets,) array
-            Array of coefficients for inlet fluid temperature.
-        a_out : (nOutlets, nOutlets,) array
-            Array of coefficients for outlet fluid temperature.
-        a_b : (nOutlets, nSegments,) array
-            Array of coefficients for borehole wall temperatures.
+        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)
-
-        # Coefficient matrices from continuity condition:
-        # [b_u]*[T_{f,u}](z=0) = [b_d]*[T_{f,d}](z=0) + [b_b]*[T_b]
-        b_d, b_u, b_b = self._continuity_condition(
-            m_flow_borehole, cp_f, nSegments, segment_ratios=segment_ratios)
-        b_u_m1 = np.linalg.inv(b_u)
+        m_flow_pipe = self._m_flow_pipe[-self.nPipes:]
+        cp_pipe = self._cp_pipe[-self.nPipes:]
 
         if self.config == 'parallel':
-            # Intermediate coefficient matrices:
-            # [T_{f,d}](z=0) = [c_in]*[T_{f,in}]
-            c_in = np.ones((self.nPipes, 1))
-
-            # Intermediate coefficient matrices:
-            # [T_{f,out}] = d_u*[T_{f,u}](z=0)
-            mcp = self._m_flow_pipe[-self.nPipes:]*self._cp_pipe[-self.nPipes:]
-            d_u = np.reshape(mcp/np.sum(mcp), (1, -1))
-
-            # Final coefficient matrices for continuity at depth (z = H):
-            # [a_out][T_{f,out}] = [a_in]*[T_{f,in}] + [a_b]*[T_b]
-            a_in = d_u @ b_u_m1 @ b_d @ c_in
-            a_out = np.array([[1.0]])
-            a_b = d_u @ b_u_m1 @ b_b
-
+            # 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':
-            # Intermediate coefficient matrices:
-            # [T_{f,d}](z=0) = [c_in]*[T_{f,in}] + [c_u]*[T_{f,u}](z=0)
-            c_in = np.eye(self.nPipes, M=1)
-            c_u = np.eye(self.nPipes, k=-1)
-
-            # Intermediate coefficient matrices:
-            # [d_u]*[T_{f,u}](z=0) = [d_in]*[T_{f,in}] + [d_b]*[T_b]
-            d_u = b_u - b_d @ c_u
-            d_in = b_d @ c_in
-            d_b = b_b
-            d_u_m1 = np.linalg.inv(d_u)
-
-            # Intermediate coefficient matrices:
-            # [T_{f,out}] = e_u*[T_{f,u}](z=0)
-            e_u = np.eye(self.nPipes, M=1, k=-self.nPipes+1).T
-
-            # Final coefficient matrices for continuity at depth (z = H):
-            # [a_out][T_{f,out}] = [a_in]*[T_{f,in}] + [a_b]*[T_b]
-            a_in = e_u @ d_u_m1 @ d_in
-            a_out = np.array([[1.0]])
-            a_b = e_u @ d_u_m1 @ d_b
+            # 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}' "
-                                      f"not implemented.")
+            raise NotImplementedError(
+                f"Configuration '{self.config}' not implemented.")
 
-        return a_in, a_out, a_b
+        return c_fu[np.newaxis, :]
 
     def _continuity_condition_head(
             self, m_flow_borehole, cp_f, nSegments, segment_ratios=None):
@@ -1853,9 +1829,8 @@ def _continuity_condition_head(
 
             .. math::
 
-                \\mathbf{T_f}(z=0) = \\mathbf{a_{in}} \\mathbf{T_{f,in}}
-                + \\mathbf{a_{out}} \\mathbf{T_{f,out}}
-                + \\mathbf{a_{b}} \\mathbf{T_{b}}
+                \\mathbf{T_fd}(z=0) = \\mathbf{c_{in}} \\mathbf{T_{f,in}}
+                + \\mathbf{c_{fu}} \\mathbf{T_{fu}}(z=0)
 
         Parameters
         ----------
@@ -1873,12 +1848,10 @@ def _continuity_condition_head(
 
         Returns
         -------
-        a_in : (2*nPipes, nInlets,) array
+        c_in : (nPipes, nInlets,) array
             Array of coefficients for inlet fluid temperature.
-        a_out : (2*nPipes, nOutlets,) array
-            Array of coefficients for outlet fluid temperature.
-        a_b : (2*nPipes, nSegments,) array
-            Array of coefficients for borehole wall temperature.
+        c_fu : (nPipes, nPipes,) array
+            Array of coefficients for upward fluid temperatures.
 
         """
         # Check if model variables need to be updated
@@ -1886,116 +1859,21 @@ def _continuity_condition_head(
             m_flow_borehole, cp_f, nSegments, segment_ratios)
 
         if self.config == 'parallel':
-            a_in = np.vstack((np.ones((self.nPipes, self.nInlets)),
-                              np.zeros((self.nPipes, self.nInlets))))
-            a_out = np.vstack((np.zeros((self.nPipes, self.nOutlets)),
-                               np.ones((self.nPipes, self.nOutlets))))
-            a_b = np.zeros((2*self.nPipes, nSegments))
-
+            # 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':
-            # Coefficient matrices from continuity condition:
-            # [b_u]*[T_{f,u}](z=0) = [b_d]*[T_{f,d}](z=0) + [b_b]*[T_b]
-            b_d, b_u, b_b = self._continuity_condition(
-                m_flow_borehole, cp_f, nSegments,
-                segment_ratios=segment_ratios)
-
-            # Intermediate coefficient matrices:
-            # [T_{f,d}](z=0) = [c_in]*[T_{f,in}] + [c_u]*[T_{f,u}](z=0)
-            c_in = np.eye(self.nPipes, M=1)
-            c_u = np.eye(self.nPipes, k=-1)
-
-            # Intermediate coefficient matrices:
-            # [d_u]*[T_{f,u}](z=0) = [d_in]*[T_{f,in}] + [d_b]*[T_b]
-            d_u = b_u - b_d @ c_u
-            d_in = b_d @ c_in
-            d_b = b_b
-            d_u_m1 = np.linalg.inv(d_u)
-
-            # Intermediate coefficient matrices:
-            # [T_f](z=0) = [e_d]*[T_{f,d}](z=0) + [e_u]*[T_{f,u}](z=0)
-            e_d = np.eye(2*self.nPipes, M=self.nPipes)
-            e_u = np.eye(2*self.nPipes, M=self.nPipes, k=-self.nPipes)
-
-            # Final coefficient matrices for temperatures at depth (z = 0):
-            # [T_f](z=0) = [a_in]*[T_{f,in}]+[a_out]*[T_{f,out}]+[a_b]*[T_b]
-            a_in = e_d @ (c_in + c_u @ d_u_m1 @ d_in) + e_u @ d_u_m1 @ d_in
-            a_out = np.zeros((2*self.nPipes, self.nOutlets))
-            a_b = e_d @ c_u @ d_u_m1 @ d_b + e_u @ d_u_m1 @ d_b
+            # 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 a_in, a_out, a_b
-
-    def _continuity_condition(
-            self, m_flow_borehole, cp_f, nSegments, segment_ratios=None):
-        """
-        Build coefficient matrices to evaluate fluid temperatures in downward
-        and upward flowing pipes at depth (z = 0).
+            raise NotImplementedError(
+                f"Configuration '{self.config}' not implemented.")
 
-        Returns coefficients for the relation:
-
-            .. math::
-
-                \\mathbf{a_{u}} \\mathbf{T_{f,u}}(z=0) =
-                + \\mathbf{a_{d}} \\mathbf{T_{f,d}}(z=0)
-                + \\mathbf{a_{b}} \\mathbf{T_{b}}
-
-        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
-        -------
-        a_d : (nPipes, nPipes,) array
-            Array of coefficients for fluid temperature in downward flowing
-            pipes.
-        a_u : (nPipes, nPipes,) array
-            Array of coefficients for fluid temperature in upward flowing
-            pipes.
-        a_b : (nPipes, nSegments,) array
-            Array of coefficients for borehole wall temperature.
-
-        """
-        # Load coefficients
-        sumA = self._sumA
-        V = self._V
-        Vm1 = self._Vm1
-        L = self._L
-        Dm1 = self._Dm1
-
-        # Matrix exponential at depth (z = H)
-        H = self.b.H
-        E = np.real(V @ np.diag(np.exp(L*H)) @ Vm1)
-
-        # Coefficient matrix for borehole wall temperatures
-        IIm1 = np.hstack((np.eye(self.nPipes), -np.eye(self.nPipes)))
-        a_b = np.zeros((self.nPipes, nSegments))
-        z = self.b._segment_edges(nSegments, segment_ratios=segment_ratios)[::-1]
-        exp_Lz = np.exp(np.multiply.outer(L, z))
-        dexp_Lz = exp_Lz[:,:-1] - exp_Lz[:,1:]
-        a_b = np.real(((IIm1 @ V @ Dm1) * (Vm1 @ sumA)) @ dexp_Lz)
-
-
-        # Configuration-specific inlet and outlet coefficient matrices
-        IZER = np.vstack((np.eye(self.nPipes),
-                          np.zeros((self.nPipes, self.nPipes))))
-        ZERI = np.vstack((np.zeros((self.nPipes, self.nPipes)),
-                          np.eye(self.nPipes)))
-        a_u = IIm1 @ E @ ZERI
-        a_d = -IIm1 @ E @ IZER
-
-        return a_d, a_u, a_b
+        return c_in[:, np.newaxis], c_fu
 
     def _general_solution(
             self, z, m_flow_borehole, cp_f, nSegments, segment_ratios=None):
@@ -2057,7 +1935,7 @@ def _general_solution(
             nSegments, segment_ratios=segment_ratios)
         dz = np.maximum(np.subtract.outer(z_array, z_edges), 0.)
         exp_Lz = np.exp(np.multiply.outer(L, dz)).transpose(1, 0, 2)
-        dexp_Lz = exp_Lz[:,:,1:] - exp_Lz[:,:,:-1]
+        dexp_Lz = exp_Lz[:, :, 1:] - exp_Lz[:, :, :-1]
         a_b = np.real(((V @ Dm1) * (Vm1 @ sumA)) @ dexp_Lz)
         # Remove first dimension if z is a scalar
         if np.isscalar(z):
@@ -2154,9 +2032,9 @@ class IndependentMultipleUTube(MultipleUTube):
     Contains information regarding the physical dimensions and thermal
     characteristics of the pipes and the grout material, as well as methods to
     evaluate fluid temperatures and heat extraction rates based on the work of
-    Cimmino [#Cimmino2016b]_ for boreholes with any number of U-tubes. Internal
-    borehole thermal resistances are evaluated using the multipole method of
-    Claesson and Hellstrom [#Independent-Claesson2011b]_.
+    Cimmino [#Cimmino2016b]_, [#Cimmino2024b]_ for boreholes with any number of
+    U-tubes. Internal borehole thermal resistances are evaluated using the
+    multipole method of Claesson and Hellstrom [#Independent-Claesson2011b]_.
 
     Attributes
     ----------
@@ -2199,11 +2077,13 @@ 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
+       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):
@@ -2247,7 +2127,7 @@ def update_thermal_resistances(self, R_fp):
         self._initialize_stored_coefficients()
         return
 
-    def _continuity_condition_base(
+    def _pipe_connectivity(
             self, m_flow_borehole, cp_f, nSegments, segment_ratios=None):
         """
         Equation that satisfies equal fluid temperatures in both legs of
@@ -2257,8 +2137,7 @@ def _continuity_condition_base(
 
             .. math::
 
-                \\mathbf{a_{out}} T_{f,out} = \\mathbf{a_{in}} T_{f,in}
-                + \\mathbf{a_{b}} \\mathbf{T_b}
+                T_{f,out} = \\mathbf{c_{fu}} \\mathbf{T_{fu}}(z=0)
 
         Parameters
         ----------
@@ -2276,27 +2155,21 @@ def _continuity_condition_base(
 
         Returns
         -------
-        a_in : (nOutlets, nInlets,) array
-            Array of coefficients for inlet fluid temperature.
-        a_out : (nOutlets, nOutlets,) array
-            Array of coefficients for outlet fluid temperature.
-        a_b : (nOutlets, nSegments,) array
-            Array of coefficients for borehole wall temperatures.
+        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=segment_ratios)
+            m_flow_borehole, cp_f, nSegments, segment_ratios)
 
-        # Coefficient matrices from continuity condition:
-        # [b_u]*[T_{f,u}](z=0) = [b_d]*[T_{f,d}](z=0) + [b_b]*[T_b]
-        a_in, a_out, a_b = self._continuity_condition(
-            m_flow_borehole, cp_f, nSegments, segment_ratios=segment_ratios)
+        # All upward flowing pipes are connceted to their respective outlet
+        c_fu = np.eye(self.nPipes)
 
-        return a_in, a_out, a_b
+        return c_fu
 
     def _continuity_condition_head(
-            self, m_flow_borehole, cp, nSegments, segment_ratios=None):
+            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
@@ -2306,9 +2179,8 @@ def _continuity_condition_head(
 
             .. math::
 
-                \\mathbf{T_f}(z=0) = \\mathbf{a_{in}} \\mathbf{T_{f,in}}
-                + \\mathbf{a_{out}} \\mathbf{T_{f,out}}
-                + \\mathbf{a_{b}} \\mathbf{T_{b}}
+                \\mathbf{T_fd}(z=0) = \\mathbf{c_{in}} \\mathbf{T_{f,in}}
+                + \\mathbf{c_{fu}} \\mathbf{T_{fu}}(z=0)
 
         Parameters
         ----------
@@ -2326,23 +2198,22 @@ def _continuity_condition_head(
 
         Returns
         -------
-        a_in : (2*nPipes, nInlets,) array
+        c_in : (nPipes, nInlets,) array
             Array of coefficients for inlet fluid temperature.
-        a_out : (2*nPipes, nOutlets,) array
-            Array of coefficients for outlet fluid temperature.
-        a_b : (2*nPipes, nSegments,) array
-            Array of coefficients for borehole wall 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, nSegments, segment_ratios)
+            m_flow_borehole, cp_f, nSegments, segment_ratios)
 
-        a_in = np.eye(2*self.nPipes, M=self.nPipes, k=0)
-        a_out = np.eye(2*self.nPipes, M=self.nPipes, k=-self.nPipes)
-        a_b = np.zeros((2*self.nPipes, nSegments))
+        # All downward flowing pipes are connected to their respective inlet
+        c_in = np.eye(self.nPipes)
+        # None of the upward flowing pipes are connected to another pipe
+        c_fu = np.zeros((self.nPipes, self.nPipes))
 
-        return a_in, a_out, a_b
+        return c_in, c_fu
 
     def _format_inputs(self, m_flow_borehole, cp_f, nSegments, segment_ratios):
         """
diff --git a/tests/pipes_test.py b/tests/pipes_test.py
index f209626..67328b1 100644
--- a/tests/pipes_test.py
+++ b/tests/pipes_test.py
@@ -153,44 +153,44 @@ def test_multipole(J, expected):
     ('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.]), np.array([[4.34676755, 3.07354134], [4.25566624, 3.13435325]])),
-    ('single_Utube', None, np.array([1., 2., 3., 1.]), np.array([65., 75.]), np.array([[4.41754093, 3.49949295], [4.35173147, 3.54346564]])),
-    ('single_Utube', np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), np.array([65., 75.]), np.array([[4.47282344, 3.70132653], [4.4261555, 3.72955733]])),
+    ('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]])),
     # 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.09862044, 4.09862044, 3.07439258, 3.07439258])),
-    ('double_Utube_parallel', None, 1., np.array([65., 75.]), np.array([[3.87525104, 3.87525104, 2.20313908, 2.20313908], [3.73265141, 3.73265141, 2.26719823, 2.26719823]])),
-    ('double_Utube_parallel', None, np.array([1., 2., 3., 1.]), np.array([65., 75.]), np.array([[4.00464852, 4.00464852, 2.84788608, 2.84788608], [3.90522192, 3.90522192, 2.89301847, 2.89301847]])),
-    ('double_Utube_parallel', np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), np.array([65., 75.]), np.array([[4.09862044, 4.09862044, 3.07439258, 3.07439258], [4.03210009, 4.03210009, 3.09222735, 3.09222735]])),
+    ('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 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.49394782, 3.15796672, 3.7037625, 2.90601555])),
-    ('double_Utube_series', None, 1., np.array([65., 75.]), np.array([[4.36908096, 2.53231146, 3.13441957, 2.03763963], [4.28094228, 2.49706752, 3.19348974, 2.0612353]])),
-    ('double_Utube_series', None, np.array([1., 2., 3., 1.]), np.array([65., 75.]), np.array([[4.44022419, 2.94528677, 3.54323578, 2.65057213], [4.37608674, 2.92420008, 3.58625745, 2.66472128]])),
-    ('double_Utube_series', np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), np.array([65., 75.]), np.array([[4.49394782, 3.15796672, 3.7037625, 2.90601555], [4.44772426, 3.15021634, 3.73124558, 2.90769612]])),
+    ('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]])),
     # 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.]), np.array([[3.15203088, 2.18408362], [2.96401382, 2.15051705]])),
-    ('coaxial_annular_in', None, np.array([1., 2., 3., 1.]), np.array([65., 75.]), np.array([[3.4176666, 2.73205968], [3.2920645, 2.7081367]])),
-    ('coaxial_annular_in', np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), np.array([65., 75.]), np.array([[3.57440691, 3.0183749], [3.49639377, 2.99877788]])),
+    ('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 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.]), np.array([[4.50649998, 2.92933532], [4.44976224, 3.02086677]])),
-    ('coaxial_annular_out', None, np.array([1., 2., 3., 1.]), np.array([65., 75.]), np.array([[4.62416027, 3.50307539], [4.58402116, 3.57860389]])),
-    ('coaxial_annular_out', np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), np.array([65., 75.]), np.array([[4.67083279, 3.72110166], [4.63666726, 3.76190344]])),
+    ('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 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.]), np.array([[3.77463898, 1.10936191], [2.8355954, 1.48415788]])),
-    ('coaxial_no_grout', None, np.array([1., 2., 3., 1.]), np.array([65., 140.]), np.array([[3.95937381, 2.30543206], [3.2799902, 1.60137175]])),
-    ('coaxial_no_grout', np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), np.array([65., 140.]), np.array([[4.07193759, 2.7536371], [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]])),
     ])
 def test_temperature(
         pipe_fixture, segment_ratios, T_b, z, expected, request):
@@ -216,11 +216,11 @@ def test_temperature(
     # 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.691532668379643),
+    ('double_Utube_parallel', np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 2.7404064201268588),
     # Double U-tube (Series)
-    ('double_Utube_series', None, 1., 1.8983711735742064),
-    ('double_Utube_series', None, np.array([1., 2., 3., 1.]), 2.4755999700741573),
-    ('double_Utube_series', np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 2.744018786582172),
+    ('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),
     # 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),
@@ -257,11 +257,11 @@ def test_outlet_temperature(
     # 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.861469307697622),
+    ('double_Utube_parallel', np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 6.923635874226969),
     # Double U-tube (Series)
-    ('double_Utube_series', None, 1., 5.8644202354664365),
-    ('double_Utube_series', None, np.array([1., 2., 3., 1.]), 6.60884044665609),
-    ('double_Utube_series', np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 6.955005421937585),
+    ('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),
     # 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),
@@ -297,11 +297,11 @@ def test_inlet_temperature(pipe_fixture, segment_ratios, T_b, expected, request)
     # 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([-342.52693756, -711.1206948, -251.40681559, -447.57582464])),
+    ('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 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([-319.92216781, -677.35565178, -267.23882139, -463.25689612])),
+    ('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])),
     # 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])),
@@ -336,13 +336,13 @@ def test_borehole_heat_extraction_rate(
     ('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),
     # Double U-tube (Parallel)
-    ('double_Utube_parallel', None, 1., -2501.14645849),
+    ('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.]), -1836.0481583644691),
+    ('double_Utube_parallel', np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), -1797.176236436576),
     # Double U-tube (Series)
-    ('double_Utube_series', None, 1., -2466.89213085),
-    ('double_Utube_series', None, np.array([1., 2., 3., 1.]), -2007.7910405893485),
-    ('double_Utube_series', np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), -1794.303127215247),
+    ('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),
     # 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),
@@ -377,13 +377,13 @@ def test_fluid_heat_extraction_rate(
     ('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),
     # Double U-tube (Parallel)
-    ('double_Utube_parallel', None, 1., -2501.14645849),
+    ('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.]), -1836.0481583644691),
+    ('double_Utube_parallel', np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), -1797.176236436576),
     # Double U-tube (Series)
-    ('double_Utube_series', None, 1., -2466.89213085),
-    ('double_Utube_series', None, np.array([1., 2., 3., 1.]), -2007.7910405893485),
-    ('double_Utube_series', np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), -1794.303127215247),
+    ('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),
     # 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),
@@ -418,19 +418,19 @@ def test_total_heat_extraction_rate(
 @pytest.mark.parametrize(
     "pipe_fixture, m_flow, T_f_in, segment_ratios, T_b, z, expected", [
     # Double U-tube
-    ('double_Utube_independent', np.array([0.2, 0.15]), np.array([5., -1.]), None, 1., 65., np.array([4.33561246, -0.53401739,  3.03985865, 0.28974217])),
-    ('double_Utube_independent', np.array([0.2, 0.15]), np.array([5., -1.]), None, np.array([1., 2., 3., 1.]), 65., np.array([4.40351925, -0.44632268,  3.43990994, 0.77984857])),
-    ('double_Utube_independent', np.array([0.2, 0.15]), np.array([5., -1.]), np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 65., np.array([4.45571533, -0.3797212, 3.59641415, 0.96525637])),
-    ('double_Utube_independent', np.array([0.2, 0.15]), np.array([5., -1.]), None, 1., np.array([65., 75.]), np.array([[4.33561246, -0.53401739, 3.03985865, 0.28974217], [4.2430049, -0.47133142, 3.10196228, 0.25289182]])),
-    ('double_Utube_independent', np.array([0.2, 0.15]), np.array([5., -1.]), None, np.array([1., 2., 3., 1.]), np.array([65., 75.]), np.array([[4.40351925, -0.44632268, 3.43990994, 0.77984857], [4.33450267, -0.35334337, 3.48624558, 0.72533714]])),
-    ('double_Utube_independent', np.array([0.2, 0.15]), np.array([5., -1.]), np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), np.array([65., 75.]), np.array([[4.45571533, -0.3797212, 3.59641415, 0.96525637], [4.40441296, -0.26381686, 3.62735271, 0.89129191]])),
+    ('double_Utube_independent', np.array([0.2, 0.15]), np.array([5., -1.]), None, 1., 65., np.array([4.33561246, -0.53401739, 3.03985865, 0.28974217])),
+    ('double_Utube_independent', np.array([0.2, 0.15]), np.array([5., -1.]), None, np.array([1., 2., 3., 1.]), 65., np.array([4.40351925, -0.44632268, 3.43990994, 0.77984857])),
+    ('double_Utube_independent', np.array([0.2, 0.15]), np.array([5., -1.]), np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 65., np.array([4.45572396, -0.37988472, 3.63010142, 1.00960735])),
+    ('double_Utube_independent', np.array([0.2, 0.15]), np.array([5., -1.]), None, 1., np.array([65., 75., 150.]), np.array([[4.33561246, -0.53401739, 3.03985865, 0.28974217], [4.2430049, -0.47133142, 3.10196228, 0.25289182], [3.61910357, -0.06417533, 3.61910357, -0.06417533]])),
+    ('double_Utube_independent', np.array([0.2, 0.15]), np.array([5., -1.]), None, np.array([1., 2., 3., 1.]), np.array([65., 75., 150.]), np.array([[4.40351925, -0.44632268, 3.43990994, 0.77984857], [4.33450267, -0.35334337, 3.48624558, 0.72533714], [3.86172238, 0.2372006, 3.86172238, 0.2372006]])),
+    ('double_Utube_independent', np.array([0.2, 0.15]), np.array([5., -1.]), np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), np.array([65., 75., 150.]), np.array([[4.45572396, -0.37988472, 3.63010142, 1.00960735], [4.40442447, -0.26400763, 3.66189091, 0.93718219], [3.99089068, 0.39795939, 3.99089068, 0.39795939]])),
     # Triple U-tube
     ('triple_Utube_independent', np.array([0.2, 0.15, 0.10]), np.array([5., -1., 7.]), None, 1., 65., np.array([4.42233734, -0.45836746, 5.35867251, 3.21820626, 0.54732768, 2.80410243])),
     ('triple_Utube_independent', np.array([0.2, 0.15, 0.10]), np.array([5., -1., 7.]), None, np.array([1., 2., 3., 1.]), 65., np.array([4.47430559, -0.3726503, 5.47798812, 3.56864252, 0.98391184, 3.36084568])),
-    ('triple_Utube_independent', np.array([0.2, 0.15, 0.10]), np.array([5., -1., 7.]), np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 65., np.array([4.51735491, -0.31260448, 5.56070713, 3.7101146, 1.15077001, 3.55786056])),
-    ('triple_Utube_independent', np.array([0.2, 0.15, 0.10]), np.array([5., -1., 7.]), None, 1., np.array([65., 75.]), np.array([[4.42233734, -0.45836746, 5.35867251, 3.21820626, 0.54732768, 2.80410243], [4.33960106, -0.38478969, 5.148416, 3.27807117, 0.49886893, 2.90891002]])),
-    ('triple_Utube_independent', np.array([0.2, 0.15, 0.10]), np.array([5., -1., 7.]), None, np.array([1., 2., 3., 1.]), np.array([65., 75.]), np.array([[4.47430559, -0.3726503, 5.47798812, 3.56864252, 0.98391184, 3.36084568], [4.41115886, -0.27214867, 5.30491836, 3.61292754, 0.92062819, 3.45008461]])),
-    ('triple_Utube_independent', np.array([0.2, 0.15, 0.10]), np.array([5., -1., 7.]), np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), np.array([65., 75.]), np.array([[4.51735491, -0.31260448, 5.56070713, 3.7101146, 1.15077001, 3.55786056], [4.4690976, -0.1924684, 5.41547435, 3.74095158, 1.07113754, 3.62431085]])),
+    ('triple_Utube_independent', np.array([0.2, 0.15, 0.10]), np.array([5., -1., 7.]), np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 65., np.array([4.51695038, -0.3115071, 5.56146368, 3.73728694, 1.18959709, 3.6163127])),
+    ('triple_Utube_independent', np.array([0.2, 0.15, 0.10]), np.array([5., -1., 7.]), None, 1., np.array([65., 75., 150.]), np.array([[4.42233734, -0.45836746, 5.35867251, 3.21820626, 0.54732768, 2.80410243], [4.33960106, -0.38478969, 5.148416, 3.27807117, 0.49886893, 2.90891002], [3.76853744, 0.10098669,  3.84628411, 3.76853744, 0.10098669, 3.84628411]])),
+    ('triple_Utube_independent', np.array([0.2, 0.15, 0.10]), np.array([5., -1., 7.]), None, np.array([1., 2., 3., 1.]), np.array([65., 75., 150.]), np.array([[4.47430559, -0.3726503, 5.47798812, 3.56864252, 0.98391184, 3.36084568], [4.41115886, -0.27214867, 5.30491836, 3.61292754, 0.92062819, 3.45008461], [3.9754362, 0.37402984, 4.20147852, 3.9754362, 0.37402984, 4.20147852]])),
+    ('triple_Utube_independent', np.array([0.2, 0.15, 0.10]), np.array([5., -1., 7.]), np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), np.array([65., 75., 150.]), np.array([[4.51695038, -0.3115071, 5.56146368, 3.73728694, 1.18959709, 3.6163127], [4.46863202, -0.19118263, 5.41633523, 3.76847629, 1.11127473, 3.68580011], [4.08533665, 0.51957984, 4.39206115, 4.08533665, 0.51957984, 4.39206115]])),
     ])
 def test_temperature_independent(
         pipe_fixture, m_flow, T_f_in, segment_ratios, T_b, z, expected,
@@ -451,11 +451,11 @@ def test_temperature_independent(
     # Double U-tube
     ('double_Utube_independent', np.array([0.2, 0.15]), np.array([5., -1.]), None, 1., np.array([2.66975268, 0.50433911])),
     ('double_Utube_independent', np.array([0.2, 0.15]), np.array([5., -1.]), None, np.array([1., 2., 3., 1.]), np.array([3.07152772, 0.9760115])),
-    ('double_Utube_independent', np.array([0.2, 0.15]), np.array([5., -1.]), np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), np.array([3.25838257, 1.19473512])),
+    ('double_Utube_independent', np.array([0.2, 0.15]), np.array([5., -1.]), np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), np.array([3.28701848, 1.23027312])),
     # Triple U-tube
     ('triple_Utube_independent', np.array([0.2, 0.15, 0.10]), np.array([5., -1., 7.]), None, 1., np.array([2.85759687, 0.84655363, 2.20811155])),
     ('triple_Utube_independent', np.array([0.2, 0.15, 0.10]), np.array([5., -1., 7.]), None, np.array([1., 2., 3., 1.]), np.array([3.22566425, 1.26691386, 2.70350023])),
-    ('triple_Utube_independent', np.array([0.2, 0.15, 0.10]), np.array([5., -1., 7.]), np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), np.array([3.39731878, 1.46195771, 2.93172173])),
+    ('triple_Utube_independent', np.array([0.2, 0.15, 0.10]), np.array([5., -1., 7.]), np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), np.array([3.42207691, 1.49331263, 2.97388805])),
     ])
 def test_outlet_temperature_independent(
         pipe_fixture, m_flow, T_f_in, segment_ratios, T_b, expected, request):
@@ -475,11 +475,11 @@ def test_outlet_temperature_independent(
     # Double U-tube
     ('double_Utube_independent', np.array([0.2, 0.15]), np.array([-3000., 2000.]), None, 1., np.array([7.40595748, -3.59946781])),
     ('double_Utube_independent', np.array([0.2, 0.15]), np.array([-3000., 2000.]), None, np.array([1., 2., 3., 1.]), np.array([8.15037424, -2.85931237])),
-    ('double_Utube_independent', np.array([0.2, 0.15]), np.array([-3000., 2000.]), np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), np.array([8.49653478, -2.51603365])),
+    ('double_Utube_independent', np.array([0.2, 0.15]), np.array([-3000., 2000.]), np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), np.array([8.54973576, -2.4604302])),
     # Triple U-tube
     ('triple_Utube_independent', np.array([0.2, 0.15, 0.10]), np.array([-3000., 2000., -2500.]), None, 1., np.array([7.87321014, -2.88443189, 8.87646527])),
     ('triple_Utube_independent', np.array([0.2, 0.15, 0.10]), np.array([-3000., 2000., -2500.]), None, np.array([1., 2., 3., 1.]), np.array([8.62102769, -2.14384786, 9.60745667])),
-    ('triple_Utube_independent', np.array([0.2, 0.15, 0.10]), np.array([-3000., 2000., -2500.]), np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), np.array([8.96946474, -1.80024805, 9.94472697])),
+    ('triple_Utube_independent', np.array([0.2, 0.15, 0.10]), np.array([-3000., 2000., -2500.]), np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), np.array([9.02076799, -1.74491241, 10.00533859])),
     ])
 def test_inlet_temperature_independent(
         pipe_fixture, m_flow, Q_f, segment_ratios, T_b, expected, request):
@@ -499,11 +499,11 @@ def test_inlet_temperature_independent(
     # Double U-tube
     ('double_Utube_independent', np.array([0.2, 0.15]), np.array([5., -1.]), None, 1., -956.00963184),
     ('double_Utube_independent', np.array([0.2, 0.15]), np.array([5., -1.]), None, np.array([1., 2., 3., 1.]), np.array([-311.65264343, -12.22888682, 289.43868666, -320.65369703])),
-    ('double_Utube_independent', np.array([0.2, 0.15]), np.array([5., -1.]), np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), np.array([-136.63705522, -61.48218337, 405.20011405, -214.40320201])),
+    ('double_Utube_independent', np.array([0.2, 0.15]), np.array([5., -1.]), np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), np.array([-138.6403424, -69.24200957, 394.70262906, -218.85448613])),
     # Triple U-tube
     ('triple_Utube_independent', np.array([0.2, 0.15, 0.10]), np.array([5., -1., 7.]), None, 1., -2508.09408199),
     ('triple_Utube_independent', np.array([0.2, 0.15, 0.10]), np.array([5., -1., 7.]), None, np.array([1., 2., 3., 1.]), np.array([-754.8446691, -351.75861521, 43.23971775, -704.23081162])),
-    ('triple_Utube_independent', np.array([0.2, 0.15, 0.10]), np.array([5., -1., 7.]), np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), np.array([-321.07744291, -563.32294806, -1.60256316, -449.20675261])),
+    ('triple_Utube_independent', np.array([0.2, 0.15, 0.10]), np.array([5., -1., 7.]), np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), np.array([-323.63535083, -573.51851525, -16.07605399, -455.56959127])),
     ])
 def test_borehole_heat_extraction_rate_independent(
         pipe_fixture, m_flow, T_f_in, segment_ratios, T_b, expected, request):
@@ -523,11 +523,11 @@ def test_borehole_heat_extraction_rate_independent(
     # Double U-tube
     ('double_Utube_independent', np.array([0.2, 0.15]), np.array([5., -1.]), None, 1., np.array([-1853.37094997, 897.36131814])),
     ('double_Utube_independent', np.array([0.2, 0.15]), np.array([5., -1.]), None, np.array([1., 2., 3., 1.]), np.array([-1533.81766537, 1178.72112475])),
-    ('double_Utube_independent', np.array([0.2, 0.15]), np.array([5., -1.]), np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), np.array([-1385.20195806, 1309.19311228])),
+    ('double_Utube_independent', np.array([0.2, 0.15]), np.array([5., -1.]), np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), np.array([-1362.42628578, 1330.39207674])),
     # Triple U-tube
     ('triple_Utube_independent', np.array([0.2, 0.15, 0.10]), np.array([5., -1., 7.]), None, 1., np.array([-1703.96836927, 1101.4975216, -1905.62323433])),
     ('triple_Utube_independent', np.array([0.2, 0.15, 0.10]), np.array([5., -1., 7.]), None, np.array([1., 2., 3., 1.]), np.array([-1411.22460372, 1352.24883733, -1708.61861179])),
-    ('triple_Utube_independent', np.array([0.2, 0.15, 0.10]), np.array([5., -1., 7.]), np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), np.array([-1274.69852553, 1468.59548314, -1617.86019738])),
+    ('triple_Utube_independent', np.array([0.2, 0.15, 0.10]), np.array([5., -1., 7.]), np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), np.array([-1255.0070539, 1487.29916766, -1601.09162509])),
     ])
 def test_fluid_heat_extraction_rate_independent(
         pipe_fixture, m_flow, T_f_in, segment_ratios, T_b, expected, request):
@@ -545,13 +545,13 @@ def test_fluid_heat_extraction_rate_independent(
 @pytest.mark.parametrize(
     "pipe_fixture, m_flow, T_f_in, segment_ratios, T_b, expected", [
     # Double U-tube
-    ('double_Utube_independent', np.array([0.2, 0.15]), np.array([5., -1.]), None, 1., -956.0096318353349),
-    ('double_Utube_independent', np.array([0.2, 0.15]), np.array([5., -1.]), None, np.array([1., 2., 3., 1.]), -355.0965406202829),
-    ('double_Utube_independent', np.array([0.2, 0.15]), np.array([5., -1.]), np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), -76.00884577783381),
+    ('double_Utube_independent', np.array([0.2, 0.15]), np.array([5., -1.]), None, 1., -956.0096318353369),
+    ('double_Utube_independent', np.array([0.2, 0.15]), np.array([5., -1.]), None, np.array([1., 2., 3., 1.]), -355.0965406202847),
+    ('double_Utube_independent', np.array([0.2, 0.15]), np.array([5., -1.]), np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), -32.03420904269501),
     # Triple U-tube
-    ('triple_Utube_independent', np.array([0.2, 0.15, 0.10]), np.array([5., -1., 7.]), None, 1., -2508.094081991294),
-    ('triple_Utube_independent', np.array([0.2, 0.15, 0.10]), np.array([5., -1., 7.]), None, np.array([1., 2., 3., 1.]), -1767.5943781810597),
-    ('triple_Utube_independent', np.array([0.2, 0.15, 0.10]), np.array([5., -1., 7.]), np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), -1423.9632397638252),
+    ('triple_Utube_independent', np.array([0.2, 0.15, 0.10]), np.array([5., -1., 7.]), None, 1., -2508.094081991296),
+    ('triple_Utube_independent', np.array([0.2, 0.15, 0.10]), np.array([5., -1., 7.]), None, np.array([1., 2., 3., 1.]), -1767.5943781810613),
+    ('triple_Utube_independent', np.array([0.2, 0.15, 0.10]), np.array([5., -1., 7.]), np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), -1368.7995113394977),
     ])
 def test_total_heat_extraction_rate_independent(
         pipe_fixture, m_flow, T_f_in, segment_ratios, T_b, expected, request):