Repository: MassimoCimmino/pygfunction (language: python)

Issue:
Enhance nSegment definition to be int or list
In every current description of `boreSegment`s in _pygfunction_, the number of segments is always declared as an integer and a uniform segment length is used across all segments. However, what if a user wants to define a number of segments on a per borehole basis? This issue seeks to address the current inability to define number of segments on a per borehole basis.  
Enhance segment definition
This is an enhancement to the definition of the variable `nSegments` to now allow a list in addition to a positive integer to be defined. If `nSegments` is defined as a list, then the length of the list must be equivalent to the number of boreholes and each borehole must have a segment length greater than or equal to 1. 

related to #133 
NewFeature unequal segment lengths
Each borehole should have a certain number of segments to be split into. Each segment should have a length. 

Issue #133 allows for boreholes to have unequal number of segments per borehole. Additionally, each segment length should be definable from a high level API. This issue seeks to address the definition of segment lengths. 

This new feature would allow for `SegmentLengths` of each segment to be defined in addition to the number of segments at each borehole being different. 
Enhance segment definition
This is an enhancement to the definition of the variable `nSegments` to now allow a list in addition to a positive integer to be defined. If `nSegments` is defined as a list, then the length of the list must be equivalent to the number of boreholes and each borehole must have a segment length greater than or equal to 1. 

related to #133 
NewFeature unequal segment lengths
Each borehole should have a certain number of segments to be split into. Each segment should have a length. 

Issue #133 allows for boreholes to have unequal number of segments per borehole. Additionally, each segment length should be definable from a high level API. This issue seeks to address the definition of segment lengths. 

This new feature would allow for `SegmentLengths` of each segment to be defined in addition to the number of segments at each borehole being different. 
NewFeature unequal segment lengths
Each borehole should have a certain number of segments to be split into. Each segment should have a length. 

Issue #133 allows for boreholes to have unequal number of segments per borehole. Additionally, each segment length should be definable from a high level API. This issue seeks to address the definition of segment lengths. 

This new feature would allow for `SegmentLengths` of each segment to be defined in addition to the number of segments at each borehole being different. 


Discussion / hints:

Fix patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index a87b6531..68094f99 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -5,6 +5,7 @@
 ### New features 
 
 * [Issue 36](https://github.com/MassimoCimmino/pygfunction/issues/36) - Added a `Coaxial` class to the `pipes` module to model boreholes with coaxial pipes.
+* [Issue 135](https://github.com/MassimoCimmino/pygfunction/issues/135) - Added functionality for non-uniform discretization of the segments along the boreholes. This increases the accuracy of *g*-function calculations for the same number of segments when compared to a uniform discretization. Segment lengths are defined using the `segment_ratios` option in the `gFunction` class. A `discretize` function is added to the `utilities` module to generate borehole discretizations using an expanding mesh.
 * [Issue 146](https://github.com/MassimoCimmino/pygfunction/issues/146) - Added new solver `'equivalent'` to the `gFunction` class. This solver uses hierarchical agglomerative clustering to identify groups of boreholes that are expected to have similar borehole wall temperatures and heat extraction rates. Each group of boreholes is represented by a single equivalent borehole. The FLS solution is adapted to evaluate thermal interactions between groups of boreholes. This greatly reduces the number of evaluations of the FLS solution and the size of the system of equations to evaluate the g-function.
 
 ### Enhancements
diff --git a/doc/source/example_index.rst b/doc/source/example_index.rst
index 4cc7fca8..b77ffec0 100644
--- a/doc/source/example_index.rst
+++ b/doc/source/example_index.rst
@@ -15,6 +15,7 @@ Examples
    examples/uniform_temperature
    examples/comparison_gfunction_solvers
    examples/unequal_segments
+   examples/discretize_boreholes
    examples/load_aggregation
    examples/comparison_load_aggregation
    examples/fluid_temperature
diff --git a/doc/source/examples/discretize_boreholes.rst b/doc/source/examples/discretize_boreholes.rst
new file mode 100644
index 00000000..feac9fee
--- /dev/null
+++ b/doc/source/examples/discretize_boreholes.rst
@@ -0,0 +1,22 @@
+.. examples:
+
+********************************************************************
+Calculation of g-Functions computed with non-uniform segment lengths
+********************************************************************
+
+This example demonstrates the use of the :func:`utilities <../modules/utilities>`
+module to determine discretized segment ratios along a borehole.
+
+The following script computes g-Functions for a field of 6x4 boreholes
+utilizing the MIFT and UBWT boundary conditions with a 48 segments per borehole
+and equal segment lengths. The MIFT and UBWT g-functions are computed with only
+8 segments per borehole and non-uniform segment lengths. RMSE values are
+compared. It is shown that g-functions can be calculated accurately using a
+small number of segments.
+
+The script is located in:
+`pygfunction/examples/discretize_boreholes.py`
+
+.. literalinclude:: ../../../pygfunction/examples/discretize_boreholes.py
+   :language: python
+   :linenos:
\ No newline at end of file
diff --git a/pygfunction/boreholes.py b/pygfunction/boreholes.py
index 1b78425b..5aa7f288 100644
--- a/pygfunction/boreholes.py
+++ b/pygfunction/boreholes.py
@@ -84,14 +84,6 @@ def position(self):
         pos : tuple
             Position (x, y) (in meters) of the borehole.
 
-        Raises
-        ------
-        SomeError
-
-        See Also
-        --------
-        OtherModules
-
         Examples
         --------
         >>> b1 = gt.boreholes.Borehole(H=150., D=4., r_b=0.075, x=5., y=0.)
@@ -102,7 +94,7 @@ def position(self):
         pos = (self.x, self.y)
         return pos
 
-    def segments(self, nSegments):
+    def segments(self, nSegments, segment_ratios=None):
         """
         Split a borehole into segments.
 
@@ -110,6 +102,12 @@ def segments(self, nSegments):
         ----------
         nSegments : int
             Number of segments.
+        segment_ratios : array, optional
+            Ratio of the borehole length represented by each segment. The
+            sum of ratios must be equal to 1. The shape of the array is of
+            (nSegments,). If segment_ratios==None, segments of equal lengths
+            are considered.
+            Default is None.
 
         Returns
         -------
@@ -122,16 +120,19 @@ def segments(self, nSegments):
         >>> b1.segments(5)
 
         """
+        if segment_ratios is None:
+            segment_ratios = np.full(nSegments, 1. / nSegments)
+        z = self._segment_edges(nSegments, segment_ratios=segment_ratios)[:-1]
         boreSegments = []
         for i in range(nSegments):
             # Divide borehole into segments of equal length
-            H = self.H / nSegments
+            H = segment_ratios[i] * self.H
             # Buried depth of the i-th segment
-            D = self.D + i * self.H / nSegments * np.cos(self.tilt)
+            D = self.D + z[i] * np.cos(self.tilt)
             # x-position
-            x = self.x + i * self.H / nSegments * np.sin(self.tilt) * np.cos(self.orientation)
+            x = self.x + z[i] * np.sin(self.tilt) * np.cos(self.orientation)
             # y-position
-            y = self.y + i * self.H / nSegments * np.sin(self.tilt) * np.sin(self.orientation)
+            y = self.y + z[i] * np.sin(self.tilt) * np.sin(self.orientation)
             # Add to list of segments
             boreSegments.append(
                 Borehole(H, D, self.r_b, x, y,
@@ -139,6 +140,73 @@ def segments(self, nSegments):
                          orientation=self.orientation))
         return boreSegments
 
+    def _segment_edges(self, nSegments, segment_ratios=None):
+        """
+        Linear coordinates of the segment edges.
+
+        Parameters
+        ----------
+        nSegments : int
+            Number of segments.
+        segment_ratios : array, optional
+            Ratio of the borehole length represented by each segment. The
+            sum of ratios must be equal to 1. The shape of the array is of
+            (nSegments,). If segment_ratios==None, segments of equal lengths
+            are considered.
+            Default is None.
+
+        Returns
+        -------
+        z : array
+            Coordinates of the segment edges, with z=0. corresponding to the
+            borehole head and z=H corresponding to the bottom edge of the
+            borehole. The shape of the array is of (nSegments + 1,).
+
+        Examples
+        --------
+        >>> b1 = gt.boreholes.Borehole(H=150., D=4., r_b=0.075, x=5., y=0.)
+        >>> b1._segment_edges(5)
+
+        """
+        if segment_ratios is None:
+            segment_ratios = np.full(nSegments, 1. / nSegments)
+        z = np.concatenate(([0.], np.cumsum(segment_ratios))) * self.H
+        return z
+
+    def _segment_midpoints(self, nSegments, segment_ratios=None):
+        """
+        Linear coordinates of the segment midpoints.
+
+        Parameters
+        ----------
+        nSegments : int
+            Number of segments.
+        segment_ratios : array, optional
+            Ratio of the borehole length represented by each segment. The
+            sum of ratios must be equal to 1. The shape of the array is of
+            (nSegments,). If segment_ratios==None, segments of equal lengths
+            are considered.
+            Default is None.
+
+        Returns
+        -------
+        z : array
+            Coordinates of the segment midpoints, with z=0. corresponding to
+            the borehole head and z=H corresponding to the bottom edge of the
+            borehole. The shape of the array is of (nSegments,).
+
+        Examples
+        --------
+        >>> b1 = gt.boreholes.Borehole(H=150., D=4., r_b=0.075, x=5., y=0.)
+        >>> b1._segment_midpoints(5)
+
+        """
+        if segment_ratios is None:
+            segment_ratios = np.full(nSegments, 1. / nSegments)
+        z = self._segment_edges(nSegments, segment_ratios=segment_ratios)[:-1] \
+            + segment_ratios * self.H / 2
+        return z
+
 
 class _EquivalentBorehole(object):
     """
@@ -253,7 +321,7 @@ def position(self):
         """
         return (self.x, self.y)
 
-    def segments(self, nSegments):
+    def segments(self, nSegments, segment_ratios=None):
         """
         Split an equivalent borehole into segments.
 
@@ -261,7 +329,12 @@ def segments(self, nSegments):
         ----------
         nSegments : int
             Number of segments.
-
+        segment_ratios : array, optional
+            Ratio of the borehole length represented by each segment. The
+            sum of ratios must be equal to 1. The shape of the array is of
+            (nSegments,). If segment_ratios==None, segments of equal lengths
+            are considered.
+            Default is None.
         Returns
         -------
         boreSegments : list
@@ -273,7 +346,10 @@ def segments(self, nSegments):
         >>> b1.segments(5)
 
         """
-        return [_EquivalentBorehole((self.H/nSegments, self.D+i*self.H/nSegments, self.r_b, self.x, self.y)) for i in range(nSegments)]
+        if segment_ratios is None:
+            segment_ratios = np.full(nSegments, 1. / nSegments)
+        z = self._segment_edges(nSegments, segment_ratios=segment_ratios)[:-1]
+        return [_EquivalentBorehole((segment_ratios[i] * self.H, z[i] + self.D, self.r_b, self.x, self.y)) for i in range(nSegments)]
 
     def unique_distance(self, target, disTol=0.01):
         """
@@ -337,6 +413,73 @@ def unique_distance(self, target, disTol=0.01):
         
         return np.array(dis), np.array(wDis)
 
+    def _segment_edges(self, nSegments, segment_ratios=None):
+        """
+        Linear coordinates of the segment edges.
+
+        Parameters
+        ----------
+        nSegments : int
+            Number of segments.
+        segment_ratios : array, optional
+            Ratio of the borehole length represented by each segment. The
+            sum of ratios must be equal to 1. The shape of the array is of
+            (nSegments,). If segment_ratios==None, segments of equal lengths
+            are considered.
+            Default is None.
+
+        Returns
+        -------
+        z : array
+            Coordinates of the segment edges, with z=0. corresponding to the
+            borehole head and z=H corresponding to the bottom edge of the
+            borehole. The shape of the array is of (nSegments + 1,).
+
+        Examples
+        --------
+        >>> b1 = gt.boreholes.Borehole(H=150., D=4., r_b=0.075, x=5., y=0.)
+        >>> b1._segment_edges(5)
+
+        """
+        if segment_ratios is None:
+            segment_ratios = np.full(nSegments, 1. / nSegments)
+        z = np.concatenate(([0.], np.cumsum(segment_ratios))) * self.H
+        return z
+
+    def _segment_midpoints(self, nSegments, segment_ratios=None):
+        """
+        Linear coordinates of the segment midpoints.
+
+        Parameters
+        ----------
+        nSegments : int
+            Number of segments.
+        segment_ratios : array, optional
+            Ratio of the borehole length represented by each segment. The
+            sum of ratios must be equal to 1. The shape of the array is of
+            (nSegments,). If segment_ratios==None, segments of equal lengths
+            are considered.
+            Default is None.
+
+        Returns
+        -------
+        z : array
+            Coordinates of the segment midpoints, with z=0. corresponding to
+            the borehole head and z=H corresponding to the bottom edge of the
+            borehole. The shape of the array is of (nSegments,).
+
+        Examples
+        --------
+        >>> b1 = gt.boreholes.Borehole(H=150., D=4., r_b=0.075, x=5., y=0.)
+        >>> b1._segment_midpoints(5)
+
+        """
+        if segment_ratios is None:
+            segment_ratios = np.full(nSegments, 1. / nSegments)
+        z = self._segment_edges(nSegments, segment_ratios=segment_ratios)[:-1] \
+            + segment_ratios * self.H / 2
+        return z
+
 
 def find_duplicates(boreField, disp=False):
     """
diff --git a/pygfunction/examples/comparison_gfunction_solvers.py b/pygfunction/examples/comparison_gfunction_solvers.py
index 71024461..38cf92b8 100644
--- a/pygfunction/examples/comparison_gfunction_solvers.py
+++ b/pygfunction/examples/comparison_gfunction_solvers.py
@@ -40,7 +40,7 @@ def main():
     alpha = 1.0e-6      # Ground thermal diffusivity (m2/s)
 
     # g-Function calculation options
-    options = {'nSegments':12, 'disp':True}
+    options = {'nSegments': 12, 'disp': True}
 
     # Geometrically expanding time vector.
     dt = 100*3600.                  # Time step
@@ -76,7 +76,6 @@ def main():
     t3 = tic()
     t_equivalent = t3 - t2
 
-
     # -------------------------------------------------------------------------
     # Plot results
     # -------------------------------------------------------------------------
@@ -104,9 +103,10 @@ def main():
     ax.plot(lntts, np.abs(gfunc_equivalent.gFunc - gfunc_detailed.gFunc),
             '--', label='equivalent')
     ax.legend()
-    ax.set_title("Absolute error relative to the 'detailed' solver (Field of {} by {} boreholes)".format(N_1, N_2))
+    ax.set_title("Absolute error relative to the 'detailed' solver "
+                 "(Field of {} by {} boreholes)".format(N_1, N_2))
     # Adjust to plot window
-    plt.tight_layout()
+    fig.tight_layout()
 
     # Draw relative error
     # Configure figure and axes
@@ -117,14 +117,16 @@ def main():
     ax.set_ylabel(r'Relative error')
     gt.utilities._format_axes(ax)
     # Relative error
-    ax.plot(lntts, (gfunc_similarities.gFunc - gfunc_detailed.gFunc) / gfunc_detailed.gFunc,
+    gFunc_ref = gfunc_detailed.gFunc  # reference g-function
+    ax.plot(lntts, (gfunc_similarities.gFunc - gFunc_ref) / gFunc_ref,
             '-', label='similarities')
-    ax.plot(lntts, (gfunc_equivalent.gFunc - gfunc_detailed.gFunc) / gfunc_detailed.gFunc,
+    ax.plot(lntts, (gfunc_equivalent.gFunc - gFunc_ref) / gFunc_ref,
             '--', label='equivalent')
     ax.legend()
-    ax.set_title("Relative error relative to the 'detailed' solver (Field of {} by {} boreholes)".format(N_1, N_2))
+    ax.set_title("Relative error relative to the 'detailed' solver "
+                 "(Field of {} by {} boreholes)".format(N_1, N_2))
     # Adjust to plot window
-    plt.tight_layout()
+    fig.tight_layout()
 
     # -------------------------------------------------------------------------
     # Borehole field (Second bore field)
@@ -147,7 +149,6 @@ def main():
     t3 = tic()
     t_equivalent = t3 - t2
 
-
     # -------------------------------------------------------------------------
     # Plot results
     # -------------------------------------------------------------------------
@@ -171,9 +172,10 @@ def main():
     ax.plot(lntts, np.abs(gfunc_equivalent.gFunc - gfunc_similarities.gFunc),
             label='equivalent')
     ax.legend()
-    ax.set_title("Absolute error relative to the 'similarities' solver (Field of {} by {} boreholes)".format(N_1, N_2))
+    ax.set_title("Absolute error relative to the 'similarities' solver "
+                 "(Field of {} by {} boreholes)".format(N_1, N_2))
     # Adjust to plot window
-    plt.tight_layout()
+    fig.tight_layout()
 
     # Draw relative error
     # Configure figure and axes
@@ -187,9 +189,10 @@ def main():
     ax.plot(lntts, (gfunc_equivalent.gFunc - gfunc_similarities.gFunc) / gfunc_similarities.gFunc,
             label='equivalent')
     ax.legend()
-    ax.set_title("Relative error relative to the 'similarities' solver (Field of {} by {} boreholes)".format(N_1, N_2))
+    ax.set_title("Relative error relative to the 'similarities' solver "
+                 "(Field of {} by {} boreholes)".format(N_1, N_2))
     # Adjust to plot window
-    plt.tight_layout()
+    fig.tight_layout()
 
     return
 
diff --git a/pygfunction/examples/discretize_boreholes.py b/pygfunction/examples/discretize_boreholes.py
new file mode 100644
index 00000000..67747186
--- /dev/null
+++ b/pygfunction/examples/discretize_boreholes.py
@@ -0,0 +1,177 @@
+# -*- coding: utf-8 -*-
+""" Example of g-function calculation using non-uniform segment lengths along
+    the boreholes.
+
+    The g-functions of a field of 6x4 boreholes are calculated for two
+    boundary conditions : (1) a uniform borehole wall temperature along the
+    boreholes equal for all boreholes, and (2) an equal inlet fluid
+    temperature into the boreholes. g-Functions using 8 segments in a
+    non-uniform deiscretization are compared to reference g-functions
+    calculated using 48 segments of equal lengths. It is shown that g-functions
+    can be calculated accurately using a small number of segments.
+"""
+
+import pygfunction as gt
+from numpy import pi
+import matplotlib.pyplot as plt
+import numpy as np
+
+
+def main():
+    # -------------------------------------------------------------------------
+    # Simulation parameters
+    # -------------------------------------------------------------------------
+
+    # Borehole dimensions
+    D = 4.0             # Borehole buried depth (m)
+    H = 150.0           # Borehole length (m)
+    r_b = 0.075         # Borehole radius (m)
+    B = 7.5             # Borehole spacing (m)
+
+    # Pipe dimensions
+    r_out = 0.0211      # Pipe outer radius (m)
+    r_in = 0.0147       # Pipe inner radius (m)
+    D_s = 0.052         # Shank spacing (m)
+    epsilon = 1.0e-6    # Pipe roughness (m)
+
+    # Pipe positions
+    # Single U-tube [(x_in, y_in), (x_out, y_out)]
+    pos_pipes = [(-D_s, 0.), (D_s, 0.)]
+
+    # Ground properties
+    alpha = 1.0e-6      # Ground thermal diffusivity (m2/s)
+    k_s = 2.0           # Ground thermal conductivity (W/m.K)
+
+    # Grout properties
+    k_g = 1.0           # Grout thermal conductivity (W/m.K)
+
+    # Pipe properties
+    k_p = 0.4           # Pipe thermal conductivity (W/m.K)
+
+    # Fluid properties
+    m_flow_borehole = 0.25  # Total fluid mass flow rate per borehole (kg/s)
+    # The fluid is propylene-glycol (20 %) at 20 degC
+    fluid = gt.media.Fluid('MPG', 20.)
+    cp_f = fluid.cp     # Fluid specific isobaric heat capacity (J/kg.K)
+    rho_f = fluid.rho   # Fluid density (kg/m3)
+    mu_f = fluid.mu     # Fluid dynamic viscosity (kg/m.s)
+    k_f = fluid.k       # Fluid thermal conductivity (W/m.K)
+
+    # g-Function calculation options
+
+    # Number of segments used in the reference calculation with uniform
+    # discretization
+    nSegments_uniform = 48
+    options_uniform = {'nSegments': nSegments_uniform,
+                       'disp': True}
+    # Number of segments used in the calculation with non-uniform
+    # discretization
+    nSegments_unequal = 8
+    segment_ratios = gt.utilities.segment_ratios(
+        nSegments_unequal, end_length_ratio=0.02)
+    options_unequal = {'nSegments': nSegments_unequal,
+                       'segment_ratios': segment_ratios,
+                       'disp':True}
+
+    # Geometrically expanding time vector.
+    dt = 100*3600.                  # Time step
+    tmax = 3000. * 8760. * 3600.    # Maximum time
+    Nt = 25                         # Number of time steps
+    ts = H**2/(9.*alpha)            # Bore field characteristic time
+    time = gt.utilities.time_geometric(dt, tmax, Nt)
+
+    # -------------------------------------------------------------------------
+    # Borehole field
+    # -------------------------------------------------------------------------
+
+    # Field of 6x4 (n=24) boreholes
+    N_1 = 6
+    N_2 = 4
+    boreField = gt.boreholes.rectangle_field(N_1, N_2, B, B, H, D, r_b)
+    gt.boreholes.visualize_field(boreField)
+    nBoreholes = len(boreField)
+
+    # -------------------------------------------------------------------------
+    # Initialize pipe model
+    # -------------------------------------------------------------------------
+
+    # Pipe thermal resistance
+    R_p = gt.pipes.conduction_thermal_resistance_circular_pipe(
+        r_in, r_out, k_p)
+    # Fluid to inner pipe wall thermal resistance (Single U-tube)
+    m_flow_pipe = m_flow_borehole
+    h_f = gt.pipes.convective_heat_transfer_coefficient_circular_pipe(
+        m_flow_pipe, r_in, mu_f, rho_f, k_f, cp_f, epsilon)
+    R_f = 1.0/(h_f*2*pi*r_in)
+
+    # Single U-tube, same for all boreholes in the bore field
+    UTubes = []
+    for borehole in boreField:
+        SingleUTube = gt.pipes.SingleUTube(
+            pos_pipes, r_in, r_out, borehole, k_s, k_g, R_f + R_p)
+        UTubes.append(SingleUTube)
+    m_flow_network = m_flow_borehole*nBoreholes
+
+    # Network of boreholes connected in parallel
+    network = gt.networks.Network(
+        boreField, UTubes, m_flow_network=m_flow_network, cp_f=cp_f)
+
+    # -------------------------------------------------------------------------
+    # Evaluate the g-functions for the borefield
+    # -------------------------------------------------------------------------
+
+    # Compute g-function for the converged MIFT case with equal number of
+    # segments per borehole, and equal segment lengths along the boreholes
+    gfunc_MIFT_uniform = gt.gfunction.gFunction(
+        network, alpha, time=time, boundary_condition='MIFT',
+        options=options_uniform)
+
+    # Calculate the g-function for uniform borehole wall temperature
+    gfunc_UBWT_uniform = gt.gfunction.gFunction(
+        boreField, alpha, time=time, boundary_condition='UBWT',
+        options=options_uniform)
+
+    # Compute g-function for the MIFT case with equal number of segments per
+    # borehole, and non-uniform segment lengths along the boreholes
+    gfunc_MIFT_unequal = gt.gfunction.gFunction(
+        network, alpha, time=time, boundary_condition='MIFT',
+        options=options_unequal)
+
+    # Calculate the g-function for uniform borehole wall temperature
+    gfunc_UBWT_unequal = gt.gfunction.gFunction(
+        boreField, alpha, time=time, boundary_condition='UBWT',
+        options=options_unequal)
+
+    # Compute the rmse between the reference cases and the discretized
+    # (predicted) cases
+    RMSE_MIFT = RMSE(gfunc_MIFT_uniform.gFunc, gfunc_MIFT_unequal.gFunc)
+    print('RMSE (MIFT) = {0:.5f}'.format(RMSE_MIFT))
+    RMSE_UBWT = RMSE(gfunc_UBWT_uniform.gFunc, gfunc_UBWT_unequal.gFunc)
+    print('RMSE (UBWT) = {0:.5f}'.format(RMSE_UBWT))
+
+    # -------------------------------------------------------------------------
+    # Plot g-functions
+    # -------------------------------------------------------------------------
+
+    ax = gfunc_MIFT_uniform.visualize_g_function().axes[0]
+    ax.plot(np.log(time / ts), gfunc_UBWT_uniform.gFunc)
+    ax.plot(np.log(time / ts), gfunc_MIFT_unequal.gFunc, 'o')
+    ax.plot(np.log(time / ts), gfunc_UBWT_unequal.gFunc, 'o')
+    ax.legend(
+        ['Equal inlet temperature (uniform segments)',
+         'Uniform borehole wall temperature (uniform segments)',
+         'Equal inlet temperature (non-uniform segments)',
+         'Uniform borehole wall temperature (non-uniform segments)'])
+    plt.tight_layout()
+
+    return
+
+
+def RMSE(reference, predicted):
+    rmse = np.linalg.norm(predicted - reference) / len(reference)
+    return rmse
+
+
+# Main function
+if __name__ == '__main__':
+    main()
diff --git a/pygfunction/examples/unequal_segments.py b/pygfunction/examples/unequal_segments.py
index 4aade5f1..941d97b3 100644
--- a/pygfunction/examples/unequal_segments.py
+++ b/pygfunction/examples/unequal_segments.py
@@ -47,14 +47,18 @@ def main():
     # Evaluate g-functions with different segment options
     # -------------------------------------------------------------------------
 
-    # Calculate g-function with equal number of segments
-    nSegments = 12
-    options = {'nSegments': nSegments, 'disp': True}
+    # Calculate g-function with equal number of segments for all boreholes and
+    # uniform segment lengths
+    nSegments = 24
+    options = {'nSegments': nSegments,
+               'disp': True}
 
     gfunc_equal = gt.gfunction.gFunction(
         boreField, alpha, time=time, options=options)
 
-    # Calculate g-function with unequal number of segments
+    # Calculate g-function with predefined number of segments for each
+    # borehole, the segment lengths will be uniform along each borehole, but
+    # the number of segments for the boreholes can be unequal.
 
     # Boreholes 12, 14 and 18 have more segments than the others and their
     # heat extraction rate profiles are plotted.
@@ -67,20 +71,51 @@ def main():
     gfunc_unequal = gt.gfunction.gFunction(
         boreField, alpha, time=time, options=options)
 
+    # Calculate g-function with equal number of segments for each borehole,
+    # unequal segment lengths along the length of the borehole defined by
+    # segment ratios. The segment ratios for each borehole are the same.
+
+    nSegments = 8
+    # Define the segment ratios for each borehole in each segment
+    # the segment lengths are defined top to bottom left to right
+    segment_ratios = np.array([0.05, 0.10, 0.10, 0.25, 0.25, 0.10, 0.10, 0.05])
+    options = {'nSegments': nSegments,
+               'segment_ratios': segment_ratios,
+               'disp': True,
+               'profiles': True}
+
+    g_func_predefined = gt.gfunction.gFunction(
+        boreField, alpha, time=time, options=options)
+
     # -------------------------------------------------------------------------
     # Plot g-functions
     # -------------------------------------------------------------------------
 
     ax = gfunc_equal.visualize_g_function().axes[0]
     ax.plot(np.log(time/ts), gfunc_unequal.gFunc, 'r-.')
-    ax.legend(['Equal number of segments', 'Unequal number of segments'])
+    ax.plot(np.log(time/ts), g_func_predefined.gFunc, 'k-.')
+    ax.legend(['Equal number of segments',
+               'Unequal number of segments',
+               'Unequal segment lengths'])
     plt.tight_layout()
 
     # Heat extraction rate profiles
-    gfunc_unequal.visualize_heat_extraction_rates(
+    fig = gfunc_unequal.visualize_heat_extraction_rates(
+        iBoreholes=[18, 12, 14])
+    fig.suptitle('Heat extraction rates (unequal number of segments)')
+    fig.tight_layout()
+    fig = g_func_predefined.visualize_heat_extraction_rates(
+        iBoreholes=[18, 12, 14])
+    fig.suptitle('Heat extraction rates (unequal segment lengths)')
+    fig.tight_layout()
+    fig = gfunc_unequal.visualize_heat_extraction_rate_profiles(
         iBoreholes=[18, 12, 14])
-    gfunc_unequal.visualize_heat_extraction_rate_profiles(
+    fig.suptitle('Heat extraction rate profiles (unequal number of segments)')
+    fig.tight_layout()
+    fig = g_func_predefined.visualize_heat_extraction_rate_profiles(
         iBoreholes=[18, 12, 14])
+    fig.suptitle('Heat extraction rate profiles (unequal segment lengths)')
+    fig.tight_layout()
 
     return
 
diff --git a/pygfunction/gfunction.py b/pygfunction/gfunction.py
index 0425440c..ca8b32d7 100644
--- a/pygfunction/gfunction.py
+++ b/pygfunction/gfunction.py
@@ -24,7 +24,6 @@ class gFunction(object):
     modeled as a series of finite line source segments, as proposed in
     [#gFunction-CimBer2014]_. Different boundary conditions and solution
     methods are implemented.
-    
 
     Attributes
     ----------
@@ -86,6 +85,12 @@ class gFunction(object):
                 Number of line segments used per borehole, or list of number of
                 line segments used for each borehole.
                 Default is 12.
+            segment_ratios : array or list of arrays, optional
+                Ratio of the borehole length represented by each segment. The
+                sum of ratios must be equal to 1. The shape of the array is of
+                (nSegments,) or list of (nSegments[i],). If
+                segment_ratios==None, segments of equal lengths are considered.
+                Default is None.
             disp : bool, optional
                 Set to true to print progression messages.
                 Default is False.
@@ -538,7 +543,10 @@ def _heat_extraction_rates(self, iBoreholes):
                 # heat extraction rate.
                 i0 = self.solver._i0Segments[i]
                 i1 = self.solver._i1Segments[i]
-                Q_t.append(np.mean(self.solver.Q_b[i0:i1,:], axis=0))
+                segment_ratios = self.solver.segment_ratios[i]
+                Q_t.append(
+                    np.sum(self.solver.Q_b[i0:i1,:]*segment_ratios[:,np.newaxis],
+                           axis=0))
         return Q_t
 
     def _heat_extraction_rate_profiles(self, time, iBoreholes):
@@ -590,11 +598,12 @@ def _heat_extraction_rate_profiles(self, time, iBoreholes):
                                     axis=1)(time).flatten()
                 if self.solver.nBoreSegments[i] > 1:
                     # Borehole length ratio at the mid-depth of each segment
-                    z_ratio = np.linspace(
-                        start=0.5/self.solver.nBoreSegments[i],
-                        stop=1-0.5/self.solver.nBoreSegments[i],
-                        num=self.solver.nBoreSegments[i])
-                    z.append(self.solver.boreholes[i].D + self.solver.boreholes[i].H*z_ratio)
+                    segment_ratios = self.solver.segment_ratios[i]
+                    z.append(
+                        self.solver.boreholes[i].D \
+                        + self.solver.boreholes[i]._segment_midpoints(
+                            self.solver.nBoreSegments[i],
+                            segment_ratios=segment_ratios))
                     Q_b.append(Q_bi)
                 else:
                     # If there is only one segment, the heat extraction rate is
@@ -633,7 +642,10 @@ def _temperatures(self, iBoreholes):
                 # borehole wall temperature.
                 i0 = self.solver._i0Segments[i]
                 i1 = self.solver._i1Segments[i]
-                T_b.append(np.mean(self.solver.T_b[i0:i1,:], axis=0))
+                segment_ratios = self.solver.segment_ratios[i]
+                T_b.append(
+                    np.sum(self.solver.T_b[i0:i1,:]*segment_ratios[:,np.newaxis],
+                           axis=0))
         return T_b
 
     def _temperature_profiles(self, time, iBoreholes):
@@ -696,11 +708,13 @@ def _temperature_profiles(self, time, iBoreholes):
                                     axis=1)(time).flatten()
                 if self.solver.nBoreSegments[i] > 1:
                     # Borehole length ratio at the mid-depth of each segment
-                    z_ratio = np.linspace(
-                        start=0.5/self.solver.nBoreSegments[i],
-                        stop=1-0.5/self.solver.nBoreSegments[i],
-                        num=self.solver.nBoreSegments[i])
-                    z.append(self.solver.boreholes[i].D + self.solver.boreholes[i].H*z_ratio)
+
+                    segment_ratios = self.solver.segment_ratios[i]
+                    z.append(
+                        self.solver.boreholes[i].D \
+                        + self.solver.boreholes[i]._segment_midpoints(
+                            self.solver.nBoreSegments[i],
+                            segment_ratios=segment_ratios))
                     T_b.append(T_bi)
                 else:
                     # If there is only one segment, the temperature is
@@ -770,7 +784,7 @@ def _check_inputs(self):
         assert type(self.method) is str and self.method in acceptable_methods, \
             "Method \'{}\' is not an acceptable method. \n" \
             "Please provide one of the following inputs : {}".format(
-                self.boundary_condition, acceptable_methods)
+                self.method, acceptable_methods)
         return
 
 
@@ -855,7 +869,8 @@ def uniform_heat_extraction(boreholes, time, alpha, use_similarities=True,
     return gFunc.gFunc
 
 
-def uniform_temperature(boreholes, time, alpha, nSegments=12, kind='linear',
+def uniform_temperature(boreholes, time, alpha, nSegments=12,
+                        segment_ratios=None, kind='linear',
                         use_similarities=True, disTol=0.01, tol=1.0e-6,
                         dtype=np.double, disp=False, **kwargs):
     """
@@ -878,6 +893,12 @@ def uniform_temperature(boreholes, time, alpha, nSegments=12, kind='linear',
         Number of line segments used per borehole, or list of number of
         line segments used for each borehole.
         Default is 12.
+    segment_ratios : array or list of arrays, optional
+        Ratio of the borehole length represented by each segment. The
+        sum of ratios must be equal to 1. The shape of the array is of
+        (nSegments,) or list of (nSegments[i],). If segment_ratios==None,
+        segments of equal lengths are considered.
+        Default is None.
     kind : string, optional
         Interpolation method used for segment-to-segment thermal response
         factors. See documentation for scipy.interpolate.interp1d.
@@ -927,6 +948,7 @@ def uniform_temperature(boreholes, time, alpha, nSegments=12, kind='linear',
     boundary_condition = 'UBWT'
     # Build options dict
     options = {'nSegments':nSegments,
+               'segment_ratios':segment_ratios,
                'disp':disp,
                'kind':kind,
                'disTol':disTol,
@@ -948,8 +970,9 @@ def uniform_temperature(boreholes, time, alpha, nSegments=12, kind='linear',
 
 def equal_inlet_temperature(
         boreholes, UTubes, m_flow_borehole, cp_f, time, alpha,
-        kind='linear', nSegments=12, use_similarities=True,
-        disTol=0.01, tol=1.0e-6, dtype=np.double, disp=False, **kwargs):
+        kind='linear', nSegments=12, segment_ratios=None,
+        use_similarities=True, disTol=0.01, tol=1.0e-6, dtype=np.double,
+        disp=False, **kwargs):
     """
     Evaluate the g-function with equal inlet fluid temperatures.
 
@@ -976,6 +999,12 @@ def equal_inlet_temperature(
         Number of line segments used per borehole, or list of number of
         line segments used for each borehole.
         Default is 12.
+    segment_ratios : array or list of arrays, optional
+        Ratio of the borehole length represented by each segment. The
+        sum of ratios must be equal to 1. The shape of the array is of
+        (nSegments,) or list of (nSegments[i],). If segment_ratios==None,
+        segments of equal lengths are considered.
+        Default is None.
     kind : string, optional
         Interpolation method used for segment-to-segment thermal response
         factors. See documentation for scipy.interpolate.interp1d.
@@ -1018,6 +1047,7 @@ def equal_inlet_temperature(
     boundary_condition = 'MIFT'
     # Build options dict
     options = {'nSegments':nSegments,
+               'segment_ratios':segment_ratios,
                'disp':disp,
                'kind':kind,
                'disTol':disTol,
@@ -1039,8 +1069,9 @@ def equal_inlet_temperature(
 
 def mixed_inlet_temperature(network, m_flow_network, cp_f,
                             time, alpha, kind='linear', nSegments=12,
-                            use_similarities=True, disTol=0.01, tol=1.0e-6,
-                            dtype=np.double, disp=False, **kwargs):
+                            segment_ratios=None, use_similarities=True,
+                            disTol=0.01, tol=1.0e-6, dtype=np.double,
+                            disp=False, **kwargs):
     """
     Evaluate the g-function with mixed inlet fluid temperatures.
 
@@ -1069,6 +1100,12 @@ def mixed_inlet_temperature(network, m_flow_network, cp_f,
         Number of line segments used per borehole, or list of number of
         line segments used for each borehole.
         Default is 12.
+    segment_ratios : array or list of arrays, optional
+        Ratio of the borehole length represented by each segment. The
+        sum of ratios must be equal to 1. The shape of the array is of
+        (nSegments,) or list of (nSegments[i],). If segment_ratios==None,
+        segments of equal lengths are considered.
+        Default is None.
     kind : string, optional
         Interpolation method used for segment-to-segment thermal response
         factors. See documentation for scipy.interpolate.interp1d.
@@ -1129,6 +1166,7 @@ def mixed_inlet_temperature(network, m_flow_network, cp_f,
     boundary_condition = 'MIFT'
     # Build options dict
     options = {'nSegments':nSegments,
+               'segment_ratios':segment_ratios,
                'disp':disp,
                'kind':kind,
                'disTol':disTol,
@@ -1177,6 +1215,12 @@ class _BaseSolver(object):
         Number of line segments used per borehole, or list of number of
         line segments used for each borehole.
         Default is 12.
+    segment_ratios : array or list of arrays, optional
+        Ratio of the borehole length represented by each segment. The
+        sum of ratios must be equal to 1. The shape of the array is of
+        (nSegments,) or list of (nSegments[i],). If segment_ratios==None,
+        segments of equal lengths are considered.
+        Default is None.
     disp : bool, optional
         Set to true to print progression messages.
         Default is False.
@@ -1195,19 +1239,35 @@ class _BaseSolver(object):
 
     """
     def __init__(self, boreholes, network, time, boundary_condition,
-                 nSegments=12, disp=False, profiles=False,
+                 nSegments=12, segment_ratios=None, disp=False, profiles=False,
                  kind='linear', dtype=np.double, **other_options):
         self.boreholes = boreholes
         self.network = network
         # Convert time to a 1d array
         self.time = np.atleast_1d(time).flatten()
         self.boundary_condition = boundary_condition
-        # Find indices of first and last segments along boreholes
         nBoreholes = len(self.boreholes)
+        # Format number of segments and segment ratios
         if type(nSegments) is int:
             self.nBoreSegments = [nSegments] * nBoreholes
         else:
             self.nBoreSegments = nSegments
+        if isinstance(segment_ratios, np.ndarray):
+            segment_ratios = [segment_ratios] * nBoreholes
+        elif segment_ratios is None:
+            segment_ratios = [np.full(n, 1./n) for n in self.nBoreSegments]
+        self.segment_ratios = segment_ratios
+        # Shortcut for segment_ratios comparisons
+        self._equal_segment_ratios = \
+            (np.all(np.array(self.nBoreSegments, dtype=np.uint) == self.nBoreSegments[0])
+             and np.all([np.allclose(segment_ratios, self.segment_ratios[0]) for segment_ratios in self.segment_ratios]))
+        # Boreholes with a uniform discretization
+        self._uniform_segment_ratios = [
+            np.allclose(segment_ratios,
+                        segment_ratios[0:1],
+                        rtol=1e-6)
+            for segment_ratios in self.segment_ratios]
+        # Find indices of first and last segments along boreholes
         self._i0Segments = [sum(self.nBoreSegments[0:i])
                             for i in range(nBoreholes)]
         self._i1Segments = [sum(self.nBoreSegments[0:(i + 1)])
@@ -1283,7 +1343,7 @@ def solve(self, time, alpha):
         # Segment lengths
         H_b = self.segment_lengths()
         if self.boundary_condition == 'MIFT':
-            Hb_individual = np.array([b.H/nSegments for (b, nSegments) in zip(self.boreholes, self.nBoreSegments) for i in range(nSegments)])
+            Hb_individual = np.array([b.H for b in self.boreSegments], dtype=self.dtype)
         H_tot = np.sum(H_b)
         if self.disp: print('Building and solving the system of equations ...',
                             end='')
@@ -1364,7 +1424,8 @@ def solve(self, time, alpha):
                     a_in, a_b = self.network.coefficients_borehole_heat_extraction_rate(
                             self.network.m_flow_network,
                             self.network.cp_f,
-                            self.nBoreSegments)
+                            self.nBoreSegments,
+                            segment_ratios=self.segment_ratios)
                     k_s = self.network.p[0].k_s
                     A = np.block(
                         [[h_dt,
@@ -1437,15 +1498,10 @@ def borehole_segments(self):
 
         """
         boreSegments = []  # list for storage of boreSegments
-        for j in range(len(self.boreholes)):
-            for i in range(self.nBoreSegments[j]):
-                b = self.boreholes[j]
-                # Divide borehole into segments of equal length
-                H_b = b.H / self.nBoreSegments[j]
-                # Buried depth of the i-th segment
-                D = b.D + i * b.H / self.nBoreSegments[j]
-                # Add to list of segments
-                boreSegments.append(Borehole(H_b, D, b.r_b, b.x, b.y))
+        for b, nSegments, segment_ratios in zip(self.boreholes, self.nBoreSegments, self.segment_ratios):
+            segments = b.segments(nSegments, segment_ratios=segment_ratios)
+            boreSegments.extend(segments)
+
         return boreSegments
 
     def temporal_superposition(self, h_ij, Q_reconstructed):
@@ -1571,6 +1627,16 @@ def _check_inputs(self):
             "Data type \'{}\' is not an acceptable data type. \n" \
             "Please provide one of the following inputs : {}".format(
                 self.dtype, acceptable_dtypes)
+        # Check segment ratios
+        for j in range(len(self.boreholes)):
+            assert len(self.segment_ratios[j]) == self.nBoreSegments[j], \
+                "The length of the segment ratios vectors must correspond to " \
+                "the number of segments, check borehole {}.".format(j)
+            error = np.abs(1. - np.sum(self.segment_ratios[j]))
+            assert(error < 1.0e-6), \
+                "Defined segment ratios must add up to 1. " \
+                ", check borehole {}.".format(j)
+
         return
 
 
@@ -1613,6 +1679,12 @@ class _Detailed(_BaseSolver):
         Number of line segments used per borehole, or list of number of
         line segments used for each borehole.
         Default is 12.
+    segment_ratios : array or list of arrays, optional
+        Ratio of the borehole length represented by each segment. The
+        sum of ratios must be equal to 1. The shape of the array is of
+        (nSegments,) or list of (nSegments[i],). If segment_ratios==None,
+        segments of equal lengths are considered.
+        Default is None.
     disp : bool, optional
         Set to true to print progression messages.
         Default is False.
@@ -1702,7 +1774,7 @@ def thermal_response_factors(self, time, alpha, kind='linear'):
 
         for i in range(nBoreholes):
             # Segments of the receiving borehole
-            b2 = self.boreholes[i].segments(self.nBoreSegments[i])
+            b2 = self.boreholes[i].segments(self.nBoreSegments[i], segment_ratios=self.segment_ratios[i])
             # -----------------------------------------------------------------
             # Segment-to-segment thermal response factors for same-borehole
             # thermal interactions
@@ -1722,7 +1794,7 @@ def thermal_response_factors(self, time, alpha, kind='linear'):
                 # Segments of the emitting borehole
                 b1 = [seg
                       for (b, iBor) in zip(self.boreholes[i+1:], range(i+1, nBoreholes))
-                      for seg in b.segments(self.nBoreSegments[iBor])]
+                      for seg in b.segments(self.nBoreSegments[iBor], segment_ratios=self.segment_ratios[iBor])]
                 h = finite_line_source(time, alpha, b1, b2)
                 # Broadcast values to h_ij matrix
                 for j in range(i+1, nBoreholes):
@@ -1790,6 +1862,12 @@ class _Similarities(_BaseSolver):
         Number of line segments used per borehole, or list of number of
         line segments used for each borehole.
         Default is 12.
+    segment_ratios : array or list of arrays, optional
+        Ratio of the borehole length represented by each segment. The
+        sum of ratios must be equal to 1. The shape of the array is of
+        (nSegments,) or list of (nSegments[i],). If segment_ratios==None,
+        segments of equal lengths are considered.
+        Default is None.
     disp : bool, optional
         Set to true to print progression messages.
         Default is False.
@@ -1914,6 +1992,8 @@ def thermal_response_factors(self, time, alpha, kind='linear'):
                 self._map_segment_pairs(
                     i_pair, j_pair, k_pair, [(n, n) for n in group], [0])
             # Evaluate FLS at all time steps
+            H1 = H1.reshape(1, -1)
+            H2 = H2.reshape(1, -1)
             D1 = D1.reshape(1, -1)
             D2 = D2.reshape(1, -1)
             dis = self.boreholes[i].r_b
@@ -1938,17 +2018,20 @@ def thermal_response_factors(self, time, alpha, kind='linear'):
                     self.borehole_to_borehole_indices[n])
             # Evaluate FLS at all time steps
             dis = np.reshape(self.borehole_to_borehole_distances[n], (-1, 1))
+            H1 = H1.reshape(1, -1)
+            H2 = H2.reshape(1, -1)
             D1 = D1.reshape(1, -1)
             D2 = D2.reshape(1, -1)
             h = finite_line_source_vectorized(time, alpha, dis, H1, D1, H2, D2)
             # Broadcast values to h_ij matrix
             h_ij[j_segment, i_segment, 1:] = h[l_segment, k_segment, :]
             if (self._compare_boreholes(self.boreholes[j], self.boreholes[i]) and
-                self.nBoreSegments[i] == self.nBoreSegments[j]):
+                self.nBoreSegments[i] == self.nBoreSegments[j] and
+                self._uniform_segment_ratios[i] and
+                self._uniform_segment_ratios[j]):
                 h_ij[i_segment, j_segment, 1:] = h[l_segment, k_segment, :]
             else:
-                h_ij[i_segment, j_segment, 1:] = (h[l_segment, k_segment, :].T \
-                    * segment_lengths[j_segment]/segment_lengths[i_segment]).T
+                h_ij[i_segment, j_segment, 1:] = (h * H2.T / H1.T)[l_segment, k_segment, :]
 
         # Return 2d array if time is a scalar
         if np.isscalar(time):
@@ -2146,7 +2229,9 @@ def _find_axial_borehole_pairs(self, boreholes):
                     # Add the borehole to the group if a similar borehole is
                     # found
                     if (self._compare_boreholes(boreholes[i], boreholes[m]) and
-                        self.nBoreSegments[i] == self.nBoreSegments[m]):
+                        (self._equal_segment_ratios or
+                         (self.nBoreSegments[i] == self.nBoreSegments[m] and
+                          np.allclose(self.segment_ratios[i], self.segment_ratios[m], rtol=self.tol)))):
                         borehole_to_self[k].append(i)
                         break
                 else:
@@ -2162,13 +2247,19 @@ def _find_axial_borehole_pairs(self, boreholes):
                         # Add the pair (or the reciprocal pair) to a group
                         # if a similar one is found
                         if (compare_pairs(pair0, pair_ref) and
-                            self.nBoreSegments[i] == self.nBoreSegments[m] and
-                            self.nBoreSegments[j] == self.nBoreSegments[n]):
+                            (self._equal_segment_ratios or
+                             (self.nBoreSegments[i] == self.nBoreSegments[m] and
+                              self.nBoreSegments[j] == self.nBoreSegments[n] and
+                              np.allclose(self.segment_ratios[i], self.segment_ratios[m], rtol=self.tol) and
+                              np.allclose(self.segment_ratios[j], self.segment_ratios[n], rtol=self.tol)))):
                             borehole_to_borehole[k].append((i, j))
                             break
                         elif (compare_pairs(pair1, pair_ref) and
-                              self.nBoreSegments[j] == self.nBoreSegments[m] and
-                              self.nBoreSegments[i] == self.nBoreSegments[n]):
+                              (self._equal_segment_ratios or
+                               (self.nBoreSegments[j] == self.nBoreSegments[m] and
+                                self.nBoreSegments[i] == self.nBoreSegments[n] and
+                                np.allclose(self.segment_ratios[j], self.segment_ratios[m], rtol=self.tol) and
+                                np.allclose(self.segment_ratios[i], self.segment_ratios[n], rtol=self.tol)))):
                             borehole_to_borehole[k].append((j, i))
                             break
                     # If no similar pairs are known, append the groups
@@ -2266,11 +2357,11 @@ def _map_axial_segment_pairs(self, i, j,
 
         Returns
         -------
-        H1 : float
+        H1 : array
             Length of the emitting segments.
         D1 : array
             Array of buried depths of the emitting segments.
-        H2 : float
+        H2 : array
             Length of the receiving segments.
         D2 : array
             Array of buried depths of the receiving segments.
@@ -2298,11 +2389,11 @@ def _map_axial_segment_pairs(self, i, j,
             # Find segment pairs for the image FLS solution
             compare_pairs = self._compare_image_pairs
         # Dive both boreholes into segments
-        segments1 = borehole1.segments(self.nBoreSegments[i])
-        segments2 = borehole2.segments(self.nBoreSegments[j])
-        # Segments have equal lengths
-        H1 = segments1[0].H
-        H2 = segments2[0].H
+        segments1 = borehole1.segments(self.nBoreSegments[i], segment_ratios=self.segment_ratios[i])
+        segments2 = borehole2.segments(self.nBoreSegments[j], segment_ratios=self.segment_ratios[j])
+        # Prepare lists of segment lengths
+        H1 = []
+        H2 = []
         # Prepare lists of segment buried depths
         D1 = []
         D2 = []
@@ -2340,12 +2431,14 @@ def _map_axial_segment_pairs(self, i, j,
                 # depths
                 else:
                     k_pair[p] = nPairs
+                    H1.append(segments1[ii].H)
+                    H2.append(segments2[jj].H)
                     D1.append(segments1[ii].D)
                     D2.append(segments2[jj].D)
                     unique_pairs.append((ii, jj))
                     nPairs += 1
                 p += 1
-        return H1, np.array(D1), H2, np.array(D2), i_pair, j_pair, k_pair
+        return np.array(H1), np.array(D1), np.array(H2), np.array(D2), i_pair, j_pair, k_pair
 
     def _map_segment_pairs(self, i_pair, j_pair, k_pair, borehole_to_borehole,
                            borehole_to_borehole_indices):
@@ -2452,6 +2545,12 @@ class _Equivalent(_BaseSolver):
         Number of line segments used per borehole, or list of number of
         line segments used for each borehole.
         Default is 12.
+    segment_ratios : array or list of arrays, optional
+        Ratio of the borehole length represented by each segment. The
+        sum of ratios must be equal to 1. The shape of the array is of
+        (nSegments,) or list of (nSegments[i],). If segment_ratios==None,
+        segments of equal lengths are considered.
+        Default is None.
     disp : bool, optional
         Set to true to print progression messages.
         Default is False.
@@ -2522,6 +2621,8 @@ def initialize(self, disTol=0.01, tol=1.0e-6, kClusters=1, **kwargs):
         # Initialize groups for equivalent boreholes
         nSources = self.find_groups()
         self.nBoreSegments = [self.nBoreSegments[0]] * self.nEqBoreholes
+        self.segment_ratios = [self.segment_ratios[0]] * self.nEqBoreholes
+        self.boreSegments = self.borehole_segments()
         self._i0Segments = [sum(self.nBoreSegments[0:i]) 
                             for i in range(self.nEqBoreholes)]
         self._i1Segments = [sum(self.nBoreSegments[0:(i + 1)])
@@ -2568,6 +2669,7 @@ def thermal_response_factors(self, time, alpha, kind='linear'):
         tic = tim.time()
         # Initialize segment-to-segment response factors
         h_ij = np.zeros((self.nSources, self.nSources, nt+1), dtype=self.dtype)
+        segment_lengths = self.segment_lengths()
 
         # ---------------------------------------------------------------------
         # Segment-to-segment thermal response factors for borehole-to-borehole
@@ -2583,6 +2685,8 @@ def thermal_response_factors(self, time, alpha, kind='linear'):
             dis, wDis = self._find_unique_distances(self.dis, indices)
             H1, D1, H2, D2, i_pair, j_pair, k_pair = \
                 self._map_axial_segment_pairs(i, j)
+            H1 = H1.reshape(1, -1)
+            H2 = H2.reshape(1, -1)
             D1 = D1.reshape(1, -1)
             D2 = D2.reshape(1, -1)
             N2 = np.array([[self.boreholes[j].nBoreholes for (i, j) in indices]]).T
@@ -2596,9 +2700,8 @@ def thermal_response_factors(self, time, alpha, kind='linear'):
                 j_segment = self._i0Segments[j] + j_pair
                 h_ij[j_segment, i_segment, 1:] = h[k, k_pair, :]
                 if not i == j:
-                    w_ratio = self.wBoreholes[j] / self.wBoreholes[i]
-                    H_ratio = self.boreholes[j].H / self.boreholes[i].H
-                    h_ij[i_segment, j_segment, 1:] = h[k, k_pair, :] * H_ratio * w_ratio
+                    h_ij[i_segment, j_segment, 1:] = (h[k, k_pair, :].T \
+                        * segment_lengths[j_segment]/segment_lengths[i_segment]).T
 
         # ---------------------------------------------------------------------
         # Segment-to-segment thermal response factors for same-borehole thermal
@@ -2613,6 +2716,8 @@ def thermal_response_factors(self, time, alpha, kind='linear'):
                 self._map_axial_segment_pairs(i, i)
             # Evaluate FLS at all time steps
             dis = self.boreholes[i].r_b
+            H1 = H1.reshape(1, -1)
+            H2 = H2.reshape(1, -1)
             D1 = D1.reshape(1, -1)
             D2 = D2.reshape(1, -1)
             h = finite_line_source_vectorized(time, alpha, dis, H1, D1, H2, D2)
@@ -2733,7 +2838,8 @@ def find_groups(self, tol=1e-6):
                 pipes,
                 m_flow_network=self.network.m_flow_network,
                 cp_f=self.network.cp_f,
-                nSegments=self.nBoreSegments[0])
+                nSegments=self.nBoreSegments[0],
+                segment_ratios=self.segment_ratios[0])
 
         # Stop chrono
         toc = tim.time()
@@ -2759,8 +2865,8 @@ def segment_lengths(self):
         """
         # Borehole lengths
         H = np.array([seg.H*seg.nBoreholes
-                      for (b, nSegments) in zip(self.boreholes, self.nBoreSegments)
-                      for seg in b.segments(nSegments)],
+                      for (b, nSegments, segment_ratios) in zip(self.boreholes, self.nBoreSegments, self.segment_ratios)
+                      for seg in b.segments(nSegments, segment_ratios=segment_ratios)],
                      dtype=self.dtype)
         return H
 
@@ -3077,11 +3183,11 @@ def _map_axial_segment_pairs(self, iBor, jBor,
             # Find segment pairs for the image FLS solution
             compare_pairs = self._compare_image_pairs
         # Dive both boreholes into segments
-        segments1 = borehole1.segments(self.nBoreSegments[iBor])
-        segments2 = borehole2.segments(self.nBoreSegments[jBor])
-        # Segments have equal lengths
-        H1 = segments1[0].H
-        H2 = segments2[0].H
+        segments1 = borehole1.segments(self.nBoreSegments[iBor], segment_ratios=self.segment_ratios[iBor])
+        segments2 = borehole2.segments(self.nBoreSegments[jBor], segment_ratios=self.segment_ratios[jBor])
+        # Prepare lists of segment lengths
+        H1 = []
+        H2 = []
         # Prepare lists of segment buried depths
         D1 = []
         D2 = []
@@ -3119,12 +3225,14 @@ def _map_axial_segment_pairs(self, iBor, jBor,
                 # depths
                 else:
                     k_pair[p] = nPairs
+                    H1.append(segments1[i].H)
+                    H2.append(segments2[j].H)
                     D1.append(segments1[i].D)
                     D2.append(segments2[j].D)
                     unique_pairs.append((i, j))
                     nPairs += 1
                 p += 1
-        return H1, np.array(D1), H2, np.array(D2), i_pair, j_pair, k_pair
+        return np.array(H1), np.array(D1), np.array(H2), np.array(D2), i_pair, j_pair, k_pair
 
     def _check_solver_specific_inputs(self):
         """
@@ -3140,6 +3248,8 @@ def _check_solver_specific_inputs(self):
             "The precision increment 'kClusters' should be a positive int."
         assert np.all(np.array(self.nBoreSegments, dtype=np.uint) == self.nBoreSegments[0]), \
             "Solver 'equivalent' can only handle equal numbers of segments."
+        assert np.all([np.allclose(segment_ratios, self.segment_ratios[0]) for segment_ratios in self.segment_ratios]), \
+            "Solver 'equivalent' can only handle identical segment_ratios for all boreholes."
         if self.boundary_condition == 'MIFT':
             assert np.all(np.array(self.network.c, dtype=np.int) == -1), \
                 "Solver 'equivalent' is only valid for parallel-connected " \
diff --git a/pygfunction/heat_transfer.py b/pygfunction/heat_transfer.py
index 99ec1a11..c10bf0fc 100644
--- a/pygfunction/heat_transfer.py
+++ b/pygfunction/heat_transfer.py
@@ -110,10 +110,10 @@ def finite_line_source(
         elif isinstance(time, (np.floating, float)):
             # Lower bound of integration
             a = 1.0 / np.sqrt(4.0*alpha*time)
-            h = quad(f, a, np.inf)[0]
+            h = 0.5 / H2 * quad(f, a, np.inf)[0]
         else:
             h = np.stack(
-                [quad(f, 1.0 / np.sqrt(4.0*alpha*t), np.inf)[0] for t in time],
+                [0.5 / H2 * quad(f, 1.0 / np.sqrt(4.0*alpha*t), np.inf)[0] for t in time],
                 axis=-1)
     else:
         # Unpack parameters
@@ -237,10 +237,10 @@ def finite_line_source_vectorized(
     if isinstance(time, (np.floating, float)):
         # Lower bound of integration
         a = 1.0 / np.sqrt(4.0*alpha*time)
-        h = quad_vec(f, a, np.inf)[0]
+        h = 0.5 / H2 * quad_vec(f, a, np.inf)[0]
     else:
         h = np.stack(
-            [quad_vec(f, 1.0 / np.sqrt(4.0*alpha*t), np.inf)[0]
+            [0.5 / H2 * quad_vec(f, 1.0 / np.sqrt(4.0*alpha*t), np.inf)[0]
              for t in time],
             axis=-1)
     return h
@@ -301,6 +301,8 @@ def finite_line_source_equivalent_boreholes_vectorized(
         Lengths of the receiving heat sources.
     D2 : float or array
         Buried depths of the receiving heat sources.
+    N2 : float or array,
+        Number of segments represented by the receiving heat sources.
     reaSource : bool
         True if the real part of the FLS solution is to be included.
         Default is True.
@@ -341,10 +343,10 @@ def finite_line_source_equivalent_boreholes_vectorized(
     if isinstance(time, (np.floating, float)):
         # Lower bound of integration
         a = 1.0 / np.sqrt(4.0*alpha*time)
-        h = quad_vec(f, a, np.inf)[0]
+        h = 0.5 / (N2*H2) * quad_vec(f, a, np.inf)[0]
     else:
         h = np.stack(
-            [quad_vec(f, 1.0 / np.sqrt(4.0*alpha*t), np.inf)[0]
+            [0.5 / (N2*H2) * quad_vec(f, 1.0 / np.sqrt(4.0*alpha*t), np.inf)[0]
              for t in time],
             axis=-1)
     return h
@@ -412,7 +414,7 @@ def _finite_line_source_integrand(dis, H1, D1, H2, D2, reaSource, imgSource):
                       D2 + D1 + H1,
                       D2 + D1 + H2 + H1],
                      axis=-1)
-        f = lambda s: 0.5 / (H2*s**2) * np.exp(-dis**2*s**2) * np.inner(p, _erfint(q*s))
+        f = lambda s: s**-2 * np.exp(-dis**2*s**2) * np.inner(p, _erfint(q*s))
     elif reaSource:
         # Real FLS solution
         p = np.array([1, -1, 1, -1])
@@ -421,7 +423,7 @@ def _finite_line_source_integrand(dis, H1, D1, H2, D2, reaSource, imgSource):
                       D2 - D1 - H1,
                       D2 - D1 + H2 - H1],
                      axis=-1)
-        f = lambda s: 0.5 / (H2*s**2) * np.exp(-dis**2*s**2) * np.inner(p, _erfint(q*s))
+        f = lambda s: s**-2 * np.exp(-dis**2*s**2) * np.inner(p, _erfint(q*s))
     elif imgSource:
         # Image FLS solution
         p = np.array([1, -1, 1, -1])
@@ -430,7 +432,7 @@ def _finite_line_source_integrand(dis, H1, D1, H2, D2, reaSource, imgSource):
                       D2 + D1 + H1,
                       D2 + D1 + H2 + H1],
                      axis=-1)
-        f = lambda s: 0.5 / (H2*s**2) * np.exp(-dis**2*s**2) * np.inner(p, _erfint(q*s))
+        f = lambda s: s**-2 * np.exp(-dis**2*s**2) * np.inner(p, _erfint(q*s))
     else:
         # No heat source
         f = lambda s: 0.
@@ -455,7 +457,8 @@ def _finite_line_source_equivalent_boreholes_integrand(dis, wDis, H1, D1, H2, D2
         Lengths of the receiving heat sources.
     D2 : float or array
         Buried depths of the receiving heat sources.
-    N2 :
+    N2 : float or array,
+        Number of segments represented by the receiving heat sources.
     reaSource : bool
         True if the real part of the FLS solution is to be included.
     imgSource : bool
@@ -484,7 +487,7 @@ def _finite_line_source_equivalent_boreholes_integrand(dis, wDis, H1, D1, H2, D2
                       D2 + D1 + H1,
                       D2 + D1 + H2 + H1],
                      axis=-1)
-        f = lambda s: 0.5 / (N2*H2*s**2) * (np.exp(-dis**2*s**2) @ wDis).T * np.inner(p, _erfint(q*s))
+        f = lambda s: s**-2 * (np.exp(-dis**2*s**2) @ wDis).T * np.inner(p, _erfint(q*s))
     elif reaSource:
         # Real FLS solution
         p = np.array([1, -1, 1, -1])
@@ -493,7 +496,7 @@ def _finite_line_source_equivalent_boreholes_integrand(dis, wDis, H1, D1, H2, D2
                       D2 - D1 - H1,
                       D2 - D1 + H2 - H1],
                      axis=-1)
-        f = lambda s: 0.5 / (N2*H2*s**2) * (np.exp(-dis**2*s**2) @ wDis).T * np.inner(p, _erfint(q*s))
+        f = lambda s: s**-2 * (np.exp(-dis**2*s**2) @ wDis).T * np.inner(p, _erfint(q*s))
     elif imgSource:
         # Image FLS solution
         p = np.array([1, -1, 1, -1])
@@ -502,7 +505,7 @@ def _finite_line_source_equivalent_boreholes_integrand(dis, wDis, H1, D1, H2, D2
                       D2 + D1 + H1,
                       D2 + D1 + H2 + H1],
                      axis=-1)
-        f = lambda s: 0.5 / (N2*H2*s**2) * (np.exp(-dis**2*s**2) @ wDis).T * np.inner(p, _erfint(q*s))
+        f = lambda s: s**-2 * (np.exp(-dis**2*s**2) @ wDis).T * np.inner(p, _erfint(q*s))
     else:
         # No heat source
         f = lambda s: 0.
diff --git a/pygfunction/networks.py b/pygfunction/networks.py
index f17e5976..f4a86034 100644
--- a/pygfunction/networks.py
+++ b/pygfunction/networks.py
@@ -39,6 +39,12 @@ class Network(object):
         Number of line segments used per borehole. This parameter is used to
         initialize the coefficients if it is provided.
         Default is None.
+    segment_ratios :
+    (nSegments,) array or list of (nSegments[i],) arrays, 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.
 
     Notes
     -----
@@ -62,7 +68,8 @@ class Network(object):
 
     """
     def __init__(self, boreholes, pipes, bore_connectivity=None,
-                 m_flow_network=None, cp_f=None, nSegments=None):
+                 m_flow_network=None, cp_f=None, nSegments=None,
+                 segment_ratios=None):
         self.b = boreholes
         self.H_tot = sum([b.H for b in self.b])
         self.nBoreholes = len(boreholes)
@@ -89,10 +96,12 @@ def __init__(self, boreholes, pipes, bore_connectivity=None,
 
         # Initialize stored_coefficients
         self._initialize_coefficients_connectivity()
-        self._initialize_stored_coefficients(m_flow_network, cp_f, nSegments)
+        self._initialize_stored_coefficients(
+            m_flow_network, cp_f, nSegments, segment_ratios)
 
     def get_inlet_temperature(
-            self, T_f_in, T_b, m_flow_network, cp_f, nSegments):
+            self, T_f_in, T_b, m_flow_network, cp_f, nSegments,
+            segment_ratios=None):
         """
         Returns the inlet fluid temperatures of all boreholes.
 
@@ -113,6 +122,12 @@ def get_inlet_temperature(
             Number of borehole segments for each borehole. If an int is
             supplied, all boreholes are considered to have the same number of
             segments.
+        segment_ratios :
+        (nSegments,) array or list of (nSegments[i],) arrays, 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
         -------
@@ -122,14 +137,16 @@ def get_inlet_temperature(
         """
         # Build coefficient matrices
         a_in, a_b = self.coefficients_inlet_temperature(
-                m_flow_network, cp_f, nSegments)
+                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_f_in_borehole = a_in @ np.atleast_1d(T_f_in) + a_b @ T_b
         return T_f_in_borehole
 
-    def get_outlet_temperature(self, T_f_in, T_b, m_flow_network, cp_f, nSegments):
+    def get_outlet_temperature(
+            self, T_f_in, T_b, m_flow_network, cp_f, nSegments,
+            segment_ratios=None):
         """
         Returns the outlet fluid temperatures of all boreholes.
 
@@ -150,6 +167,12 @@ def get_outlet_temperature(self, T_f_in, T_b, m_flow_network, cp_f, nSegments):
             Number of borehole segments for each borehole. If an int is
             supplied, all boreholes are considered to have the same number of
             segments.
+        segment_ratios :
+        (nSegments,) array or list of (nSegments[i],) arrays, 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
         -------
@@ -159,7 +182,7 @@ def get_outlet_temperature(self, T_f_in, T_b, m_flow_network, cp_f, nSegments):
         """
         # Build coefficient matrices
         a_in, a_b = self.coefficients_outlet_temperature(
-                m_flow_network, cp_f, nSegments)
+                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))
@@ -167,7 +190,8 @@ def get_outlet_temperature(self, T_f_in, T_b, m_flow_network, cp_f, nSegments):
         return T_f_out
 
     def get_borehole_heat_extraction_rate(
-            self, T_f_in, T_b, m_flow_network, cp_f, nSegments):
+            self, T_f_in, T_b, m_flow_network, cp_f, nSegments,
+            segment_ratios=None):
         """
         Returns the heat extraction rates of all boreholes.
 
@@ -188,6 +212,12 @@ def get_borehole_heat_extraction_rate(
             Number of borehole segments for each borehole. If an int is
             supplied, all boreholes are considered to have the same number of
             segments.
+        segment_ratios :
+        (nSegments,) array or list of (nSegments[i],) arrays, 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
         -------
@@ -196,7 +226,7 @@ def get_borehole_heat_extraction_rate(
 
         """
         a_in, a_b = self.coefficients_borehole_heat_extraction_rate(
-                m_flow_network, cp_f, nSegments)
+                m_flow_network, cp_f, nSegments, segment_ratios=segment_ratios)
         if np.isscalar(T_b):
             T_b = np.tile(T_b, sum(self.nSegments))
         Q_b = a_in @ np.atleast_1d(T_f_in) + a_b @ T_b
@@ -204,7 +234,8 @@ def get_borehole_heat_extraction_rate(
         return Q_b
 
     def get_fluid_heat_extraction_rate(
-            self, T_f_in, T_b, m_flow_network, cp_f, nSegments):
+            self, T_f_in, T_b, m_flow_network, cp_f, nSegments,
+            segment_ratios=None):
         """
         Returns the total heat extraction rates of all boreholes.
 
@@ -225,6 +256,12 @@ def get_fluid_heat_extraction_rate(
             Number of borehole segments for each borehole. If an int is
             supplied, all boreholes are considered to have the same number of
             segments.
+        segment_ratios :
+        (nSegments,) array or list of (nSegments[i],) arrays, 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
         -------
@@ -233,7 +270,7 @@ def get_fluid_heat_extraction_rate(
 
         """
         a_in, a_b = self.coefficients_fluid_heat_extraction_rate(
-                m_flow_network, cp_f, nSegments)
+                m_flow_network, cp_f, nSegments, segment_ratios=segment_ratios)
         if np.isscalar(T_b):
             T_b = np.tile(T_b, sum(self.nSegments))
         Q_f = a_in @ np.atleast_1d(T_f_in) + a_b @ T_b
@@ -241,7 +278,8 @@ def get_fluid_heat_extraction_rate(
         return Q_f
 
     def get_network_inlet_temperature(
-            self, Q_t, T_b, m_flow_network, cp_f, nSegments):
+            self, Q_t, T_b, m_flow_network, cp_f, nSegments,
+            segment_ratios=None):
         """
         Returns the inlet fluid temperature of the network.
 
@@ -262,6 +300,12 @@ def get_network_inlet_temperature(
             Number of borehole segments for each borehole. If an int is
             supplied, all boreholes are considered to have the same number of
             segments.
+        segment_ratios :
+        (nSegments,) array or list of (nSegments[i],) arrays, 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
         -------
@@ -272,7 +316,7 @@ def get_network_inlet_temperature(
         """
         # Build coefficient matrices
         a_q, a_b = self.coefficients_network_inlet_temperature(
-                m_flow_network, cp_f, nSegments)
+                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))
@@ -282,7 +326,8 @@ def get_network_inlet_temperature(
         return T_f_in
 
     def get_network_outlet_temperature(
-            self, T_f_in, T_b, m_flow_network, cp_f, nSegments):
+            self, T_f_in, T_b, m_flow_network, cp_f, nSegments,
+            segment_ratios=None):
         """
         Returns the outlet fluid temperature of the network.
 
@@ -303,6 +348,12 @@ def get_network_outlet_temperature(
             Number of borehole segments for each borehole. If an int is
             supplied, all boreholes are considered to have the same number of
             segments.
+        segment_ratios :
+        (nSegments,) array or list of (nSegments[i],) arrays, 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
         -------
@@ -313,7 +364,7 @@ def get_network_outlet_temperature(
         """
         # Build coefficient matrices
         a_in, a_b = self.coefficients_network_outlet_temperature(
-                m_flow_network, cp_f, nSegments)
+                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))
@@ -323,7 +374,8 @@ def get_network_outlet_temperature(
         return T_f_out
 
     def get_network_heat_extraction_rate(
-            self, T_f_in, T_b, m_flow_network, cp_f, nSegments):
+            self, T_f_in, T_b, m_flow_network, cp_f, nSegments,
+            segment_ratios=None):
         """
         Returns the total heat extraction rate of the network.
 
@@ -344,6 +396,12 @@ def get_network_heat_extraction_rate(
             Number of borehole segments for each borehole. If an int is
             supplied, all boreholes are considered to have the same number of
             segments.
+        segment_ratios :
+        (nSegments,) array or list of (nSegments[i],) arrays, 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
         -------
@@ -353,7 +411,7 @@ def get_network_heat_extraction_rate(
 
         """
         a_in, a_b = self.coefficients_network_heat_extraction_rate(
-                m_flow_network, cp_f, nSegments)
+                m_flow_network, cp_f, nSegments, segment_ratios=segment_ratios)
         if np.isscalar(T_b):
             T_b = np.tile(T_b, sum(self.nSegments))
         Q_t = a_in @ np.atleast_1d(T_f_in) + a_b @ T_b
@@ -362,7 +420,8 @@ def get_network_heat_extraction_rate(
 
         return Q_t
 
-    def coefficients_inlet_temperature(self, m_flow_network, cp_f, nSegments):
+    def coefficients_inlet_temperature(
+            self, m_flow_network, cp_f, nSegments, segment_ratios=None):
         """
         Build coefficient matrices to evaluate intlet fluid temperatures of all
         boreholes.
@@ -387,6 +446,12 @@ def coefficients_inlet_temperature(self, m_flow_network, cp_f, nSegments):
             Number of borehole segments for each borehole. If an int is
             supplied, all boreholes are considered to have the same number of
             segments.
+        segment_ratios :
+        (nSegments,) array or list of (nSegments[i],) arrays, 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
         -------
@@ -400,22 +465,21 @@ def coefficients_inlet_temperature(self, m_flow_network, cp_f, nSegments):
         method_id = 0
         # Check if stored coefficients are available
         if self._check_coefficients(
-                m_flow_network, cp_f, nSegments, method_id):
+                m_flow_network, cp_f, nSegments, segment_ratios, method_id):
             a_in, a_b = self._get_stored_coefficients(method_id)
         else:
             # Update input variables
-            self._format_inputs(m_flow_network, cp_f, nSegments)
+            self._format_inputs(
+                m_flow_network, cp_f, nSegments, segment_ratios)
             # Coefficient matrices for borehole inlet temperatures:
             # [T_{f,b,in}] = [c_in]*[T_{f,n,in}] + [c_out]*[T_{f,b,out}]
             c_in = self._c_in
             c_out = self._c_out
             # Coefficient matrices for borehole outlet temperatures:
             # [T_{f,b,out}] = [A]*[T_{f,b,in}] + [B]*[T_{b}]
-            AB = list(zip(*[
-                self.p[i].coefficients_outlet_temperature(
-                    self._m_flow_borehole[i],
-                    self._cp_borehole[i],
-                    self.nSegments[i])
+            AB = list(zip(*[self.p[i].coefficients_outlet_temperature(
+                self._m_flow_borehole[i], self._cp_borehole[i],
+                self.nSegments[i], segment_ratios=self._segment_ratios[i])
                 for i in range(self.nBoreholes)]))
             A = block_diag(*AB[0])
             B = block_diag(*AB[1])
@@ -427,11 +491,13 @@ def coefficients_inlet_temperature(self, m_flow_network, cp_f, nSegments):
 
             # Store coefficients
             self._set_stored_coefficients(
-                m_flow_network, cp_f, nSegments, (a_in, a_b), method_id)
+                m_flow_network, cp_f, nSegments, segment_ratios, (a_in, a_b),
+                method_id)
 
         return a_in, a_b
 
-    def coefficients_outlet_temperature(self, m_flow_network, cp_f, nSegments):
+    def coefficients_outlet_temperature(
+            self, m_flow_network, cp_f, nSegments, segment_ratios=None):
         """
         Build coefficient matrices to evaluate outlet fluid temperatures of all
         boreholes.
@@ -456,6 +522,12 @@ def coefficients_outlet_temperature(self, m_flow_network, cp_f, nSegments):
             Number of borehole segments for each borehole. If an int is
             supplied, all boreholes are considered to have the same number of
             segments.
+        segment_ratios :
+        (nSegments,) array or list of (nSegments[i],) arrays, 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
         -------
@@ -469,11 +541,12 @@ def coefficients_outlet_temperature(self, m_flow_network, cp_f, nSegments):
         method_id = 1
         # Check if stored coefficients are available
         if self._check_coefficients(
-                m_flow_network, cp_f, nSegments, method_id):
+                m_flow_network, cp_f, nSegments, segment_ratios, method_id):
             a_in, a_b = self._get_stored_coefficients(method_id)
         else:
             # Update input variables
-            self._format_inputs(m_flow_network, cp_f, nSegments)
+            self._format_inputs(
+                m_flow_network, cp_f, nSegments, segment_ratios)
             # Coefficient matrices for borehole inlet temperatures:
             # [T_{f,b,in}] = [c_in]*[T_{f,n,in}] + [c_out]*[T_{f,b,out}]
             c_in = self._c_in
@@ -484,7 +557,8 @@ def coefficients_outlet_temperature(self, m_flow_network, cp_f, nSegments):
                 self.p[i].coefficients_outlet_temperature(
                     self._m_flow_borehole[i],
                     self._cp_borehole[i],
-                    self.nSegments[i])
+                    self.nSegments[i],
+                    segment_ratios=self._segment_ratios[i])
                 for i in range(self.nBoreholes)]))
             A = block_diag(*AB[0])
             B = block_diag(*AB[1])
@@ -496,12 +570,13 @@ def coefficients_outlet_temperature(self, m_flow_network, cp_f, nSegments):
 
             # Store coefficients
             self._set_stored_coefficients(
-                m_flow_network, cp_f, nSegments, (a_in, a_b), method_id)
+                m_flow_network, cp_f, nSegments, segment_ratios, (a_in, a_b),
+                method_id)
 
         return a_in, a_b
 
     def coefficients_network_inlet_temperature(
-            self, m_flow_network, cp_f, nSegments):
+            self, m_flow_network, cp_f, nSegments, segment_ratios=None):
         """
         Build coefficient matrices to evaluate intlet fluid temperature of the
         network.
@@ -526,6 +601,12 @@ def coefficients_network_inlet_temperature(
             Number of borehole segments for each borehole. If an int is
             supplied, all boreholes are considered to have the same number of
             segments.
+        segment_ratios :
+        (nSegments,) array or list of (nSegments[i],) arrays, 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
         -------
@@ -538,13 +619,15 @@ def coefficients_network_inlet_temperature(
         # method_id for coefficients_network_inlet_temperature is 2
         method_id = 2
         # Check if stored coefficients are available
-        if self._check_coefficients(m_flow_network, cp_f, nSegments, method_id):
+        if self._check_coefficients(
+                m_flow_network, cp_f, nSegments, segment_ratios, method_id):
             a_qf, a_b = self._get_stored_coefficients(method_id)
         else:
             # Coefficient matrices for network heat extraction rates:
             # [Q_{tot}] = [b_in]*[T_{f,n,in}] + [b_b]*[T_{b}]
             b_in, b_b = self.coefficients_network_heat_extraction_rate(
-                    m_flow_network, cp_f, nSegments)
+                    m_flow_network, cp_f, nSegments,
+                    segment_ratios=segment_ratios)
             # Coefficient matrices for network inlet temperature:
             # [T_{f,n,in}] = [a_qf]*[Q_{tot}] + [a_b]*[T_{b}]
             b_in_inv = np.linalg.inv(b_in)
@@ -553,12 +636,13 @@ def coefficients_network_inlet_temperature(
 
             # Store coefficients
             self._set_stored_coefficients(
-                m_flow_network, cp_f, nSegments, (a_qf, a_b), method_id)
+                m_flow_network, cp_f, nSegments, segment_ratios, (a_qf, a_b),
+                method_id)
 
         return a_qf, a_b
 
     def coefficients_network_outlet_temperature(
-            self, m_flow_network, cp_f, nSegments):
+            self, m_flow_network, cp_f, nSegments, segment_ratios=None):
         """
         Build coefficient matrices to evaluate outlet fluid temperature of the
         network.
@@ -583,6 +667,12 @@ def coefficients_network_outlet_temperature(
             Number of borehole segments for each borehole. If an int is
             supplied, all boreholes are considered to have the same number of
             segments.
+        segment_ratios :
+        (nSegments,) array or list of (nSegments[i],) arrays, 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
         -------
@@ -596,13 +686,13 @@ def coefficients_network_outlet_temperature(
         method_id = 3
         # Check if stored coefficients are available
         if self._check_coefficients(
-                m_flow_network, cp_f, nSegments, method_id):
+                m_flow_network, cp_f, nSegments, segment_ratios, method_id):
             a_in, a_b = self._get_stored_coefficients(method_id)
         else:
             # Coefficient matrices for borehole outlet temperatures:
             # [T_{f,b,out}] = [b_in]*[T_{f,n,in}] + [b_b]*[T_{b}]
             b_in, b_b = self.coefficients_outlet_temperature(
-                m_flow_network, cp_f, nSegments)
+                m_flow_network, cp_f, nSegments, segment_ratios=segment_ratios)
             # Coefficient matrices for network outlet temperature:
             # [T_{f,n,out}] = [a_in]*[T_{f,n,in}] + [a_b]*[T_{b}]
             mix_out = self._coefficients_mixing(m_flow_network)
@@ -611,12 +701,13 @@ def coefficients_network_outlet_temperature(
 
             # Store coefficients
             self._set_stored_coefficients(
-                m_flow_network, cp_f, nSegments, (a_in, a_b), method_id)
+                m_flow_network, cp_f, nSegments, segment_ratios, (a_in, a_b),
+                method_id)
 
         return a_in, a_b
 
     def coefficients_borehole_heat_extraction_rate(
-            self, m_flow_network, cp_f, nSegments):
+            self, m_flow_network, cp_f, nSegments, segment_ratios=None):
         """
         Build coefficient matrices to evaluate heat extraction rates of all
         boreholes segments.
@@ -641,6 +732,12 @@ def coefficients_borehole_heat_extraction_rate(
             Number of borehole segments for each borehole. If an int is
             supplied, all boreholes are considered to have the same number of
             segments.
+        segment_ratios :
+        (nSegments,) array or list of (nSegments[i],) arrays, 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
         -------
@@ -654,22 +751,24 @@ def coefficients_borehole_heat_extraction_rate(
         method_id = 4
         # Check if stored coefficients are available
         if self._check_coefficients(
-                m_flow_network, cp_f, nSegments, method_id):
+                m_flow_network, cp_f, nSegments, segment_ratios, method_id):
             a_in, a_b = self._get_stored_coefficients(method_id)
         else:
             # Update input variables
-            self._format_inputs(m_flow_network, cp_f, nSegments)
+            self._format_inputs(
+                m_flow_network, cp_f, nSegments, segment_ratios)
             # Coefficient matrices for borehole inlet temperatures:
             # [T_{f,b,in}] = [b_in]*[T_{f,n,in}] + [b_b]*[T_{b}]
             b_in, b_b = self.coefficients_inlet_temperature(
-                m_flow_network, cp_f, nSegments)
+                m_flow_network, cp_f, nSegments, segment_ratios=segment_ratios)
             # Coefficient matrices for borehole heat extraction rates:
             # [Q_{b}] = [A]*[T_{f,b,in}] + [B]*[T_{b}]
             AB = list(zip(*[
                 self.p[i].coefficients_borehole_heat_extraction_rate(
                     self._m_flow_borehole[i],
                     self._cp_borehole[i],
-                    self.nSegments[i])
+                    self.nSegments[i],
+                    segment_ratios=self._segment_ratios[i])
                 for i in range(self.nBoreholes)]))
             A = block_diag(*AB[0])
             B = block_diag(*AB[1])
@@ -680,12 +779,13 @@ def coefficients_borehole_heat_extraction_rate(
 
             # Store coefficients
             self._set_stored_coefficients(
-                m_flow_network, cp_f, nSegments, (a_in, a_b), method_id)
+                m_flow_network, cp_f, nSegments, segment_ratios, (a_in, a_b),
+                method_id)
 
         return a_in, a_b
 
     def coefficients_fluid_heat_extraction_rate(
-            self, m_flow_network, cp_f, nSegments):
+            self, m_flow_network, cp_f, nSegments, segment_ratios=None):
         """
         Build coefficient matrices to evaluate heat extraction rates of all
         boreholes.
@@ -710,6 +810,12 @@ def coefficients_fluid_heat_extraction_rate(
             Number of borehole segments for each borehole. If an int is
             supplied, all boreholes are considered to have the same number of
             segments.
+        segment_ratios :
+        (nSegments,) array or list of (nSegments[i],) arrays, 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
         -------
@@ -723,22 +829,24 @@ def coefficients_fluid_heat_extraction_rate(
         method_id = 5
         # Check if stored coefficients are available
         if self._check_coefficients(
-                m_flow_network, cp_f, nSegments, method_id):
+                m_flow_network, cp_f, nSegments, segment_ratios, method_id):
             a_in, a_b = self._get_stored_coefficients(method_id)
         else:
             # Update input variables
-            self._format_inputs(m_flow_network, cp_f, nSegments)
+            self._format_inputs(
+                m_flow_network, cp_f, nSegments, segment_ratios)
             # Coefficient matrices for borehole inlet temperatures:
             # [T_{f,b,in}] = [b_in]*[T_{f,n,in}] + [b_b]*[T_{b}]
             b_in, b_b = self.coefficients_inlet_temperature(
-                m_flow_network, cp_f, nSegments)
+                m_flow_network, cp_f, nSegments, segment_ratios=segment_ratios)
             # Coefficient matrices for fluid heat extraction rates:
             # [Q_{f}] = [A]*[T_{f,b,in}] + [B]*[T_{b}]
             AB = list(zip(*[
                 self.p[i].coefficients_fluid_heat_extraction_rate(
                     self._m_flow_borehole[i],
                     self._cp_borehole[i],
-                    self.nSegments[i])
+                    self.nSegments[i],
+                    segment_ratios=self._segment_ratios[i])
                 for i in range(self.nBoreholes)]))
             A = block_diag(*AB[0])
             B = block_diag(*AB[1])
@@ -749,12 +857,13 @@ def coefficients_fluid_heat_extraction_rate(
 
             # Store coefficients
             self._set_stored_coefficients(
-                m_flow_network, cp_f, nSegments, (a_in, a_b), method_id)
+                m_flow_network, cp_f, nSegments, segment_ratios, (a_in, a_b),
+                method_id)
 
         return a_in, a_b
 
     def coefficients_network_heat_extraction_rate(
-            self, m_flow_network, cp_f, nSegments):
+            self, m_flow_network, cp_f, nSegments, segment_ratios=None):
         """
         Build coefficient matrices to evaluate total heat extraction rate of
         the network.
@@ -779,6 +888,12 @@ def coefficients_network_heat_extraction_rate(
             Number of borehole segments for each borehole. If an int is
             supplied, all boreholes are considered to have the same number of
             segments.
+        segment_ratios :
+        (nSegments,) array or list of (nSegments[i],) arrays, 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
         -------
@@ -792,13 +907,14 @@ def coefficients_network_heat_extraction_rate(
         method_id = 6
         # Check if stored coefficients are available
         if self._check_coefficients(
-                m_flow_network, cp_f, nSegments, method_id):
+                m_flow_network, cp_f, nSegments, segment_ratios, method_id):
             a_in, a_b = self._get_stored_coefficients(method_id)
         else:
             # Coefficient matrices for fluid heat extraction rates:
             # [Q_{f}] = [b_in]*[T_{f,n,in}] + [b_b]*[T_{b}]
             b_in, b_b = self.coefficients_fluid_heat_extraction_rate(
-                    m_flow_network, cp_f, nSegments)
+                    m_flow_network, cp_f, nSegments,
+                    segment_ratios=segment_ratios)
             # The total network heat extraction rate is the sum of heat
             # extraction rates from all boreholes:
             # [Q_{tot}] = [a_in]*[T_{f,n,in}] + [a_b]*[T_{b}]
@@ -807,7 +923,8 @@ def coefficients_network_heat_extraction_rate(
 
             # Store coefficients
             self._set_stored_coefficients(
-                m_flow_network, cp_f, nSegments, (a_in, a_b), method_id)
+                m_flow_network, cp_f, nSegments, segment_ratios, (a_in, a_b),
+                method_id)
 
         return a_in, a_b
 
@@ -861,14 +978,17 @@ def _initialize_coefficients_connectivity(self):
             
         return
 
-    def _initialize_stored_coefficients(self, m_flow_network, cp_f, nSegments):
+    def _initialize_stored_coefficients(
+            self, m_flow_network, cp_f, nSegments, segment_ratios):
         nMethods = 7    # Number of class methods
         self._stored_coefficients = [() for i in range(nMethods)]
         self._stored_m_flow_cp = [np.empty(self.nInlets)
                                   for i in range(nMethods)]
         self._stored_nSegments = [np.nan for i in range(nMethods)]
+        self._stored_segment_ratios = [np.nan for i in range(nMethods)]
         self._m_flow_cp_model_variables = np.empty(self.nInlets)
         self._nSegments_model_variables = np.nan
+        self._segment_ratios_model_variables = np.nan
         self._mixing_m_flow = np.empty(self.nInlets)
         self._mixing_m_flow[:] = np.nan
         self._mix_out = np.empty((1, self.nBoreholes))
@@ -878,27 +998,29 @@ def _initialize_stored_coefficients(self, m_flow_network, cp_f, nSegments):
         # matrix coefficients.
         if m_flow_network is not None and cp_f is not None and nSegments is not None:
             self.coefficients_inlet_temperature(
-                m_flow_network, cp_f, nSegments)
+                m_flow_network, cp_f, nSegments, segment_ratios=segment_ratios)
             self.coefficients_outlet_temperature(
-                m_flow_network, cp_f, nSegments)
+                m_flow_network, cp_f, nSegments, segment_ratios=segment_ratios)
             self.coefficients_network_inlet_temperature(
-                m_flow_network, cp_f, nSegments)
+                m_flow_network, cp_f, nSegments, segment_ratios=segment_ratios)
             self.coefficients_network_outlet_temperature(
-                m_flow_network, cp_f, nSegments)
+                m_flow_network, cp_f, nSegments, segment_ratios=segment_ratios)
             self.coefficients_borehole_heat_extraction_rate(
-                m_flow_network, cp_f, nSegments)
+                m_flow_network, cp_f, nSegments, segment_ratios=segment_ratios)
             self.coefficients_fluid_heat_extraction_rate(
-                m_flow_network, cp_f, nSegments)
+                m_flow_network, cp_f, nSegments, segment_ratios=segment_ratios)
             self.coefficients_network_heat_extraction_rate(
-                m_flow_network, cp_f, nSegments)
+                m_flow_network, cp_f, nSegments, segment_ratios=segment_ratios)
 
         return
 
-    def _set_stored_coefficients(self, m_flow_network, cp_f, nSegments, coefficients,
-                                 method_id):
+    def _set_stored_coefficients(
+            self, m_flow_network, cp_f, nSegments, segment_ratios,
+            coefficients, method_id):
         self._stored_coefficients[method_id] = coefficients
         self._stored_m_flow_cp[method_id] = m_flow_network*cp_f
         self._stored_nSegments[method_id] = nSegments
+        self._stored_segment_ratios[method_id] = segment_ratios
 
         return
 
@@ -916,18 +1038,28 @@ def _check_mixing_coefficients(self, m_flow_network, tol=1e-6):
 
         return check
 
-    def _check_coefficients(self, m_flow_network, cp_f, nSegments, method_id, tol=1e-6):
+    def _check_coefficients(
+            self, m_flow_network, cp_f, nSegments, segment_ratios, method_id,
+            tol=1e-6):
         stored_m_flow_cp = self._stored_m_flow_cp[method_id]
         stored_nSegments = self._stored_nSegments[method_id]
+        stored_segment_ratios = self._stored_segment_ratios[method_id]
+        if stored_segment_ratios is None and segment_ratios is None:
+            check_ratios = True
+        elif isinstance(stored_segment_ratios, list) and isinstance(segment_ratios, list):
+            check_ratios = (np.all([len(segment_ratios[i]) == len(stored_segment_ratios[i]) for i in range(len(segment_ratios))]) and
+                            np.all([np.all(np.abs(segment_ratios[i] - stored_segment_ratios[i]) < np.abs(stored_segment_ratios[i])*tol) for i in range(len(segment_ratios))]))
+        else:
+            check_ratios = False
         if (np.all(np.abs(m_flow_network*cp_f - stored_m_flow_cp) < np.abs(stored_m_flow_cp)*tol)
-            and nSegments == stored_nSegments):
+            and nSegments == stored_nSegments and check_ratios):
             check = True
         else:
             check = False
 
         return check
 
-    def _format_inputs(self, m_flow_network, cp_f, nSegments):
+    def _format_inputs(self, m_flow_network, cp_f, nSegments, segment_ratios):
         """
         Format mass flow rate and heat capacity inputs.
         """
@@ -971,6 +1103,17 @@ def _format_inputs(self, m_flow_network, cp_f, nSegments):
         else:
             self.nSegments = nSegments
 
+        # Format segment ratios
+        if segment_ratios is None:
+            self._segment_ratios = [None] * self.nBoreholes
+        elif isinstance(segment_ratios, np.ndarray):
+            self._segment_ratios = [segment_ratios] * self.nBoreholes
+        elif isinstance(segment_ratios, list):
+            self._segment_ratios = segment_ratios
+        else:
+            raise ValueError(
+                'Incorrect format of the segment ratios list.')
+
 
 class _EquivalentNetwork(Network):
     """
@@ -1003,6 +1146,12 @@ class _EquivalentNetwork(Network):
         Number of line segments used per borehole. This parameter is used to
         initialize the coefficients if it is provided.
         Default is None.
+    segment_ratios :
+    (nSegments,) array or list of (nSegments[i],) arrays, 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.
 
     Notes
     -----
@@ -1026,7 +1175,7 @@ class _EquivalentNetwork(Network):
 
     """
     def __init__(self, equivalentBoreholes, pipes, m_flow_network=None,
-                 cp_f=None, nSegments=None):
+                 cp_f=None, nSegments=None, segment_ratios=None):
         self.b = equivalentBoreholes
         self.H_tot = sum([b.H*b.nBoreholes for b in self.b])
         self.nBoreholes = len(equivalentBoreholes)
@@ -1052,11 +1201,12 @@ def __init__(self, equivalentBoreholes, pipes, m_flow_network=None,
 
         # Initialize stored_coefficients
         self._initialize_coefficients_connectivity()
-        self._initialize_stored_coefficients(m_flow_network, cp_f, nSegments)
+        self._initialize_stored_coefficients(
+            m_flow_network, cp_f, nSegments, segment_ratios)
         return
 
     def coefficients_network_heat_extraction_rate(
-            self, m_flow_network, cp_f, nSegments):
+            self, m_flow_network, cp_f, nSegments, segment_ratios=None):
         """
         Build coefficient matrices to evaluate total heat extraction rate of
         the network.
@@ -1081,6 +1231,12 @@ def coefficients_network_heat_extraction_rate(
             Number of borehole segments for each borehole. If an int is
             supplied, all boreholes are considered to have the same number of
             segments.
+        segment_ratios :
+        (nSegments,) array or list of (nSegments[i],) arrays, 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
         -------
@@ -1094,13 +1250,14 @@ def coefficients_network_heat_extraction_rate(
         method_id = 6
         # Check if stored coefficients are available
         if self._check_coefficients(
-                m_flow_network, cp_f, nSegments, method_id):
+                m_flow_network, cp_f, nSegments, segment_ratios, method_id):
             a_in, a_b = self._get_stored_coefficients(method_id)
         else:
             # Coefficient matrices for fluid heat extraction rates:
             # [Q_{f}] = [b_in]*[T_{f,n,in}] + [b_b]*[T_{b}]
             b_in, b_b = self.coefficients_fluid_heat_extraction_rate(
-                    m_flow_network, cp_f, nSegments)
+                    m_flow_network, cp_f, nSegments,
+                    segment_ratios=segment_ratios)
             # The total network heat extraction rate is the sum of heat
             # extraction rates from all boreholes:
             # [Q_{tot}] = [a_in]*[T_{f,n,in}] + [a_b]*[T_{b}]
@@ -1109,7 +1266,8 @@ def coefficients_network_heat_extraction_rate(
 
             # Store coefficients
             self._set_stored_coefficients(
-                m_flow_network, cp_f, nSegments, (a_in, a_b), method_id)
+                m_flow_network, cp_f, nSegments, segment_ratios, (a_in, a_b),
+                method_id)
 
         return a_in, a_b
 
@@ -1144,7 +1302,7 @@ def _coefficients_mixing(self, m_flow_network):
             self._mixing_m_flow = m_flow_network
         return self._mix_out
 
-    def _format_inputs(self, m_flow_network, cp_f, nSegments):
+    def _format_inputs(self, m_flow_network, cp_f, nSegments, segment_ratios):
         """
         Format mass flow rate and heat capacity inputs.
         """
@@ -1181,14 +1339,25 @@ def _format_inputs(self, m_flow_network, cp_f, nSegments):
 
         # Format number of segments for each borehole
         nSeg = np.atleast_1d(nSegments)
-        if len(nSeg) == 1:
-            self.nSegments = [nSegments for i in range(self.nBoreholes)]
+        if len(nSeg) == 1 and not isinstance(nSegments, list):
+            self.nSegments = [nSegments] * self.nBoreholes
         elif not len(nSeg) == self.nBoreholes:
             raise ValueError(
                 'Incorrect length of number of segments list.')
         else:
             self.nSegments = nSegments
 
+        # Format segment ratios
+        if segment_ratios is None:
+            self._segment_ratios = [None] * self.nBoreholes
+        elif isinstance(segment_ratios, np.ndarray):
+            self._segment_ratios = [segment_ratios] * self.nBoreholes
+        elif isinstance(segment_ratios, list):
+            self._segment_ratios = segment_ratios
+        else:
+            raise ValueError(
+                'Incorrect format of the segment ratios list.')
+
 
 def network_thermal_resistance(network, m_flow_network, cp_f):
     """
diff --git a/pygfunction/pipes.py b/pygfunction/pipes.py
index e1faa709..50bbf3ec 100644
--- a/pygfunction/pipes.py
+++ b/pygfunction/pipes.py
@@ -42,7 +42,8 @@ def __init__(self, borehole):
         self.nInlets = 1
         self.nOutlets = 1
 
-    def get_temperature(self, z, T_f_in, T_b, m_flow_borehole, cp_f):
+    def get_temperature(
+            self, z, T_f_in, T_b, m_flow_borehole, cp_f, segment_ratios=None):
         """
         Returns the fluid temperatures of the borehole at a depth (z).
 
@@ -58,6 +59,11 @@ def get_temperature(self, z, T_f_in, T_b, m_flow_borehole, cp_f):
             Inlet mass flow rates (in kg/s) into the borehole.
         cp_f : float or (nInlets,) array
             Fluid specific isobaric heat capacity (in J/kg.degC).
+        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
         -------
@@ -69,8 +75,10 @@ def get_temperature(self, z, T_f_in, T_b, m_flow_borehole, cp_f):
         T_b = np.atleast_1d(T_b)
         nSegments = len(T_b)
         z_all = np.atleast_1d(z).flatten()
-        AB = list(zip(*[self.coefficients_temperature(
-            zi, m_flow_borehole, cp_f, nSegments) for zi in z_all]))
+        AB = list(zip(*[
+            self.coefficients_temperature(
+                zi, m_flow_borehole, cp_f, nSegments, segment_ratios=segment_ratios)
+            for zi in z_all]))
         a_in = np.stack(AB[0], axis=-1)
         a_b = np.stack(AB[1], axis=-1)
         T_f = np.einsum('ijk,j->ki', a_in, np.atleast_1d(T_f_in)) \
@@ -81,7 +89,8 @@ def get_temperature(self, z, T_f_in, T_b, m_flow_borehole, cp_f):
             T_f = T_f.flatten()
         return T_f
 
-    def get_inlet_temperature(self, Q_f, T_b, m_flow_borehole, cp_f):
+    def get_inlet_temperature(
+            self, Q_f, T_b, m_flow_borehole, cp_f, segment_ratios=None):
         """
         Returns the outlet fluid temperatures of the borehole.
 
@@ -95,6 +104,11 @@ def get_inlet_temperature(self, Q_f, T_b, m_flow_borehole, cp_f):
             Inlet mass flow rates (in kg/s) into the borehole.
         cp_f : float or (nInlets,) array
             Fluid specific isobaric heat capacity (in J/kg.degC).
+        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
         -------
@@ -107,7 +121,7 @@ def get_inlet_temperature(self, Q_f, T_b, m_flow_borehole, cp_f):
         nSegments = len(T_b)
         # Build coefficient matrices
         a_qf, a_b = self.coefficients_inlet_temperature(
-            m_flow_borehole, cp_f, nSegments)
+            m_flow_borehole, cp_f, nSegments, segment_ratios=segment_ratios)
         # Evaluate outlet temperatures
         T_f_in = a_qf @ np.atleast_1d(Q_f) + a_b @ T_b
         # Return float if Qf was supplied as scalar
@@ -115,7 +129,8 @@ def get_inlet_temperature(self, Q_f, T_b, m_flow_borehole, cp_f):
             T_f_in = T_f_in.item()
         return T_f_in
 
-    def get_outlet_temperature(self, T_f_in, T_b, m_flow_borehole, cp_f):
+    def get_outlet_temperature(
+            self, T_f_in, T_b, m_flow_borehole, cp_f, segment_ratios=None):
         """
         Returns the outlet fluid temperatures of the borehole.
 
@@ -129,6 +144,11 @@ def get_outlet_temperature(self, T_f_in, T_b, m_flow_borehole, cp_f):
             Inlet mass flow rates (in kg/s) into the borehole.
         cp_f : float or (nInlets,) array
             Fluid specific isobaric heat capacity (in J/kg.degC).
+        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
         -------
@@ -141,7 +161,7 @@ def get_outlet_temperature(self, T_f_in, T_b, m_flow_borehole, cp_f):
         nSegments = len(T_b)
         # Build coefficient matrices
         a_in, a_b = self.coefficients_outlet_temperature(
-            m_flow_borehole, cp_f, nSegments)
+            m_flow_borehole, cp_f, nSegments, segment_ratios=segment_ratios)
         # Evaluate outlet temperatures
         T_f_out = a_in @ np.atleast_1d(T_f_in) + a_b @ T_b
         # Return float if Tin was supplied as scalar
@@ -150,7 +170,7 @@ def get_outlet_temperature(self, T_f_in, T_b, m_flow_borehole, cp_f):
         return T_f_out
 
     def get_borehole_heat_extraction_rate(
-            self, T_f_in, T_b, m_flow_borehole, cp_f):
+            self, T_f_in, T_b, m_flow_borehole, cp_f, segment_ratios=None):
         """
         Returns the heat extraction rates of the borehole.
 
@@ -164,6 +184,11 @@ def get_borehole_heat_extraction_rate(
             Inlet mass flow rates (in kg/s) into the borehole.
         cp_f : float or (nInlets,) array
             Fluid specific isobaric heat capacity (in J/kg.degC).
+        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
         -------
@@ -175,7 +200,7 @@ def get_borehole_heat_extraction_rate(
         T_b = np.atleast_1d(T_b)
         nSegments = len(T_b)
         a_in, a_b = self.coefficients_borehole_heat_extraction_rate(
-            m_flow_borehole, cp_f, nSegments)
+            m_flow_borehole, cp_f, nSegments, segment_ratios=segment_ratios)
         Q_b = a_in @ np.atleast_1d(T_f_in) + a_b @ T_b
         # Return float if Tb was supplied as scalar
         if np.isscalar(T_b) and not np.isscalar(Q_b):
@@ -183,7 +208,7 @@ def get_borehole_heat_extraction_rate(
         return Q_b
 
     def get_fluid_heat_extraction_rate(
-            self, T_f_in, T_b, m_flow_borehole, cp_f):
+            self, T_f_in, T_b, m_flow_borehole, cp_f, segment_ratios=None):
         """
         Returns the heat extraction rates of the borehole.
 
@@ -197,6 +222,11 @@ def get_fluid_heat_extraction_rate(
             Inlet mass flow rates (in kg/s) into the borehole.
         cp_f : float or (nInlets,) array
             Fluid specific isobaric heat capacity (in J/kg.degC).
+        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
         -------
@@ -208,7 +238,7 @@ def get_fluid_heat_extraction_rate(
         T_b = np.atleast_1d(T_b)
         nSegments = len(T_b)
         a_in, a_b = self.coefficients_fluid_heat_extraction_rate(
-            m_flow_borehole, cp_f, nSegments)
+            m_flow_borehole, cp_f, nSegments, segment_ratios=segment_ratios)
         Q_f = a_in @ np.atleast_1d(T_f_in) + a_b @ T_b
         # Return float if Tb was supplied as scalar
         if np.isscalar(T_f_in) and not np.isscalar(Q_f):
@@ -216,7 +246,7 @@ def get_fluid_heat_extraction_rate(
         return Q_f
 
     def get_total_heat_extraction_rate(
-            self, T_f_in, T_b, m_flow_borehole, cp_f):
+            self, T_f_in, T_b, m_flow_borehole, cp_f, segment_ratios=None):
         """
         Returns the total heat extraction rate of the borehole.
 
@@ -230,6 +260,11 @@ def get_total_heat_extraction_rate(
             Inlet mass flow rates (in kg/s) into the borehole.
         cp_f : float or (nInlets,) array
             Fluid specific isobaric heat capacity (in J/kg.degC).
+        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
         -------
@@ -238,11 +273,12 @@ def get_total_heat_extraction_rate(
 
         """
         Q_f = self.get_fluid_heat_extraction_rate(
-            T_f_in, T_b, m_flow_borehole, cp_f)
+            T_f_in, T_b, m_flow_borehole, cp_f, segment_ratios=segment_ratios)
         Q_t = np.sum(Q_f)
         return Q_t
 
-    def coefficients_inlet_temperature(self, m_flow_borehole, cp_f, nSegments):
+    def coefficients_inlet_temperature(
+            self, m_flow_borehole, cp_f, nSegments, segment_ratios=None):
         """
         Build coefficient matrices to evaluate outlet fluid temperature.
 
@@ -261,6 +297,11 @@ def coefficients_inlet_temperature(self, m_flow_borehole, cp_f, nSegments):
             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
         -------
@@ -274,13 +315,14 @@ def coefficients_inlet_temperature(self, m_flow_borehole, cp_f, nSegments):
         method_id = 3
         # Check if stored coefficients are available
         if self._check_coefficients(
-                m_flow_borehole, cp_f, nSegments, method_id):
+                m_flow_borehole, cp_f, nSegments, segment_ratios, method_id):
             a_qf, a_b = self._get_stored_coefficients(method_id)
         else:
             # Coefficient matrices for fluid heat extraction rates:
             # [Q_{f}] = [b_in]*[T_{f,in}] + [b_b]*[T_{b}]
             b_in, b_b = self.coefficients_fluid_heat_extraction_rate(
-                m_flow_borehole, cp_f, nSegments)
+                m_flow_borehole, cp_f, nSegments,
+                segment_ratios=segment_ratios)
             b_in_m1 = np.linalg.inv(b_in)
 
             # Matrices for fluid heat extraction rates:
@@ -290,12 +332,13 @@ def coefficients_inlet_temperature(self, m_flow_borehole, cp_f, nSegments):
 
             # Store coefficients
             self._set_stored_coefficients(
-                m_flow_borehole, cp_f, nSegments, (a_qf, a_b), method_id)
+                m_flow_borehole, cp_f, nSegments, segment_ratios, (a_qf, a_b),
+                method_id)
 
         return a_qf, a_b
 
     def coefficients_outlet_temperature(
-            self, m_flow_borehole, cp_f, nSegments):
+            self, m_flow_borehole, cp_f, nSegments, segment_ratios=None):
         """
         Build coefficient matrices to evaluate outlet fluid temperature.
 
@@ -314,6 +357,11 @@ def coefficients_outlet_temperature(
             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
         -------
@@ -327,22 +375,25 @@ def coefficients_outlet_temperature(
         method_id = 4
         # Check if stored coefficients are available
         if self._check_coefficients(
-                m_flow_borehole, cp_f, nSegments, method_id):
+                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
-            if self._check_coefficients(m_flow_borehole, cp_f, nSegments, 0):
+            if self._check_coefficients(
+                    m_flow_borehole, cp_f, nSegments, segment_ratios, 0):
                 b_in, b_out, b_b = self._get_stored_coefficients(0)
             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(
-                    m_flow_borehole, cp_f, nSegments)
+                    m_flow_borehole, cp_f, nSegments,
+                    segment_ratios=segment_ratios)
 
                 # Store coefficients
                 self._set_stored_coefficients(
-                    m_flow_borehole, cp_f, nSegments, (b_in, b_out, b_b), 0)
+                    m_flow_borehole, cp_f, nSegments, segment_ratios,
+                    (b_in, b_out, b_b), 0)
 
             # Final coefficient matrices for outlet temperatures:
             # [T_{f,out}] = [a_in]*[T_{f,in}] + [a_b]*[T_b]
@@ -352,11 +403,13 @@ def coefficients_outlet_temperature(
 
             # Store coefficients
             self._set_stored_coefficients(
-                m_flow_borehole, cp_f, nSegments, (a_in, a_b), method_id)
+                m_flow_borehole, cp_f, nSegments, segment_ratios, (a_in, a_b),
+                method_id)
 
         return a_in, a_b
 
-    def coefficients_temperature(self, z, m_flow_borehole, cp_f, nSegments):
+    def coefficients_temperature(
+            self, z, m_flow_borehole, cp_f, nSegments, segment_ratios=None):
         """
         Build coefficient matrices to evaluate fluid temperatures at a depth
         (z).
@@ -378,6 +431,11 @@ def coefficients_temperature(self, z, m_flow_borehole, cp_f, nSegments):
             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
         -------
@@ -393,26 +451,29 @@ def coefficients_temperature(self, z, m_flow_borehole, cp_f, nSegments):
         # Coefficient matrices for outlet temperatures:
         # [T_{f,out}] = [b_in]*[T_{f,in}] + [b_b]*[T_b]
         b_in, b_b = self.coefficients_outlet_temperature(
-            m_flow_borehole, cp_f, nSegments)
+            m_flow_borehole, cp_f, nSegments, segment_ratios=segment_ratios)
 
         # 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, 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)
         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)
+                m_flow_borehole, cp_f, nSegments, segment_ratios=segment_ratios)
 
             # Store coefficients
             self._set_stored_coefficients(
-                m_flow_borehole, cp_f, nSegments, (c_in, c_out, c_b), 1)
+                m_flow_borehole, cp_f, nSegments, segment_ratios,
+                (c_in, c_out, c_b), 1)
 
         # Coefficient matrices from general solution:
         # [T_f](z) = [d_f0]*[T_f](0) + [d_b]*[T_b]
-        d_f0, d_b = self._general_solution(z, m_flow_borehole, cp_f, nSegments)
+        d_f0, d_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]
@@ -422,7 +483,7 @@ def coefficients_temperature(self, z, m_flow_borehole, cp_f, nSegments):
         return a_in, a_b
 
     def coefficients_borehole_heat_extraction_rate(
-            self, m_flow_borehole, cp_f, nSegments):
+            self, m_flow_borehole, cp_f, nSegments, segment_ratios=None):
         """
         Build coefficient matrices to evaluate heat extraction rates.
 
@@ -441,6 +502,11 @@ def coefficients_borehole_heat_extraction_rate(
             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
         -------
@@ -455,11 +521,13 @@ def coefficients_borehole_heat_extraction_rate(
 
         nPipes = self.nPipes
         # Check if stored coefficients are available
-        if self._check_coefficients(m_flow_borehole, cp_f, nSegments, method_id):
+        if self._check_coefficients(
+                m_flow_borehole, cp_f, nSegments, segment_ratios, method_id):
             a_in, a_b = self._get_stored_coefficients(method_id)
         else:
             # Update input variables
-            self._format_inputs(m_flow_borehole, cp_f, nSegments)
+            self._format_inputs(
+                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],
@@ -470,25 +538,29 @@ def coefficients_borehole_heat_extraction_rate(
             a_b = np.zeros((nSegments, nSegments))
             # Heat extraction rates are calculated from an energy balance on a
             # borehole segment.
-            z1 = 0.
+            z = self.b._segment_edges(nSegments, segment_ratios=segment_ratios)
+            z1 = z[0]
             aTf1, bTf1 = self.coefficients_temperature(
-                z1, m_flow_borehole, cp_f, nSegments)
+                z1, m_flow_borehole, cp_f, nSegments,
+                segment_ratios=segment_ratios)
             for i in range(nSegments):
-                z2 = (i + 1) * self.b.H / nSegments
+                z2 = z[i+1]
                 aTf2, bTf2 = self.coefficients_temperature(
-                    z2, m_flow_borehole, cp_f, nSegments)
+                    z2, m_flow_borehole, cp_f, nSegments,
+                    segment_ratios=segment_ratios)
                 a_in[i, :] = mcp @ (aTf1 - aTf2)
                 a_b[i, :] = mcp @ (bTf1 - bTf2)
                 aTf1, bTf1 = aTf2, bTf2
 
             # Store coefficients
             self._set_stored_coefficients(
-                m_flow_borehole, cp_f, nSegments, (a_in, a_b), method_id)
+                m_flow_borehole, cp_f, nSegments, segment_ratios,
+                (a_in, a_b), method_id)
 
         return a_in, a_b
 
     def coefficients_fluid_heat_extraction_rate(
-            self, m_flow_borehole, cp_f, nSegments):
+            self, m_flow_borehole, cp_f, nSegments, segment_ratios=None):
         """
         Build coefficient matrices to evaluate heat extraction rates.
 
@@ -507,6 +579,11 @@ def coefficients_fluid_heat_extraction_rate(
             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
         -------
@@ -519,16 +596,19 @@ def coefficients_fluid_heat_extraction_rate(
         # method_id for coefficients_fluid_heat_extraction_rate is 7
         method_id = 7
         # Check if stored coefficients are available
-        if self._check_coefficients(m_flow_borehole, cp_f, nSegments, method_id):
+        if self._check_coefficients(
+                m_flow_borehole, cp_f, nSegments, segment_ratios, method_id):
             a_in, a_b = self._get_stored_coefficients(method_id)
         else:
             # Update input variables
-            self._format_inputs(m_flow_borehole, cp_f, nSegments)
+            self._format_inputs(
+                m_flow_borehole, cp_f, nSegments, segment_ratios)
 
             # Coefficient matrices for outlet temperatures:
             # [T_{f,out}] = [b_in]*[T_{f,in}] + [b_b]*[T_b]
             b_in, b_b = self.coefficients_outlet_temperature(
-                m_flow_borehole, cp_f, nSegments)
+                m_flow_borehole, cp_f, nSegments,
+                segment_ratios=segment_ratios)
 
             # Intermediate matrices for fluid heat extraction rates:
             # [Q_{f}] = [c_in]*[T_{f,in}] + [c_out]*[T_{f,out}]
@@ -543,7 +623,8 @@ def coefficients_fluid_heat_extraction_rate(
 
             # Store coefficients
             self._set_stored_coefficients(
-                m_flow_borehole, cp_f, nSegments, (a_in, a_b), method_id)
+                m_flow_borehole, cp_f, nSegments, segment_ratios, (a_in, a_b),
+                method_id)
 
         return a_in, a_b
 
@@ -619,41 +700,66 @@ def _initialize_stored_coefficients(self):
         self._stored_m_flow_cp = [np.empty(self.nInlets)
                                   for i in range(nMethods)]
         self._stored_nSegments = [np.nan for i in range(nMethods)]
+        self._stored_segment_ratios = [np.nan for i in range(nMethods)]
         self._m_flow_cp_model_variables = np.empty(self.nInlets)
         self._nSegments_model_variables = np.nan
+        self._segment_ratios_model_variables = np.nan
 
 
     def _set_stored_coefficients(
-            self, m_flow_borehole, cp_f, nSegments, coefficients, method_id):
+            self, m_flow_borehole, cp_f, nSegments, segment_ratios, coefficients, method_id):
         self._stored_coefficients[method_id] = coefficients
         self._stored_m_flow_cp[method_id] = m_flow_borehole*cp_f
         self._stored_nSegments[method_id] = nSegments
+        self._stored_segment_ratios[method_id] = segment_ratios
 
     def _get_stored_coefficients(self, method_id):
         coefficients = self._stored_coefficients[method_id]
-
         return coefficients
 
-    def _check_model_variables(self, m_flow_borehole, cp_f, nSegments, tol=1e-6):
+    def _check_model_variables(
+            self, m_flow_borehole, cp_f, nSegments, segment_ratios, tol=1e-6):
         stored_m_flow_cp = self._m_flow_cp_model_variables
         stored_nSegments = self._nSegments_model_variables
+        stored_segment_ratios = self._segment_ratios_model_variables
+        if stored_segment_ratios is None and segment_ratios is None:
+            check_ratios = True
+        elif isinstance(stored_segment_ratios, np.ndarray) and isinstance(segment_ratios, np.ndarray):
+            check_ratios = (len(segment_ratios) == len(stored_segment_ratios) and
+                            np.all(np.abs(segment_ratios - stored_segment_ratios) < np.abs(stored_segment_ratios)*tol))
+        else:
+            check_ratios = False
         if (np.all(np.abs(m_flow_borehole*cp_f - stored_m_flow_cp) < np.abs(stored_m_flow_cp)*tol)
-            and nSegments == stored_nSegments):
+            and nSegments == stored_nSegments
+            and check_ratios):
             check = True
         else:
-            self._update_model_variables(m_flow_borehole, cp_f, nSegments)
+            self._update_model_variables(
+                m_flow_borehole, cp_f, nSegments, segment_ratios)
             self._m_flow_cp_model_variables = m_flow_borehole*cp_f
             self._nSegments_model_variables = nSegments
+            self._segment_ratios_model_variables = segment_ratios
             check = False
 
         return check
 
     def _check_coefficients(
-            self, m_flow_borehole, cp_f, nSegments, method_id, tol=1e-6):
+            self, m_flow_borehole, cp_f, nSegments, segment_ratios, method_id,
+            tol=1e-6):
         stored_m_flow_cp = self._stored_m_flow_cp[method_id]
         stored_nSegments = self._stored_nSegments[method_id]
-        if (np.all(np.abs(m_flow_borehole*cp_f - stored_m_flow_cp) < np.abs(stored_m_flow_cp)*tol)
-            and nSegments == stored_nSegments):
+        stored_segment_ratios = self._stored_segment_ratios[method_id]
+        if nSegments == stored_nSegments:
+            if stored_segment_ratios is None and segment_ratios is None:
+                check_ratios = True
+            elif isinstance(stored_segment_ratios, np.ndarray) and isinstance(segment_ratios, np.ndarray):
+                check_ratios = np.all(np.abs(segment_ratios - stored_segment_ratios) < np.abs(stored_segment_ratios)*tol)
+            else:
+                check_ratios = False
+        else:
+            check_ratios = False
+        if (check_ratios
+            and np.all(np.abs(m_flow_borehole*cp_f - stored_m_flow_cp) < np.abs(stored_m_flow_cp)*tol)):
             check = True
         else:
             check = False
@@ -721,7 +827,8 @@ def _check_geometry(self):
 
         return True
 
-    def _continuity_condition_base(self, m_flow_borehole, cp_f, nSegments):
+    def _continuity_condition_base(
+            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]
         """
@@ -730,7 +837,8 @@ def _continuity_condition_base(self, m_flow_borehole, cp_f, nSegments):
             'this method should return matrices for the relation: '
             '[a_out]*[T_{f,out}] = [a_in]*[T_{f,in}] + [a_b]*[T_b]')
 
-    def _continuity_condition_head(self, m_flow_borehole, cp_f, nSegments):
+    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]
         """
@@ -740,7 +848,8 @@ def _continuity_condition_head(self, m_flow_borehole, cp_f, nSegments):
             '[T_f](z=0) = [a_in]*[T_{f,in}] + [a_out]*[T_{f,out}] '
             '+ [a_b]*[T_b]')
 
-    def _general_solution(self, z, m_flow_borehole, cp_f, nSegments):
+    def _general_solution(
+            self, z, m_flow_borehole, cp_f, nSegments, segment_ratios=None):
         """ Returns coefficients for the relation
             [T_f](z) = [a_f0]*[T_f](0) + [a_b]*[T_b]
         """
@@ -749,7 +858,8 @@ def _general_solution(self, z, m_flow_borehole, cp_f, nSegments):
             'this method should return matrices for the relation: '
             '[T_f](z) = [a_f0]*[T_f](0) + [a_b]*[T_b]')
 
-    def _update_model_variables(self, m_flow_borehole, cp_f, nSegments):
+    def _update_model_variables(
+            self, m_flow_borehole, cp_f, nSegments, segment_ratios):
         """
         Evaluate common coefficients needed in other class methods.
         """
@@ -758,7 +868,8 @@ def _update_model_variables(self, m_flow_borehole, cp_f, nSegments):
             'this method should evaluate common coefficients needed in other '
             'class methods.')
 
-    def _format_inputs(self, m_flow_borehole, cp_f, nSegments):
+    def _format_inputs(
+            self, m_flow_borehole, cp_f, nSegments, segment_ratios):
         """
         Format arrays of mass flow rates and heat capacity.
         """
@@ -843,7 +954,8 @@ def __init__(self, pos, r_in, r_out, borehole, k_s, k_g, R_fp, J=2):
         # Initialize stored_coefficients
         self._initialize_stored_coefficients()
 
-    def _continuity_condition_base(self, m_flow_borehole, cp_f, nSegments):
+    def _continuity_condition_base(
+            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).
@@ -864,6 +976,11 @@ def _continuity_condition_base(self, m_flow_borehole, cp_f, nSegments):
             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
         -------
@@ -876,7 +993,8 @@ def _continuity_condition_base(self, m_flow_borehole, cp_f, nSegments):
 
         """
         # Check if model variables need to be updated
-        self._check_model_variables(m_flow_borehole, cp_f, nSegments)
+        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))
@@ -886,7 +1004,8 @@ def _continuity_condition_base(self, m_flow_borehole, cp_f, nSegments):
         a_out = np.array([[1.0]])
 
         a_b = np.zeros((self.nOutlets, nSegments))
-        z = (nSegments - np.arange(nSegments + 1)) * self.b.H / nSegments
+
+        z = self.b._segment_edges(nSegments, segment_ratios=segment_ratios)[::-1]
         F4 = self._F4(z)
         dF4 = F4[:-1] - F4[1:]
         F5 = self._F5(z)
@@ -895,7 +1014,8 @@ def _continuity_condition_base(self, m_flow_borehole, cp_f, nSegments):
 
         return a_in, a_out, a_b
 
-    def _continuity_condition_head(self, m_flow_borehole, cp_f, nSegments):
+    def _continuity_condition_head(
+            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
@@ -917,6 +1037,11 @@ def _continuity_condition_head(self, m_flow_borehole, cp_f, nSegments):
             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
         -------
@@ -929,7 +1054,8 @@ def _continuity_condition_head(self, m_flow_borehole, cp_f, nSegments):
 
         """
         # Check if model variables need to be updated
-        self._check_model_variables(m_flow_borehole, cp_f, nSegments)
+        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]])
@@ -938,7 +1064,8 @@ def _continuity_condition_head(self, m_flow_borehole, cp_f, nSegments):
 
         return a_in, a_out, a_b
 
-    def _general_solution(self, z, m_flow_borehole, cp_f, nSegments):
+    def _general_solution(
+            self, z, m_flow_borehole, cp_f, nSegments, segment_ratios=None):
         """
         General solution for fluid temperatures at a depth (z).
 
@@ -959,6 +1086,11 @@ def _general_solution(self, z, m_flow_borehole, cp_f, nSegments):
             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
         -------
@@ -969,15 +1101,17 @@ def _general_solution(self, z, m_flow_borehole, cp_f, nSegments):
 
         """
         # Check if model variables need to be updated
-        self._check_model_variables(m_flow_borehole, cp_f, nSegments)
+        self._check_model_variables(
+            m_flow_borehole, cp_f, nSegments, segment_ratios)
 
         a_f0 = np.array([[self._f1(z), self._f2(z)],
                         [-self._f2(z), self._f3(z)]])
 
         a_b = np.zeros((2*self.nPipes, nSegments))
-        N = int(np.ceil(z/self.b.H*nSegments))
-        z1 = z - np.minimum(np.arange(1, N+1)*self.b.H/nSegments, z)
-        z2 = z - np.arange(N) * self.b.H / nSegments
+        z_all = self.b._segment_edges(nSegments, segment_ratios=segment_ratios)
+        N = min(int(np.sum(z > z_all)), nSegments)
+        z1 = z - np.minimum(z_all[1:N+1], z)
+        z2 = z - z_all[:N]
         dF4 = self._F4(z2) - self._F4(z1)
         dF5 = self._F5(z2) - self._F5(z1)
         a_b[0, :N] = dF4
@@ -985,7 +1119,8 @@ def _general_solution(self, z, m_flow_borehole, cp_f, nSegments):
 
         return a_f0, a_b
 
-    def _update_model_variables(self, m_flow_borehole, cp_f, nSegments):
+    def _update_model_variables(
+            self, m_flow_borehole, cp_f, nSegments, segment_ratios):
         """
         Evaluate dimensionless resistances for Hellstrom (1991) solution.
 
@@ -997,10 +1132,14 @@ def _update_model_variables(self, m_flow_borehole, cp_f, nSegments):
             Fluid specific isobaric heat capacity (in J/kg.degC).
         nSegments : int
             Number of borehole segments.
+        segment_ratios : (nSegments,) array
+            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.
         """
 
         # Format mass flow rate and heat capacity inputs
-        self._format_inputs(m_flow_borehole, cp_f, nSegments)
+        self._format_inputs(m_flow_borehole, cp_f, nSegments, segment_ratios)
         m_flow_in = self._m_flow_in
         cp_in = self._cp_in
 
@@ -1015,7 +1154,7 @@ def _update_model_variables(self, m_flow_borehole, cp_f, nSegments):
         self._delta = 1./self._gamma \
             * (self._beta12 + 0.5*(self._beta1+self._beta2))
 
-    def _format_inputs(self, m_flow_borehole, cp_f, nSegments):
+    def _format_inputs(self, m_flow_borehole, cp_f, nSegments, segment_ratios):
         """
         Format mass flow rate and heat capacity inputs.
 
@@ -1027,6 +1166,10 @@ def _format_inputs(self, m_flow_borehole, cp_f, nSegments):
             Fluid specific isobaric heat capacity (in J/kg.degC).
         nSegments : int
             Number of borehole segments.
+        segment_ratios : (nSegments,) array
+            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.
         """
 
         # Format mass flow rate inputs
@@ -1241,7 +1384,8 @@ def __init__(self, pos, r_in, r_out, borehole, k_s,
         # Initialize stored_coefficients
         self._initialize_stored_coefficients()
 
-    def _continuity_condition_base(self, m_flow_borehole, cp_f, nSegments):
+    def _continuity_condition_base(
+            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).
@@ -1261,6 +1405,11 @@ def _continuity_condition_base(self, m_flow_borehole, cp_f, nSegments):
             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
         -------
@@ -1273,12 +1422,13 @@ def _continuity_condition_base(self, m_flow_borehole, cp_f, nSegments):
 
         """
         # Check if model variables need to be updated
-        self._check_model_variables(m_flow_borehole, cp_f, nSegments)
+        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)
+            m_flow_borehole, cp_f, nSegments, segment_ratios=segment_ratios)
         b_u_m1 = np.linalg.inv(b_u)
 
         if self.config == 'parallel':
@@ -1324,7 +1474,8 @@ def _continuity_condition_base(self, m_flow_borehole, cp_f, nSegments):
 
         return a_in, a_out, a_b
 
-    def _continuity_condition_head(self, m_flow_borehole, cp_f, nSegments):
+    def _continuity_condition_head(
+            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
@@ -1346,6 +1497,11 @@ def _continuity_condition_head(self, m_flow_borehole, cp_f, nSegments):
             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
         -------
@@ -1358,7 +1514,8 @@ def _continuity_condition_head(self, m_flow_borehole, cp_f, nSegments):
 
         """
         # Check if model variables need to be updated
-        self._check_model_variables(m_flow_borehole, cp_f, nSegments)
+        self._check_model_variables(
+            m_flow_borehole, cp_f, nSegments, segment_ratios)
 
         if self.config == 'parallel':
             a_in = np.vstack((np.ones((self.nPipes, self.nInlets)),
@@ -1371,7 +1528,8 @@ def _continuity_condition_head(self, m_flow_borehole, cp_f, nSegments):
             # 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)
+                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)
@@ -1400,7 +1558,8 @@ def _continuity_condition_head(self, m_flow_borehole, cp_f, nSegments):
 
         return a_in, a_out, a_b
 
-    def _continuity_condition(self, m_flow_borehole, cp_f, nSegments):
+    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).
@@ -1421,6 +1580,11 @@ def _continuity_condition(self, m_flow_borehole, cp_f, nSegments):
             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
         -------
@@ -1448,10 +1612,11 @@ def _continuity_condition(self, m_flow_borehole, cp_f, nSegments):
         # 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 = H - np.arange(nSegments + 1) * H / 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),
@@ -1463,7 +1628,8 @@ def _continuity_condition(self, m_flow_borehole, cp_f, nSegments):
 
         return a_d, a_u, a_b
 
-    def _general_solution(self, z, m_flow_borehole, cp_f, nSegments):
+    def _general_solution(
+            self, z, m_flow_borehole, cp_f, nSegments, segment_ratios=None):
         """
         General solution for fluid temperatures at a depth (z).
 
@@ -1484,6 +1650,11 @@ def _general_solution(self, z, m_flow_borehole, cp_f, nSegments):
             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
         -------
@@ -1494,7 +1665,8 @@ def _general_solution(self, z, m_flow_borehole, cp_f, nSegments):
 
         """
         # Check if model variables need to be updated
-        self._check_model_variables(m_flow_borehole, cp_f, nSegments)
+        self._check_model_variables(
+            m_flow_borehole, cp_f, nSegments, segment_ratios)
 
         # Load coefficients
         sumA = self._sumA
@@ -1508,14 +1680,16 @@ def _general_solution(self, z, m_flow_borehole, cp_f, nSegments):
 
         # Coefficient matrix for borehole wall temperatures
         a_b = np.zeros((2*self.nPipes, nSegments))
-        dz = z - np.minimum(z, np.arange(nSegments+1)*self.b.H/nSegments)
+        z_all = self.b._segment_edges(nSegments, segment_ratios=segment_ratios)
+        dz = z - np.minimum(z, z_all)
         exp_Lz = np.exp(np.multiply.outer(L, dz))
         dexp_Lz = exp_Lz[:,1:] - exp_Lz[:,:-1]
         a_b = np.real(((V @ Dm1) * (Vm1 @ sumA)) @ dexp_Lz)
 
         return a_f0, a_b
 
-    def _update_model_variables(self, m_flow_borehole, cp_f, nSegments):
+    def _update_model_variables(
+            self, m_flow_borehole, cp_f, nSegments, segment_ratios):
         """
         Evaluate eigenvalues and eigenvectors for the system of differential
         equations.
@@ -1528,11 +1702,15 @@ def _update_model_variables(self, m_flow_borehole, cp_f, nSegments):
             Fluid specific isobaric heat capacity (in J/kg.degC).
         nSegments : int
             Number of borehole segments.
+        segment_ratios : (nSegments,) array
+            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.
         """
 
         nPipes = self.nPipes
         # Format mass flow rate and heat capacity inputs
-        self._format_inputs(m_flow_borehole, cp_f, nSegments)
+        self._format_inputs(m_flow_borehole, cp_f, nSegments, segment_ratios)
         m_flow_pipe = self._m_flow_pipe
         cp_pipe = self._cp_pipe
 
@@ -1552,7 +1730,7 @@ def _update_model_variables(self, m_flow_borehole, cp_f, nSegments):
         self._D = np.diag(self._L)
         self._Dm1 = np.diag(1./self._L)
 
-    def _format_inputs(self, m_flow_borehole, cp_f, nSegments):
+    def _format_inputs(self, m_flow_borehole, cp_f, nSegments, segment_ratios):
         """
         Format mass flow rate and heat capacity inputs.
 
@@ -1564,6 +1742,10 @@ def _format_inputs(self, m_flow_borehole, cp_f, nSegments):
             Fluid specific isobaric heat capacity (in J/kg.degC).
         nSegments : int
             Number of borehole segments.
+        segment_ratios : (nSegments,) array
+            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.
         """
 
         nPipes = self.nPipes
@@ -1661,7 +1843,8 @@ def __init__(self, pos, r_in, r_out, borehole, k_s,
         # Initialize stored_coefficients
         self._initialize_stored_coefficients()
 
-    def _continuity_condition_base(self, m_flow_borehole, cp_f, nSegments):
+    def _continuity_condition_base(
+            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).
@@ -1681,6 +1864,11 @@ def _continuity_condition_base(self, m_flow_borehole, cp_f, nSegments):
             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
         -------
@@ -1693,16 +1881,18 @@ def _continuity_condition_base(self, m_flow_borehole, cp_f, nSegments):
 
         """
         # Check if model variables need to be updated
-        self._check_model_variables(m_flow_borehole, cp_f, nSegments)
+        self._check_model_variables(
+            m_flow_borehole, cp_f, nSegments, segment_ratios=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)
+            m_flow_borehole, cp_f, nSegments, segment_ratios=segment_ratios)
 
         return a_in, a_out, a_b
 
-    def _continuity_condition_head(self, m_flow_borehole, cp, nSegments):
+    def _continuity_condition_head(
+            self, m_flow_borehole, cp, nSegments, segment_ratios=None):
         """
         Build coefficient matrices to evaluate fluid temperatures at depth
         (z = 0). These coefficients take into account connections between
@@ -1724,6 +1914,11 @@ def _continuity_condition_head(self, m_flow_borehole, cp, nSegments):
             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
         -------
@@ -1736,7 +1931,8 @@ def _continuity_condition_head(self, m_flow_borehole, cp, nSegments):
 
         """
         # Check if model variables need to be updated
-        self._check_model_variables(m_flow_borehole, cp, nSegments)
+        self._check_model_variables(
+            m_flow_borehole, cp, 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)
@@ -1744,7 +1940,7 @@ def _continuity_condition_head(self, m_flow_borehole, cp, nSegments):
 
         return a_in, a_out, a_b
 
-    def _format_inputs(self, m_flow_borehole, cp_f, nSegments):
+    def _format_inputs(self, m_flow_borehole, cp_f, nSegments, segment_ratios):
         """
         Format mass flow rate and heat capacity inputs.
         """
diff --git a/pygfunction/utilities.py b/pygfunction/utilities.py
index 3c79f364..9084159d 100644
--- a/pygfunction/utilities.py
+++ b/pygfunction/utilities.py
@@ -1,6 +1,142 @@
 # -*- coding: utf-8 -*-
 import matplotlib.pyplot as plt
 import numpy as np
+import numpy.polynomial.polynomial as poly
+import warnings
+
+
+def segment_ratios(nSegments, end_length_ratio=0.02):
+    """
+    Discretize a borehole into segments of different lengths using a
+    geometrically expanding mesh from the provided end-length-ratio towards the
+    middle of the borehole. Eskilson (1987) [#Eskilson_1987]_ proposed that
+    segment lengths increase with a factor of sqrt(2) towards the middle of the
+    borehole. Here, the expansion factor is inferred from the provided number
+    of segments and end-length-ratio.
+
+    Parameters
+    ----------
+    nSegments : int
+        Number of line segments along the borehole.
+    end_length_ratio: float, optional
+        The ratio of the height of the borehole that accounts for the end
+        segment lengths.
+        Default is 0.02.
+
+    Returns
+    -------
+    segment_ratios : array
+        The segment ratios along the borehole, from top to bottom.
+
+    Examples
+    --------
+    >>> gt.utilities.segment_ratios(5)
+    array([0.02, 0.12, 0.72, 0.12, 0.02])
+
+    References
+    ----------
+    .. [#Eskilson_1987] Eskilson, P. (1987). Thermal analysis of heat
+        extraction boreholes. PhD Thesis. University of Lund, Department of
+        Mathematical Physics. Lund, Sweden.
+    """
+    def is_even(n):
+        "Returns True if n is even."
+        return not(n & 0x1)
+    assert nSegments >= 3 and isinstance(nSegments, int), \
+            "The number of segments `nSegments` should be greater or equal " \
+            "to 3 and of type int."
+    assert end_length_ratio > 0. and isinstance(end_length_ratio, (float, np.floating)), \
+            "The end-length-ratio `end_length_ratio` should be greater than " \
+            "0. and of type float."
+
+    # If nSegments == 3, then the middle segment is simply the remainder of the
+    # length
+    if nSegments == 3:
+        assert end_length_ratio < 0.5, \
+            "For nSegments == 3, the end-length-ratio `end_length_ratio` " \
+            "should be less than 0.5."
+        segment_ratios = np.array(
+            [end_length_ratio,
+             1 - 2 * end_length_ratio,
+             end_length_ratio])
+        return segment_ratios
+
+    # If end_length_ratio == 1 / nSegments, then the discretization is
+    # uniform
+    if np.abs(1. - nSegments * end_length_ratio) < 1e-6:
+        return np.full(nSegments, 1 / nSegments)
+
+    # Find the required constant expansion ratio to fill the borehole length
+    # from the provided end-length-ratio inwards with the provided nSegments
+    if is_even(nSegments):
+        # The ratio is a root of the polynomial expression :
+        # 0 = (1 - 2 * end_length_ratio)
+        #     - ratio * x
+        #     + 2 * end_length_ratio * x**nz
+        nz = int(nSegments / 2)
+        coefs = np.zeros(nz + 1)
+        coefs[0] = 1 - 2 * end_length_ratio
+        coefs[1] = -1
+        coefs[-1] = 2 * end_length_ratio
+        # Roots of the polynomial
+        roots = poly.Polynomial(coefs).roots()
+        # Find the correct root
+        for x in roots:
+            if np.isreal(x):
+                factor = np.real(x)
+                dz = [factor**i * end_length_ratio for i in range(nz)]
+                segment_ratios = np.concatenate(
+                    (dz,
+                     dz[::-1]))
+                if (np.abs(1. - np.sum(segment_ratios)) < 1e-6
+                    and np.all(segment_ratios > 0.)):
+                    break
+        else:
+            raise RuntimeError(
+                'utilities.segment_ratios failed to generate segment '
+                'discretization for the given input parameters : '
+                ' nSegments={}, end_length_ratio={}.'.format(
+                    nSegments, end_length_ratio))
+    else:
+        # The ratio is a root of the polynomial expression
+        # 0 = (1 - 2 * end_length_ratio) - ratio * x
+        #     + end_length_ratio * x**nz
+        #     + end_length_ratio * x**(nz + 1)
+        nz = int((nSegments - 1) / 2)
+        coefs = np.zeros(nz + 2)
+        coefs[0] = 1 - 2 * end_length_ratio
+        coefs[1] = -1
+        coefs[-2] = end_length_ratio
+        coefs[-1] = end_length_ratio
+        # Roots of the polynomial
+        roots = poly.Polynomial(coefs).roots()
+        # Find the correct root
+        for x in roots:
+            if np.isreal(x):
+                factor = np.real(x)
+                dz = [factor**i * end_length_ratio for i in range(nz)]
+                segment_ratios = np.concatenate(
+                    (dz,
+                     np.array([factor**nz]) * end_length_ratio,
+                     dz[::-1]))
+                if (np.abs(1. - np.sum(segment_ratios)) < 1e-6
+                    and np.all(segment_ratios > 0.)):
+                    break
+        else:
+            raise RuntimeError(
+                'utilities.segment_ratios failed to generate segment '
+                'discretization for the given input parameters : '
+                ' nSegments={}, end_length_ratio={}.'.format(
+                    nSegments, end_length_ratio))
+
+    if factor < 1.:
+        warnings.warn(
+            'A decreasing segment ratios discretization was found by '
+            'utilities.segment_ratios(). Better accuracy is expected for '
+            'increasing segment ratios. Consider decreasing the '
+            'end-length-ratio.')
+
+    return segment_ratios
 
 
 def time_ClaessonJaved(dt, tmax, cells_per_level=5):


Test patch:
diff --git a/tests/gfunction_test.py b/tests/gfunction_test.py
index 2df8ec35..0a36b58a 100644
--- a/tests/gfunction_test.py
+++ b/tests/gfunction_test.py
@@ -5,6 +5,147 @@
 
 import numpy as np
 
+import pygfunction.utilities
+
+
+class TestSolvers(unittest.TestCase):
+    """
+    Test cases for calculation of g-functions using various solver methods
+    available.
+    """
+
+    def setUp(self) -> None:
+        self.H = 150.  # Borehole length [m]
+        self.D = 4.  # Borehole buried depth [m]
+        self.r_b = 0.075  # Borehole radius [m]
+        self.B = 7.5  # Borehole spacing [m]
+        self.alpha = 1.0e-6  # Ground thermal diffusivity [m2/s]
+
+    def test_detailed(self):
+        from pygfunction.gfunction import gFunction
+        from pygfunction.boreholes import rectangle_field
+        from pygfunction.utilities import time_geometric, segment_ratios
+        # Calculation of the g-function at the same time values
+        N_1 = 3
+        N_2 = 2
+        boreField = rectangle_field(N_1, N_2, self.B, self.B,
+                                    self.H, self.D, self.r_b)
+
+        # Thermal properties
+        alpha = 1.0e-6  # Ground thermal diffusivity (m2/s)
+
+        # Geometrically expanding time vector.
+        dt = 100 * 3600.  # Time step
+        tmax = 3000. * 8760. * 3600.  # Maximum time
+        Nt = 10  # Number of time steps
+        time = time_geometric(dt, tmax, Nt)
+
+        # g-Function calculation options
+        options_uniform = {'nSegments': 24, 'disp': False}
+        ratios = segment_ratios(8)
+        options_unequal = {'nSegments': 8,
+                          'segment_ratios': ratios,
+                          'disp': False}
+        method = 'detailed'
+
+        gfunc_UBWT_uniform = gFunction(
+            boreField, alpha, time=time, options=options_uniform,
+            method=method).gFunc
+
+        gfunc_UBWT_unequal = gFunction(
+            boreField, alpha, time=time, options=options_unequal,
+            method=method).gFunc
+
+        self.assertTrue(np.allclose(gfunc_UBWT_uniform, gfunc_UBWT_unequal,
+                                    rtol=1e-2, atol=1e-6),
+                        msg='Incorrect g-function for the detailed solver test.'
+                            'A converged solution is compared to a discretized'
+                            'solution using the UBWT boundary condition.')
+
+    def test_similarities(self):
+        from pygfunction.gfunction import gFunction
+        from pygfunction.boreholes import rectangle_field
+        from pygfunction.utilities import time_geometric, segment_ratios
+        # Calculation of the g-function at the same time values
+        N_1 = 3
+        N_2 = 2
+        boreField = rectangle_field(N_1, N_2, self.B, self.B,
+                                    self.H, self.D, self.r_b)
+
+        # Thermal properties
+        alpha = 1.0e-6  # Ground thermal diffusivity (m2/s)
+
+        # Geometrically expanding time vector.
+        dt = 100 * 3600.  # Time step
+        tmax = 3000. * 8760. * 3600.  # Maximum time
+        Nt = 10  # Number of time steps
+        time = time_geometric(dt, tmax, Nt)
+
+        # g-Function calculation options
+        options_uniform = {'nSegments': 24, 'disp': False}
+        ratios = segment_ratios(8)
+        options_unequal = {'nSegments': 8,
+                          'segment_ratios': ratios,
+                          'disp': False}
+        method = 'similarities'
+
+        gfunc_UBWT_uniform = gFunction(
+            boreField, alpha, time=time, options=options_uniform,
+            method=method).gFunc
+
+        gfunc_UBWT_unequal = gFunction(
+            boreField, alpha, time=time, options=options_unequal,
+            method=method).gFunc
+
+        self.assertTrue(np.allclose(gfunc_UBWT_uniform, gfunc_UBWT_unequal,
+                                    rtol=1e-2, atol=1e-6),
+                        msg='Incorrect g-function for the similarities solver '
+                            'test.'
+                            'A converged solution is compared to a discretized'
+                            'solution using the UT boundary condition.')
+
+    def test_equivalent(self):
+        from pygfunction.gfunction import gFunction
+        from pygfunction.boreholes import rectangle_field
+        from pygfunction.utilities import time_geometric, segment_ratios
+        # Calculation of the g-function at the same time values
+        N_1 = 3
+        N_2 = 2
+        boreField = rectangle_field(N_1, N_2, self.B, self.B,
+                                    self.H, self.D, self.r_b)
+
+        # Thermal properties
+        alpha = 1.0e-6  # Ground thermal diffusivity (m2/s)
+
+        # Geometrically expanding time vector.
+        dt = 100 * 3600.  # Time step
+        tmax = 3000. * 8760. * 3600.  # Maximum time
+        Nt = 10  # Number of time steps
+        time = time_geometric(dt, tmax, Nt)
+
+        # g-Function calculation options
+        options_uniform = {'nSegments': 24, 'disp': False}
+        ratios = segment_ratios(8)
+        options_unequal = {'nSegments': 8,
+                          'segment_ratios': ratios,
+                          'disp': False}
+        method = 'equivalent'
+
+        gfunc_UBWT_uniform = gFunction(
+            boreField, alpha, time=time, options=options_uniform,
+            method=method).gFunc
+
+        gfunc_UBWT_unequal = gFunction(
+            boreField, alpha, time=time, options=options_unequal,
+            method=method).gFunc
+
+        self.assertTrue(np.allclose(gfunc_UBWT_uniform, gfunc_UBWT_unequal,
+                                    rtol=1e-2, atol=1e-6),
+                        msg='Incorrect g-function for the equivalent solver '
+                            'test.'
+                            'A converged solution is compared to a discretized'
+                            'solution using the UT boundary condition.')
+
 
 class TestUniformHeatExtractionRate(unittest.TestCase):
     """ Test cases for calculation of g-functions using uniform heat extraction
@@ -219,13 +360,13 @@ def setUp(self):
         self.epsilon = 1.0e-06  # Pipe roughness [m]
         self.m_flow_network = 0.25  # Fluid mass flor rate in network [kg/s]
 
-    def test_unequal_segments(self, rel_tol=1.0e-3):
+    def test_unequal_segments(self, rel_tol=1.0e-2):
         """ Tests the value of the g-function of six boreholes in series with
             unequal numbers of segments.
         """
         from pygfunction.gfunction import gFunction
         from pygfunction.boreholes import rectangle_field
-        from pygfunction.utilities import time_geometric
+        from pygfunction.utilities import time_geometric, segment_ratios
         from pygfunction.pipes import \
             conduction_thermal_resistance_circular_pipe, \
             convective_heat_transfer_coefficient_circular_pipe, \
@@ -269,25 +410,64 @@ def test_unequal_segments(self, rel_tol=1.0e-3):
             UTubes.append(UTube)
         network = Network(
             boreField, UTubes, bore_connectivity=bore_connectivity,
-            m_flow_network=m_flow_pipe, cp_f=fluid.cp, nSegments=nSegments)
+            m_flow_network=m_flow_pipe, cp_f=fluid.cp)
 
         # g-Function calculation option for uniform borehole segment lengths
         # in the field by defining nSegments as an integer >= 1
         options = {'nSegments': nSegments[0]}
         gfunc = gFunction(network, self.alpha, time=time, options=options)
-        g_ref = gfunc.gFunc
+        g_reference = gfunc.gFunc
 
         # g-Function calculation with nSegments passed as list, where the list
         # is of the same length as boreField and each borehole is defined to
         # have >= 1 segment
         # Note: nSegments[i] pertains to boreField[i]
-        options = {'nSegments':nSegments}
-        g = gFunction(network, self.alpha, time=time, options=options).gFunc
-
-        self.assertTrue(np.allclose(g, g_ref, rtol=rel_tol, atol=1e-6),
-                        msg='Incorrect values of the g-function of six '
-                            'boreholes for mixed inlet temperature and '
-                            'unequal numbers of segments.')
+        options = {'nSegments': nSegments}
+        g_unequal_nSegments = gFunction(
+            network, self.alpha, time=time, options=options).gFunc
+
+        self.assertTrue(
+            np.allclose(
+                g_unequal_nSegments, g_reference, rtol=rel_tol, atol=1e-6),
+            msg='Incorrect values of the g-function of six '
+                'boreholes for mixed inlet temperature and '
+                'unequal numbers of segments.')
+
+        # Compute g-function with predefined segment lengths
+        nSegments = 8
+        # Define the segment ratios for each borehole in each segment
+        # the segment lengths are defined top to bottom left to right
+        ratios = np.array(
+            [0.05, 0.10, 0.10, 0.25, 0.25, 0.10, 0.10, 0.05])
+        options = {'nSegments': nSegments,
+                   'segment_ratios': ratios,
+                   'disp': False}
+
+        g_predefined_ratios = gFunction(
+            network, self.alpha, time=time, options=options).gFunc
+
+        self.assertTrue(
+            np.allclose(
+                g_predefined_ratios, g_reference, rtol=rel_tol, atol=1e-6),
+            msg='Incorrect values of the g-function of six '
+                'boreholes for mixed inlet temperature and '
+                'unequal numbers of segments.')
+
+        # Compute g-function with discretized segment lengths
+        ratios = segment_ratios(8)
+        options = {'nSegments': nSegments,
+                   'segment_ratios': ratios,
+                   'disp': False}
+
+        g_generated_ratios = gFunction(
+            network, self.alpha, time=time, options=options).gFunc
+
+        self.assertTrue(
+            np.allclose(
+                g_generated_ratios, g_reference, rtol=rel_tol, atol=1e-6),
+            msg='Incorrect values of the g-function of six '
+                'boreholes for mixed inlet temperature and '
+                'unequal numbers of segments.')
 
 
 class TestEqualInletTemperature(unittest.TestCase):
@@ -348,7 +528,7 @@ def test_one_borehole_twenty_four_segments(self, rel_tol=1.0e-4):
                 UTube._Rd[1,1] = Rd_00
                 UTube._Rd[0,1] = Rd_01
                 UTube._Rd[1,0] = Rd_01
-                UTube._update_model_variables(self.m_flow, self.cp, nSegments)
+                UTube._update_model_variables(self.m_flow, self.cp, nSegments, None)
             g[i] = equal_inlet_temperature(boreField, UTubes, self.m_flow,
                                            self.cp, time, self.alpha,
                                            nSegments=nSegments)
@@ -400,7 +580,7 @@ def test_four_by_four_twenty_four_segments(self, rel_tol=1.0e-4):
             UTube._Rd[1,1] = Rd_00
             UTube._Rd[0,1] = Rd_01
             UTube._Rd[1,0] = Rd_01
-            UTube._update_model_variables(self.m_flow, self.cp, nSegments)
+            UTube._update_model_variables(self.m_flow, self.cp, nSegments, None)
         g = equal_inlet_temperature(boreField, UTubes, self.m_flow,
                                        self.cp, time, self.alpha,
                                        nSegments=nSegments)