Repository: MassimoCimmino/pygfunction (language: python)

Issue:
g-Functions for variable fluid mass flow rate and reversible flow direction
I have more a question, which may result in an enhancment: Is it already possible to simulate a borefield with given heat extraction rate and different network flow values per time step?

my script is based on the example:
https://pygfunction.readthedocs.io/en/stable/examples/fluid_temperature_multiple_boreholes.html#simulation-of-fluid-temperatures-in-a-field-of-multiple-boreholes

until now I used constant network flow, and the algorithm works really well eg:
_m_flow_network_w = 0.35* N1_w * N2_w_

(milesteps of the script):
_LoadAgg_w = gt.load_aggregation.ClaessonJaved(dt, tmax)
boreField_w = gt.boreholes.rectangle_field(N1_w, N2_w, B, B, H_w, D, r_b)
network_w = gt.networks.Network(boreField_w, UTubes_w, m_flow_network=m_flow_network_w, cp_f=cp_fluid)
gFunc_w = gt.gfunction.gFunction(network_w, alpha, time=time_req, boundary_condition='MIFT', options=options)
LoadAgg_w.initialize(gFunc_w.gFunc/(2*pi*k_s))
for i, (t) in enumerate(time):
    LoadAgg_w.next_time_step(t)
    Tf_in_w[i] = network_w.get_network_inlet_temperature(P_tot_w[i], Tb_w[i], m_flow_network_w, cp_fluid, nSegments=1)_

Now, I want to implement 3 different mass flow rates  (e.g. 0.1, 0.25 and 0.5 kg/s per borehole), depending on actual power value (P_tot_w). The output should be inlet and outlet temperature of the field.
as the gfunctions itself and hence the initialization of the Load Aggregation is dependent on the network flow,  I think it do not work properly  only changing the flow value in "_get_network_inlet_temperature_"... 
is there already a possibility to use several network flow steps?
thank you


Discussion / hints:
Hello @fuma01.

We found in an earlier study that *g*-functions are not strongly dependent on mass flow rate (as long as it is above a threshold value) [1]. For parallel-connected boreholes, the variation is often negligible [2]. This means that the same *g*-function (at nominal mass flow rate) can be used for a simulation with variable mass flow rate.

One way to verify is to evaluate and compare *g*-functions at your mininum an maximum mass flow rates up to a time equal to the maximum simulation time. If there is no significant difference, then all you need is to use the variable mass flow rate when calling `network.get_<>` class methods.

It is also possible to consider different *g*-functions with different mass flow rates using a non-stationary combination [3]. This is not supported by pygfunction and would need to be adapted to time-marching simulations before it can be implemented.

-------------------------------------
[1] Cimmino, M. (2019). Semi-Analytical Method for *g*-Function Calculation of bore fields with series-and parallel-connected boreholes. *Science and Technology for the Built Environment, 25*(8), 1007-1022.
[2] Cimmino, M. (2015). The effects of borehole thermal resistances and fluid flow rate on the g-functions of geothermal bore fields. *International Journal of Heat and Mass Transfer, 91*, 1119-1127.
[3] Beaudry, G., Pasquier, P., & Marcotte, D. (2021). A fast convolution-based method to simulate time-varying flow rates in closed-loop and standing column well ground heat exchangers. *Renewable Energy, 174*, 55-72.
Dear Massimo,

Thank you for the quick response!
Indeed the g-function is not very sensitiv to flow values, at least if the field stays in turbulent mode (< 0.6 % in my simulation) but a bit more in low flow mode (up to 10 % difference)
I testet it with three flow rates (0.1, 0.3 and 0.5 kg/s), which represents my target range and plottet the g-function and the difference to 0.5 kg/s, related to Eskilson time:
![GFUNConFLOW](https://user-images.githubusercontent.com/43027994/177823645-46e2a253-767f-4652-9b37-e2aa4fd2c1d0.jpg)
I use hourly time steps up to 20 years. rectangular borefield 7x9 with 150 m length and 5 m spacing

now I am thinking about, if the error will be significant when using variable mass flow rate just at"network.get_" method, as suggested? 
also, I am wondering, that the computation time needs 5 times longer with switching flow rates @ network.get.

regards from Vienna,
Martin

 
### Accuracy
The accuracy of your simulation will depend on your simulation objectives and how often the flow rate is 0.1 kg/s.

A difference in *g*-functions at low `ln(t/ts)` translates to an error in borehole wall temperature variations in the short-term (i.e. how quickly the temperature varies following a heat pulse). The variation of fluid temperatures is a combination of two factors: (1) the temperature variation at the borehole wall, and (2) the effective borehole thermal resistance. The second also varies with mass flow and is captured by the `Network` class.

If you are simulating for heat pump energy use, then an error in temperature affects the performance of the heat pump. If you are simulating to size the field in peak conditions then I figure the peak coincides with the maximum mass flow rate. The magnitude of these impacts strongly depend on your simulation case.

Here, I would suggest running and comparing simulations using *g*-functions evaluated with the 0.1 and 0.5 kg/s mass flow rates. This would give an estimation of the expected magnitude of the errors caused by using only one *g*-function.

Factor in that *g*-functions do not capture short-term effects (i.e. the heat capacity of the grout and the travel time of the fluid). This is a limitation of steady-state resistance models such as the one implemented in the `Network` class.

### Calculation time
When the mass flow rate changes, the `Network` class re-evaluates all borehole coefficients. the class only remembers coefficients for the latest mass flow rate.

At the moment, this is done to all boreholes individually even if all borehole are geometrically identical.
in my current investigation, I have 5 heat pumps parallel to one (or two) BHE fields. Input is the heating/cooling power demand of the building and solar module. Then the field size is adapted, that the fluid temperatures are in specified range. The output is the electric power of the heat pumps (calculated with characteristic curve) in proportion to solar and BHE power and relation to field size.
In dependence of the requested power load the needed number of heat pumps are activated. In current simulation the mass flow rate is set constant, hence the delta T is varying. The idea was to set a variably flow rate and get more or less constant delta T, which will also affect input temperature of the heat pump and hence the COP (or the needed electric power).

now, I understand, that the mass flow rate is implementet in the g-function, which is calculated in advance. By switching laminar/turbulent flow mode, this will probably not be neclectable. As a workaround, I will make two simulations at minimum and maximum constant flow rate, to see the efffects. 
The increase of calculation time by changing flow rate later in network class, is too high for me...

I think, it would be exciting if an implementation of "different g-functions with different mass flow rates using a non-stationary combination [3]" in pygfunction would be implemented in future releases! 

Fix patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index f74c9122..29968344 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -4,6 +4,7 @@
 
 ### New features
 
+* [Issue 215](https://github.com/MassimoCimmino/pygfunction/issues/215) - Implemented variable fluid mass flow rate g-functions. Bore fields with series-connected boreholes and reversible flow direction can now be simulated.
 * [Issue 282](https://github.com/MassimoCimmino/pygfunction/issues/282) - Enabled the use of negative mass flow rates in `Pipe` and `Network` classes to model reversed flow direction.
 
 ## Version 2.2.3 (2024-07-01)
diff --git a/doc/source/example_index.rst b/doc/source/example_index.rst
index 17e869c3..046e89d9 100644
--- a/doc/source/example_index.rst
+++ b/doc/source/example_index.rst
@@ -21,6 +21,7 @@ Examples
    examples/fluid_temperature
    examples/equal_inlet_temperature
    examples/fluid_temperature_multiple_boreholes
+   examples/fluid_temperature_reversible_flow_direction
    examples/multiple_independent_Utubes
    examples/multipole_temperature
    examples/bore_field_thermal_resistance
diff --git a/doc/source/examples/fluid_temperature_reversible_flow_direction.rst b/doc/source/examples/fluid_temperature_reversible_flow_direction.rst
new file mode 100644
index 00000000..74a9d85b
--- /dev/null
+++ b/doc/source/examples/fluid_temperature_reversible_flow_direction.rst
@@ -0,0 +1,30 @@
+.. examples:
+
+*******************************************************************************************************
+Simulation of fluid temperatures in a field of series-connected boreholes and reversible flow direction
+*******************************************************************************************************
+
+This example demonstrates the use of the
+:doc:`networks <../modules/networks>` module to predict the fluid temperature variations
+in a bore field with series-connected boreholes and reversible flow direction.
+
+The variable fluid mass flow rates g-functions of a bore field are first calculated
+using the mixed inlet fluid temperature boundary condition [1]_. Then, the effective borehole
+wall temperature variations are calculated using the load aggregation scheme of Claesson and
+Javed [2]_. g-Functions used in the temporal superposition are interpolated with regards to
+the fluid mass flow rate at the moment of heat extraction.
+
+The script is located in:
+`pygfunction/examples/fluid_temperature_reversible_flow_direction.py`
+
+.. literalinclude:: ../../../examples/fluid_temperature_reversible_flow_direction.py
+   :language: python
+   :linenos:
+
+.. rubric:: References
+.. [1] Cimmino, M. (2024). g-Functions for fields of
+   series- and parallel-connected boreholes with variable fluid mass flow
+   rate and reversible flow direction. Renewable Energy 228: 120661.
+.. [2] Claesson, J., & Javed, S. (2011). A load-aggregation method to calculate
+   extraction temperatures of borehole heat exchangers. ASHRAE Transactions,
+   118 (1): 530–539.
diff --git a/examples/discretize_boreholes.py b/examples/discretize_boreholes.py
index cd1768be..95b98eda 100644
--- a/examples/discretize_boreholes.py
+++ b/examples/discretize_boreholes.py
@@ -114,8 +114,7 @@ def main():
     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)
+    network = gt.networks.Network(boreField, UTubes)
 
     # -------------------------------------------------------------------------
     # Evaluate the g-functions for the borefield
@@ -124,8 +123,8 @@ def main():
     # 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)
+        network, alpha, time=time, m_flow_network=m_flow_network, cp_f=cp_f,
+        boundary_condition='MIFT', options=options_uniform)
 
     # Calculate the g-function for uniform borehole wall temperature
     gfunc_UBWT_uniform = gt.gfunction.gFunction(
@@ -135,8 +134,8 @@ def main():
     # 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)
+        network, alpha, time=time, m_flow_network=m_flow_network, cp_f=cp_f,
+        boundary_condition='MIFT', options=options_unequal)
 
     # Calculate the g-function for uniform borehole wall temperature
     gfunc_UBWT_unequal = gt.gfunction.gFunction(
diff --git a/examples/equal_inlet_temperature.py b/examples/equal_inlet_temperature.py
index 4f21bb57..b39cd65a 100644
--- a/examples/equal_inlet_temperature.py
+++ b/examples/equal_inlet_temperature.py
@@ -93,13 +93,11 @@ def main():
     # 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)
+        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 = gt.networks.Network(
-        boreField, UTubes, m_flow_network=m_flow_network, cp_f=cp_f,
-        nSegments=nSegments)
+    m_flow_network = m_flow_borehole * nBoreholes
+    network = gt.networks.Network(boreField, UTubes)
 
     # -------------------------------------------------------------------------
     # Evaluate the g-functions for the borefield
@@ -107,23 +105,26 @@ def main():
 
     # Calculate the g-function for uniform heat extraction rate
     gfunc_uniform_Q = gt.gfunction.gFunction(
-        boreField, alpha, time=time, boundary_condition='UHTR', options=options)
+        boreField, alpha, time=time, boundary_condition='UHTR',
+        options=options)
 
     # Calculate the g-function for uniform borehole wall temperature
     gfunc_uniform_T = gt.gfunction.gFunction(
-        boreField, alpha, time=time, boundary_condition='UBWT', options=options)
+        boreField, alpha, time=time, boundary_condition='UBWT',
+        options=options)
 
     # Calculate the g-function for equal inlet fluid temperature
     gfunc_equal_Tf_in = gt.gfunction.gFunction(
-        network, alpha, time=time, boundary_condition='MIFT', options=options)
+        network, alpha, time=time, m_flow_network=m_flow_network, cp_f=cp_f,
+        boundary_condition='MIFT', options=options)
 
     # -------------------------------------------------------------------------
     # Plot g-functions
     # -------------------------------------------------------------------------
 
     ax = gfunc_uniform_Q.visualize_g_function().axes[0]
-    ax.plot(np.log(time/ts), gfunc_uniform_T.gFunc, 'k--')
-    ax.plot(np.log(time/ts), gfunc_equal_Tf_in.gFunc, 'r-.')
+    ax.plot(np.log(time / ts), gfunc_uniform_T.gFunc, 'k--')
+    ax.plot(np.log(time / ts), gfunc_equal_Tf_in.gFunc, 'r-.')
     ax.legend(['Uniform heat extraction rate',
                'Uniform borehole wall temperature',
                'Equal inlet temperature'])
diff --git a/examples/fluid_temperature_multiple_boreholes.py b/examples/fluid_temperature_multiple_boreholes.py
index 93b8c748..4b8685c5 100644
--- a/examples/fluid_temperature_multiple_boreholes.py
+++ b/examples/fluid_temperature_multiple_boreholes.py
@@ -55,7 +55,8 @@ def main():
 
     # Fluid properties
     m_flow_borehole = 0.25      # Total fluid mass flow rate per borehole (kg/s)
-    m_flow_network = m_flow_borehole*N_1*N_2    # Total fluid mass flow rate (kg/s)
+    # Total fluid mass flow rate (kg/s)
+    m_flow_network = m_flow_borehole * N_1 * N_2
     # 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)
@@ -102,8 +103,7 @@ def main():
             nPipes=2, config='parallel')
         UTubes.append(UTube)
     # Build a network object from the list of UTubes
-    network = gt.networks.Network(
-        boreField, UTubes, m_flow_network=m_flow_network, cp_f=cp_f)
+    network = gt.networks.Network(boreField, UTubes)
 
     # -------------------------------------------------------------------------
     # Calculate g-function
@@ -113,17 +113,17 @@ def main():
     time_req = LoadAgg.get_times_for_simulation()
     # Calculate g-function
     gFunc = gt.gfunction.gFunction(
-        network, alpha, time=time_req, boundary_condition='MIFT',
-        options=options)
+        network, alpha, time=time_req, m_flow_network=m_flow_network,
+        cp_f=cp_f, boundary_condition='MIFT', options=options)
     # Initialize load aggregation scheme
-    LoadAgg.initialize(gFunc.gFunc/(2*pi*k_s))
+    LoadAgg.initialize(gFunc.gFunc / (2 * pi * k_s))
 
     # -------------------------------------------------------------------------
     # Simulation
     # -------------------------------------------------------------------------
 
     # Evaluate heat extraction rate
-    Q_tot = nBoreholes*synthetic_load(time/3600.)
+    Q_tot = nBoreholes*synthetic_load(time / 3600.)
 
     T_b = np.zeros(Nt)
     T_f_in = np.zeros(Nt)
@@ -192,7 +192,6 @@ def main():
     T_f = UTubes[0].get_temperature(
         z, T_f_in[it], T_b[it], m_flow_borehole, cp_f)
 
-
     # Configure figure and axes
     fig = gt.utilities._initialize_figure()
 
diff --git a/examples/fluid_temperature_reversible_flow_direction.py b/examples/fluid_temperature_reversible_flow_direction.py
new file mode 100644
index 00000000..e722f413
--- /dev/null
+++ b/examples/fluid_temperature_reversible_flow_direction.py
@@ -0,0 +1,231 @@
+# -*- coding: utf-8 -*-
+""" Example of simulation of a geothermal system with multiple series-
+    connected boreholes and reversible flow direction.
+
+    The g-function of a bore field is calculated for boundary condition of
+    mixed inlet fluid temperature into the boreholes and a combination of two
+    fluid mass flow rates in two opposing flow direction. Then, the borehole
+    wall temperature variations resulting from a time-varying load profile
+    are simulated using the aggregation method of Claesson and Javed (2012).
+    Predicted outlet fluid temperatures of double U-tube borehole are
+    evaluated.
+
+"""
+import matplotlib.pyplot as plt
+import numpy as np
+from scipy.constants import pi
+
+import pygfunction as gt
+
+
+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)
+
+    # Bore field geometry (rectangular array)
+    N_1 = 1             # Number of boreholes in the x-direction (columns)
+    N_2 = 6             # Number of boreholes in the y-direction (rows)
+    B = 7.5             # Borehole spacing, in both directions (m)
+
+    # Pipe dimensions
+    r_out = 33.6e-3/2   # Pipe outer radius (m)
+    r_in = 27.4e-3/2    # Pipe inner radius (m)
+    D_s = 0.11/2        # 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 = 2.5/2.2e6   # Ground thermal diffusivity (m2/s)
+    k_s = 2.5           # Ground thermal conductivity (W/m.K)
+    T_g = 10.           # Undisturbed ground temperatue (degC)
+
+    # Grout properties
+    k_g = 1.5           # Grout thermal conductivity (W/m.K)
+
+    # Pipe properties
+    k_p = 0.42          # Pipe thermal conductivity (W/m.K)
+
+    # Fluid properties
+    m_flow_borehole_min = 0.5   # Minimum fluid mass flow rate per borehole (kg/s)
+    m_flow_borehole_max = 1.2   # Maximum fluid mass flow rate per borehole (kg/s)
+    m_flow_borehole = np.array([
+        -m_flow_borehole_max,
+        -m_flow_borehole_min,
+        m_flow_borehole_min,
+        m_flow_borehole_max
+        ])
+    nModes = len(m_flow_borehole)
+    # Total fluid mass flow rate (kg/s)
+    m_flow_network = m_flow_borehole
+    # 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
+    options = {'approximate_FLS': True,
+               'nSegments': 8,
+               'disp': True}
+
+    # Simulation parameters
+    dt = 600.                   # Time step (s)
+    tmax = 48 * 3600.           # Maximum time (s)
+    Nt = int(np.ceil(tmax/dt))  # Number of time steps
+    time = dt * np.arange(1, Nt+1)
+    time_h = time / 3600.
+
+    # Load aggregation scheme
+    LoadAgg = gt.load_aggregation.ClaessonJaved(dt, tmax, nSources=nModes)
+
+    # -------------------------------------------------------------------------
+    # Initialize bore field and pipe models
+    # -------------------------------------------------------------------------
+
+    # The field is a retangular array
+    boreField = gt.boreholes.rectangle_field(N_1, N_2, B, B, H, D, r_b)
+    nBoreholes = len(boreField)
+    H_tot = np.sum([b.H for b in boreField])
+
+    # Boreholes are connected in series
+    bore_connectivity = [i-1 for i in range(nBoreholes)]
+
+    # 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 (Double U-tube in parallel)
+    m_flow_pipe = np.max(np.abs(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)
+
+    # Double U-tube (parallel), same for all boreholes in the bore field
+    UTubes = []
+    for borehole in boreField:
+        UTube = gt.pipes.SingleUTube(
+            pos_pipes, r_in, r_out, borehole, k_s, k_g, R_f + R_p)
+        UTubes.append(UTube)
+    # Build a network object from the list of UTubes
+    network = gt.networks.Network(
+        boreField, UTubes, bore_connectivity=bore_connectivity)
+
+    # -------------------------------------------------------------------------
+    # Calculate g-function
+    # -------------------------------------------------------------------------
+
+    # Get time values needed for g-function evaluation
+    time_req = LoadAgg.get_times_for_simulation()
+    # Calculate g-function
+    gFunc = gt.gfunction.gFunction(
+        network, alpha, time=time_req, boundary_condition='MIFT',
+        m_flow_network=m_flow_network, cp_f=cp_f, options=options)
+    # Initialize load aggregation scheme
+    LoadAgg.initialize(gFunc.gFunc / (2 * pi * k_s))
+
+    # -------------------------------------------------------------------------
+    # Simulation
+    # -------------------------------------------------------------------------
+
+    # Sinusoidal variation between -80 and 80 W/m with minimum of 10 W/m
+    Q_min = -80. * H * nBoreholes
+    Q_max = 80. * H * nBoreholes
+    Q_small = 10. * H * nBoreholes
+    Q_avg = 0.5 * (Q_min + Q_max)
+    Q_amp = 0.5 * (Q_max - Q_min)
+    Q_tot = -Q_avg - Q_amp * np.sin(2 * np.pi * time_h / 24)
+    # Mass flow rate is proportional to the absolute load
+    m_flow_min = np.min(np.abs(m_flow_network))
+    m_flow_max = np.max(np.abs(m_flow_network))
+    m_flow = np.sign(Q_tot) * np.maximum(
+        m_flow_min,
+        np.minimum(m_flow_max, m_flow_max * np.abs(Q_tot) / Q_max))
+    # The minimum heat load is 10 W/m and the mass flow rate is 0 when Q_tot=0
+    Q_tot[np.abs(Q_tot) < Q_small] = 0.
+    m_flow[np.abs(Q_tot) < Q_small] = 0.
+
+    T_b = np.zeros(Nt)
+    T_f_in = np.zeros(Nt)
+    T_f_out = np.zeros(Nt)
+    for i, (t, Q_i) in enumerate(zip(time, Q_tot)):
+        # Increment time step by (1)
+        LoadAgg.next_time_step(t)
+
+        # Solve for temperatures and heat extraction rates (variable mass flow rate)
+        if np.abs(Q_tot[i]) >= Q_small:
+            xi = np.zeros(nModes)
+            iMode = np.argmax(m_flow[i] <= m_flow_network)
+            if iMode > 0:
+                xi[iMode-1] = (m_flow_network[iMode] - m_flow[i]) / (m_flow_network[iMode] - m_flow_network[iMode-1])
+                xi[iMode] = (m_flow[i] - m_flow_network[iMode-1]) / (m_flow_network[iMode] - m_flow_network[iMode-1])
+            else:
+                xi[0] = 1.
+            LoadAgg.set_current_load(Q_tot[i] / H_tot * xi)
+            deltaT_b = LoadAgg.temporal_superposition()
+            T_b[i] = T_g - deltaT_b @ xi
+            T_f_in[i] = network.get_network_inlet_temperature(
+                Q_tot[i], T_b[i], m_flow[i], cp_f, nSegments=1, segment_ratios=None)
+            T_f_out[i] = network.get_network_outlet_temperature(
+                T_f_in[i], T_b[i], m_flow[i], cp_f, nSegments=1, segment_ratios=None)
+        else:
+            T_b[i] = np.nan
+            T_f_in[i] = np.nan
+            T_f_out[i] = np.nan
+
+    # -------------------------------------------------------------------------
+    # Plot hourly heat extraction rates and temperatures
+    # -------------------------------------------------------------------------
+
+    # Configure figure and axes
+    fig = gt.utilities._initialize_figure()
+
+    ax1 = fig.add_subplot(221)
+    # Axis labels
+    ax1.set_xlabel(r'Time [hours]')
+    ax1.set_ylabel(r'Total heat extraction rate [W/m]')
+    gt.utilities._format_axes(ax1)
+
+    # Plot heat extraction rates
+    hours = np.arange(1, Nt+1) * dt / 3600.
+    ax1.plot(hours, Q_tot / H_tot)
+
+    ax2 = fig.add_subplot(222)
+    # Axis labels
+    ax2.set_xlabel(r'Time [hours]')
+    ax2.set_ylabel(r'Fluid mass flow rate [kg/s]')
+    gt.utilities._format_axes(ax2)
+
+    # Plot temperatures
+    ax2.plot(hours, m_flow)
+
+    ax3 = fig.add_subplot(223)
+    # Axis labels
+    ax3.set_xlabel(r'Time [hours]')
+    ax3.set_ylabel(r'Temperature [degC]')
+    gt.utilities._format_axes(ax3)
+
+    # Plot temperatures
+    ax3.plot(hours, T_b, label='Borehole wall')
+    ax3.plot(hours, T_f_out, '-.',
+             label='Outlet')
+    ax3.legend()
+
+    # Adjust to plot window
+    plt.tight_layout()
+
+    return
+
+
+# Main function
+if __name__ == '__main__':
+    main()
diff --git a/examples/mixed_inlet_conditions.py b/examples/mixed_inlet_conditions.py
index d21ad34b..b34fec2f 100644
--- a/examples/mixed_inlet_conditions.py
+++ b/examples/mixed_inlet_conditions.py
@@ -50,7 +50,8 @@ def main():
     k_p = 0.4           # Pipe thermal conductivity (W/m.K)
 
     # Fluid properties
-    m_flow_network = 0.25   # Total fluid mass flow rate in network (kg/s)
+    # Total fluid mass flow rate in network (kg/s)
+    m_flow_network = np.array([-0.25, 0.5])   
     # All boreholes are in series
     m_flow_borehole = m_flow_network
     # The fluid is propylene-glycol (20 %) at 20 degC
@@ -98,7 +99,7 @@ def main():
     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
+    m_flow_pipe = np.max(np.abs(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)
@@ -110,8 +111,7 @@ def main():
             pos_pipes, r_in, r_out, borehole, k_s, k_g, R_f + R_p)
         UTubes.append(SingleUTube)
     network = gt.networks.Network(
-        boreField, UTubes, bore_connectivity=bore_connectivity,
-        m_flow_network=m_flow_network, cp_f=cp_f, nSegments=nSegments)
+        boreField, UTubes, bore_connectivity=bore_connectivity)
 
     # -------------------------------------------------------------------------
     # Evaluate the g-functions for the borefield
@@ -124,16 +124,20 @@ def main():
 
     # Calculate the g-function for mixed inlet fluid conditions
     gfunc_equal_Tf_mixed = gt.gfunction.gFunction(
-        network, alpha, time=time, boundary_condition='MIFT', options=options,
-        method=method)
+        network, alpha, time=time, m_flow_network=m_flow_network, cp_f=cp_f,
+        boundary_condition='MIFT', options=options, method=method)
 
     # -------------------------------------------------------------------------
     # Plot g-functions
     # -------------------------------------------------------------------------
 
     ax = gfunc_Tb.visualize_g_function().axes[0]
-    ax.plot(np.log(time/ts), gfunc_equal_Tf_mixed.gFunc, 'r-.')
-    ax.legend(['Uniform temperature', 'Mixed inlet temperature'])
+    ax.plot(np.log(time/ts), gfunc_equal_Tf_mixed.gFunc[0, 0, :], 'C1')
+    ax.plot(np.log(time/ts), gfunc_equal_Tf_mixed.gFunc[1, 1, :], 'C2')
+    ax.legend([
+        'Uniform temperature',
+        f'Mixed inlet temperature (m_flow={m_flow_network[0]} kg/s)',
+        f'Mixed inlet temperature (m_flow={m_flow_network[1]} kg/s)'])
     plt.tight_layout()
 
     # For the mixed inlet fluid temperature condition, draw the temperatures
diff --git a/pygfunction/gfunction.py b/pygfunction/gfunction.py
index abb73d65..e3877954 100644
--- a/pygfunction/gfunction.py
+++ b/pygfunction/gfunction.py
@@ -80,6 +80,27 @@ class gFunction(object):
 
         If not given, chosen to be 'UBWT' if a list of boreholes is provided
         or 'MIFT' if a Network object is provided.
+    m_flow_borehole : (nInlets,) array or (nMassFlow, nInlets,) array, optional
+        Fluid mass flow rate into each circuit of the network. If a
+        (nMassFlow, nInlets,) array is supplied, the
+        (nMassFlow, nMassFlow,) variable mass flow rate g-functions
+        will be evaluated using the method of Cimmino (2024)
+        [#gFunction-Cimmin2024]_. Only required for the 'MIFT' boundary
+        condition. Only one of 'm_flow_borehole' and 'm_flow_network' can be
+        provided.
+        Default is None.
+    m_flow_network : float or (nMassFlow,) array, optional
+        Fluid mass flow rate into the network of boreholes. If an array
+        is supplied, the (nMassFlow, nMassFlow,) variable mass flow
+        rate g-functions will be evaluated using the method of Cimmino
+        (2024) [#gFunction-Cimmin2024]_. Only required for the 'MIFT' boundary
+        condition. Only one of 'm_flow_borehole' and 'm_flow_network' can be
+        provided.
+        Default is None.
+    cp_f : float, optional
+        Fluid specific isobaric heat capacity (in J/kg.degC). Only required
+        for the 'MIFT' boundary condition.
+        Default is None.
     options : dict, optional
         A dictionary of solver options. All methods accept the following
         generic options:
@@ -184,6 +205,8 @@ class gFunction(object):
     - The g-function is linearized for times `t < r_b**2 / (25 * self.alpha)`.
       The g-function value is then interpolated between 0 and its value at the
       threshold.
+    - If the 'MIFT' boundary condition is used, only one of the
+      'm_flow_borehole' or 'm_flow_network' can be supplied.
 
     References
     ----------
@@ -210,14 +233,22 @@ class gFunction(object):
        finite line source solution to model thermal interactions between
        geothermal boreholes. International Communications in Heat and Mass
        Transfer, 127, 105496.
+    .. [#gFunction-Cimmin2024] Cimmino, M. (2024). g-Functions for fields of
+       series- and parallel-connected boreholes with variable fluid mass flow
+       rate and reversible flow direction. Renewable Energy, 228, 120661.
 
     """
     def __init__(self, boreholes_or_network, alpha, time=None,
-                 method='equivalent', boundary_condition=None, options={}):
+                 method='equivalent', boundary_condition=None,
+                 m_flow_borehole=None, m_flow_network=None,
+                 cp_f=None, options={}):
         self.alpha = alpha
         self.time = time
         self.method = method
         self.boundary_condition = boundary_condition
+        self.m_flow_borehole = m_flow_borehole
+        self.m_flow_network = m_flow_network
+        self.cp_f = cp_f
         self.options = options
 
         # Format inputs and assign default values where needed
@@ -229,15 +260,18 @@ def __init__(self, boreholes_or_network, alpha, time=None,
         if self.method.lower()=='similarities':
             self.solver = _Similarities(
                 self.boreholes, self.network, self.time,
-                self.boundary_condition, **self.options)
+                self.boundary_condition, self.m_flow_borehole,
+                self.m_flow_network, self.cp_f, **self.options)
         elif self.method.lower()=='detailed':
             self.solver = _Detailed(
                 self.boreholes, self.network, self.time,
-                self.boundary_condition, **self.options)
+                self.boundary_condition, self.m_flow_borehole,
+                self.m_flow_network, self.cp_f, **self.options)
         elif self.method.lower()=='equivalent':
             self.solver = _Equivalent(
                 self.boreholes, self.network, self.time,
-                self.boundary_condition, **self.options)
+                self.boundary_condition, self.m_flow_borehole,
+                self.m_flow_network, self.cp_f, **self.options)
         else:
             raise ValueError(f"'{method}' is not a valid method.")
 
@@ -283,10 +317,17 @@ def evaluate_g_function(self, time):
             print(60*'-')
         return self.gFunc
 
-    def visualize_g_function(self):
+    def visualize_g_function(self, which=None):
         """
         Plot the g-function of the borefield.
 
+        Parameters
+        ----------
+        which : list of tuple, optional
+            Tuples (i, j) of the variable mass flow rate g-functions to plot.
+            If None, all g-functions are plotted.
+            Default is None.
+
         Returns
         -------
         fig : figure
@@ -305,13 +346,34 @@ def visualize_g_function(self):
         # Dimensionless time (log)
         lntts = np.log(self.time/ts)
         # Draw g-function
-        ax.plot(lntts, self.gFunc)
+        if self.solver.nMassFlow == 0:
+            ax.plot(lntts, self.gFunc)
+        elif which is None:
+            for j in range(self.solver.nMassFlow):
+                for i in range(self.solver.nMassFlow):
+                    ax.plot(
+                        lntts,
+                        self.gFunc[i,j,:],
+                        label=f'$g_{{{i}{j}}}$')
+            plt.legend()
+        else:
+            if which is None:
+                which = [
+                    (i, j) for j in range(self.solver.nMassFlow)
+                    for i in range(self.solver.nMassFlow)]
+            for (i, j) in which:
+                ax.plot(
+                    lntts,
+                    self.gFunc[i,j,:],
+                    label=f'$g_{{{i}{j}}}$')
+            plt.legend()
 
         # Adjust figure to window
         plt.tight_layout()
         return fig
 
-    def visualize_heat_extraction_rates(self, iBoreholes=None, showTilt=True):
+    def visualize_heat_extraction_rates(
+            self, iBoreholes=None, showTilt=True, which=None):
         """
         Plot the time-variation of the average heat extraction rates.
 
@@ -325,6 +387,11 @@ def visualize_heat_extraction_rates(self, iBoreholes=None, showTilt=True):
         showTilt : bool
             Set to True to show borehole inclination.
             Default is True
+        which : list of int, optional
+            Indices i of the diagonal variable mass flow rate g-functions for
+            which to plot heat extraction rates.
+            If None, all diagonal g-functions are plotted.
+            Default is None.
 
         Returns
         -------
@@ -337,63 +404,120 @@ def visualize_heat_extraction_rates(self, iBoreholes=None, showTilt=True):
             iBoreholes = range(len(self.solver.boreholes))
         # Import heat extraction rates
         Q_t = self._heat_extraction_rates(iBoreholes)
-
-        # Configure figure and axes
-        fig = _initialize_figure()
-        ax1 = fig.add_subplot(121)
-        ax1.set_xlabel(r'$x$ [m]')
-        ax1.set_ylabel(r'$y$ [m]')
-        ax1.axis('equal')
-        _format_axes(ax1)
-        ax2 = fig.add_subplot(122)
-        ax2.set_xlabel(r'ln$(t/t_s)$')
-        ax2.set_ylabel(r'$\bar{Q}_b$')
-        _format_axes(ax2)
-
         # Borefield characteristic time
         ts = np.mean([b.H for b in self.solver.boreholes])**2/(9.*self.alpha)
         # Dimensionless time (log)
         lntts = np.log(self.time/ts)
-        # Plot curves for requested boreholes
-        for i, borehole in enumerate(self.solver.boreholes):
-            if i in iBoreholes:
-                # Draw heat extraction rate
-                line = ax2.plot(lntts, Q_t[iBoreholes.index(i)])
-                color = line[-1]._color
-                # Draw colored marker for borehole position
-                if showTilt:
-                    ax1.plot(
-                        [borehole.x, borehole.x + borehole.H*np.sin(borehole.tilt)*np.cos(borehole.orientation)],
-                        [borehole.y, borehole.y + borehole.H*np.sin(borehole.tilt)*np.sin(borehole.orientation)],
-                         linestyle='--',
-                         marker='None',
-                         color=color)
-                ax1.plot(borehole.x,
-                         borehole.y,
-                         linestyle='None',
-                         marker='o',
-                         color=color)
-            else:
-                # Draw black marker for borehole position
-                if showTilt:
-                    ax1.plot(
-                        [borehole.x, borehole.x + borehole.H*np.sin(borehole.tilt)*np.cos(borehole.orientation)],
-                        [borehole.y, borehole.y + borehole.H*np.sin(borehole.tilt)*np.sin(borehole.orientation)],
-                         linestyle='--',
-                         marker='None',
-                         color='k')
-                ax1.plot(borehole.x,
-                         borehole.y,
-                         linestyle='None',
-                         marker='o',
-                         color='k')
 
-        # Adjust figure to window
-        plt.tight_layout()
+        if self.solver.nMassFlow == 0:
+            # Configure figure and axes
+            fig = _initialize_figure()
+            ax1 = fig.add_subplot(121)
+            ax1.set_xlabel(r'$x$ [m]')
+            ax1.set_ylabel(r'$y$ [m]')
+            ax1.axis('equal')
+            _format_axes(ax1)
+            ax2 = fig.add_subplot(122)
+            ax2.set_xlabel(r'ln$(t/t_s)$')
+            ax2.set_ylabel(r'$\bar{Q}_b$')
+            _format_axes(ax2)
+
+            # Plot curves for requested boreholes
+            for i, borehole in enumerate(self.solver.boreholes):
+                if i in iBoreholes:
+                    # Draw heat extraction rate
+                    line = ax2.plot(lntts, Q_t[iBoreholes.index(i)])
+                    color = line[-1]._color
+                    # Draw colored marker for borehole position
+                    if showTilt:
+                        ax1.plot(
+                            [borehole.x, borehole.x + borehole.H*np.sin(borehole.tilt)*np.cos(borehole.orientation)],
+                            [borehole.y, borehole.y + borehole.H*np.sin(borehole.tilt)*np.sin(borehole.orientation)],
+                             linestyle='--',
+                             marker='None',
+                             color=color)
+                    ax1.plot(borehole.x,
+                             borehole.y,
+                             linestyle='None',
+                             marker='o',
+                             color=color)
+                else:
+                    # Draw black marker for borehole position
+                    if showTilt:
+                        ax1.plot(
+                            [borehole.x, borehole.x + borehole.H*np.sin(borehole.tilt)*np.cos(borehole.orientation)],
+                            [borehole.y, borehole.y + borehole.H*np.sin(borehole.tilt)*np.sin(borehole.orientation)],
+                             linestyle='--',
+                             marker='None',
+                             color='k')
+                    ax1.plot(borehole.x,
+                             borehole.y,
+                             linestyle='None',
+                             marker='o',
+                             color='k')
+
+            # Adjust figure to window
+            plt.tight_layout()
+        else:
+            m_flow = self.solver.m_flow
+            if which is None:
+                which = [n for n in range(self.solver.nMassFlow)]
+            for n in which:
+                # Configure figure and axes
+                fig = _initialize_figure()
+                fig.suptitle(
+                    f'Heat extraction rates for m_flow={m_flow[n]} kg/s')
+                ax1 = fig.add_subplot(121)
+                ax1.set_xlabel(r'$x$ [m]')
+                ax1.set_ylabel(r'$y$ [m]')
+                ax1.axis('equal')
+                _format_axes(ax1)
+                ax2 = fig.add_subplot(122)
+                ax2.set_xlabel(r'ln$(t/t_s)$')
+                ax2.set_ylabel(r'$\bar{Q}_b$')
+                _format_axes(ax2)
+
+                # Plot curves for requested boreholes
+                for i, borehole in enumerate(self.solver.boreholes):
+                    if i in iBoreholes:
+                        # Draw heat extraction rate
+                        line = ax2.plot(lntts, Q_t[iBoreholes.index(i)][n])
+                        color = line[-1]._color
+                        # Draw colored marker for borehole position
+                        if showTilt:
+                            ax1.plot(
+                                [borehole.x, borehole.x + borehole.H*np.sin(borehole.tilt)*np.cos(borehole.orientation)],
+                                [borehole.y, borehole.y + borehole.H*np.sin(borehole.tilt)*np.sin(borehole.orientation)],
+                                 linestyle='--',
+                                 marker='None',
+                                 color=color)
+                        ax1.plot(borehole.x,
+                                 borehole.y,
+                                 linestyle='None',
+                                 marker='o',
+                                 color=color)
+                    else:
+                        # Draw black marker for borehole position
+                        if showTilt:
+                            ax1.plot(
+                                [borehole.x, borehole.x + borehole.H*np.sin(borehole.tilt)*np.cos(borehole.orientation)],
+                                [borehole.y, borehole.y + borehole.H*np.sin(borehole.tilt)*np.sin(borehole.orientation)],
+                                 linestyle='--',
+                                 marker='None',
+                                 color='k')
+                        ax1.plot(borehole.x,
+                                 borehole.y,
+                                 linestyle='None',
+                                 marker='o',
+                                 color='k')
+
+                # Adjust figure to window
+                plt.tight_layout()
+            
         return fig
 
     def visualize_heat_extraction_rate_profiles(
-            self, time=None, iBoreholes=None, showTilt=True):
+            self, time=None, iBoreholes=None, showTilt=True, which=None):
         """
         Plot the heat extraction rate profiles at chosen time.
 
@@ -412,6 +536,11 @@ def visualize_heat_extraction_rate_profiles(
         showTilt : bool
             Set to True to show borehole inclination.
             Default is True
+        which : list of int, optional
+            Indices i of the diagonal variable mass flow rate g-functions for
+            which to plot heat extraction rates.
+            If None, all diagonal g-functions are plotted.
+            Default is None.
 
         Returns
         -------
@@ -425,58 +554,120 @@ def visualize_heat_extraction_rate_profiles(
         # Import heat extraction rate profiles
         z, Q_b = self._heat_extraction_rate_profiles(time, iBoreholes)
 
-        # Configure figure and axes
-        fig = _initialize_figure()
-        ax1 = fig.add_subplot(121)
-        ax1.set_xlabel(r'$x$ [m]')
-        ax1.set_ylabel(r'$y$ [m]')
-        ax1.axis('equal')
-        _format_axes(ax1)
-        ax2 = fig.add_subplot(122)
-        ax2.set_xlabel(r'$Q_b$')
-        ax2.set_ylabel(r'$z$ [m]')
-        ax2.invert_yaxis()
-        _format_axes(ax2)
-
-        # Plot curves for requested boreholes
-        for i, borehole in enumerate(self.solver.boreholes):
-            if i in iBoreholes:
-                # Draw heat extraction rate profile
-                line = ax2.plot(
-                    Q_b[iBoreholes.index(i)], z[iBoreholes.index(i)])
-                color = line[-1]._color
-                # Draw colored marker for borehole position
-                if showTilt:
-                    ax1.plot(
-                        [borehole.x, borehole.x + borehole.H*np.sin(borehole.tilt)*np.cos(borehole.orientation)],
-                        [borehole.y, borehole.y + borehole.H*np.sin(borehole.tilt)*np.sin(borehole.orientation)],
-                         linestyle='--',
-                         marker='None',
-                         color=color)
-                ax1.plot(borehole.x,
-                         borehole.y,
-                         linestyle='None',
-                         marker='o',
-                         color=color)
-            else:
-                # Draw black marker for borehole position
-                if showTilt:
-                    ax1.plot(
-                        [borehole.x, borehole.x + borehole.H*np.sin(borehole.tilt)*np.cos(borehole.orientation)],
-                        [borehole.y, borehole.y + borehole.H*np.sin(borehole.tilt)*np.sin(borehole.orientation)],
-                         linestyle='--',
-                         marker='None',
-                         color='k')
-                ax1.plot(borehole.x,
-                         borehole.y,
-                         linestyle='None',
-                         marker='o',
-                         color='k')
+        if self.solver.nMassFlow == 0:
+            # Configure figure and axes
+            fig = _initialize_figure()
+            ax1 = fig.add_subplot(121)
+            ax1.set_xlabel(r'$x$ [m]')
+            ax1.set_ylabel(r'$y$ [m]')
+            ax1.axis('equal')
+            _format_axes(ax1)
+            ax2 = fig.add_subplot(122)
+            ax2.set_xlabel(r'$Q_b$')
+            ax2.set_ylabel(r'$z$ [m]')
+            ax2.invert_yaxis()
+            _format_axes(ax2)
+
+            # Plot curves for requested boreholes
+            for i, borehole in enumerate(self.solver.boreholes):
+                if i in iBoreholes:
+                    # Draw heat extraction rate profile
+                    line = ax2.plot(
+                        Q_b[iBoreholes.index(i)], z[iBoreholes.index(i)])
+                    color = line[-1]._color
+                    # Draw colored marker for borehole position
+                    if showTilt:
+                        ax1.plot(
+                            [borehole.x, borehole.x + borehole.H*np.sin(borehole.tilt)*np.cos(borehole.orientation)],
+                            [borehole.y, borehole.y + borehole.H*np.sin(borehole.tilt)*np.sin(borehole.orientation)],
+                             linestyle='--',
+                             marker='None',
+                             color=color)
+                    ax1.plot(borehole.x,
+                             borehole.y,
+                             linestyle='None',
+                             marker='o',
+                             color=color)
+                else:
+                    # Draw black marker for borehole position
+                    if showTilt:
+                        ax1.plot(
+                            [borehole.x, borehole.x + borehole.H*np.sin(borehole.tilt)*np.cos(borehole.orientation)],
+                            [borehole.y, borehole.y + borehole.H*np.sin(borehole.tilt)*np.sin(borehole.orientation)],
+                             linestyle='--',
+                             marker='None',
+                             color='k')
+                    ax1.plot(borehole.x,
+                             borehole.y,
+                             linestyle='None',
+                             marker='o',
+                             color='k')
+
+            # Adjust figure to window
+            plt.tight_layout()
+        else:
+            m_flow = self.solver.m_flow
+            if which is None:
+                which = [n for n in range(self.solver.nMassFlow)]
+            for n in which:
+                # Configure figure and axes
+                fig = _initialize_figure()
+                fig.suptitle(
+                    f'Heat extraction rate profiles for m_flow={m_flow[n]} kg/s')
+                ax1 = fig.add_subplot(121)
+                ax1.set_xlabel(r'$x$ [m]')
+                ax1.set_ylabel(r'$y$ [m]')
+                ax1.axis('equal')
+                _format_axes(ax1)
+                ax2 = fig.add_subplot(122)
+                ax2.set_xlabel(r'$Q_b$')
+                ax2.set_ylabel(r'$z$ [m]')
+                ax2.invert_yaxis()
+                _format_axes(ax2)
+
+                # Plot curves for requested boreholes
+                for i, borehole in enumerate(self.solver.boreholes):
+                    if i in iBoreholes:
+                        # Draw heat extraction rate profile
+                        line = ax2.plot(
+                            Q_b[iBoreholes.index(i)][n],
+                            z[iBoreholes.index(i)])
+                        color = line[-1]._color
+                        # Draw colored marker for borehole position
+                        if showTilt:
+                            ax1.plot(
+                                [borehole.x, borehole.x + borehole.H*np.sin(borehole.tilt)*np.cos(borehole.orientation)],
+                                [borehole.y, borehole.y + borehole.H*np.sin(borehole.tilt)*np.sin(borehole.orientation)],
+                                 linestyle='--',
+                                 marker='None',
+                                 color=color)
+                        ax1.plot(borehole.x,
+                                 borehole.y,
+                                 linestyle='None',
+                                 marker='o',
+                                 color=color)
+                    else:
+                        # Draw black marker for borehole position
+                        if showTilt:
+                            ax1.plot(
+                                [borehole.x, borehole.x + borehole.H*np.sin(borehole.tilt)*np.cos(borehole.orientation)],
+                                [borehole.y, borehole.y + borehole.H*np.sin(borehole.tilt)*np.sin(borehole.orientation)],
+                                 linestyle='--',
+                                 marker='None',
+                                 color='k')
+                        ax1.plot(borehole.x,
+                                 borehole.y,
+                                 linestyle='None',
+                                 marker='o',
+                                 color='k')
+
+                # Adjust figure to window
+                plt.tight_layout()
 
-        plt.tight_layout()
         return fig
 
-    def visualize_temperatures(self, iBoreholes=None, showTilt=True):
+    def visualize_temperatures(
+            self, iBoreholes=None, showTilt=True, which=None):
         """
         Plot the time-variation of the average borehole wall temperatures.
 
@@ -489,6 +680,11 @@ def visualize_temperatures(self, iBoreholes=None, showTilt=True):
         showTilt : bool
             Set to True to show borehole inclination.
             Default is True
+        which : list of int, optional
+            Indices i of the diagonal variable mass flow rate g-functions for
+            which to plot borehole wall temperatures.
+            If None, all diagonal g-functions are plotted.
+            Default is None.
 
         Returns
         -------
@@ -501,63 +697,118 @@ def visualize_temperatures(self, iBoreholes=None, showTilt=True):
             iBoreholes = range(len(self.solver.boreholes))
         # Import temperatures
         T_b = self._temperatures(iBoreholes)
-
-        # Configure figure and axes
-        fig = _initialize_figure()
-        ax1 = fig.add_subplot(121)
-        ax1.set_xlabel(r'$x$ [m]')
-        ax1.set_ylabel(r'$y$ [m]')
-        ax1.axis('equal')
-        _format_axes(ax1)
-        ax2 = fig.add_subplot(122)
-        ax2.set_xlabel(r'ln$(t/t_s)$')
-        ax2.set_ylabel(r'$\bar{T}_b$')
-        _format_axes(ax2)
-
         # Borefield characteristic time
         ts = np.mean([b.H for b in self.solver.boreholes])**2/(9.*self.alpha)
         # Dimensionless time (log)
         lntts = np.log(self.time/ts)
-        # Plot curves for requested boreholes
-        for i, borehole in enumerate(self.solver.boreholes):
-            if i in iBoreholes:
-                # Draw borehole wall temperature
-                line = ax2.plot(lntts, T_b[iBoreholes.index(i)])
-                color = line[-1]._color
-                # Draw colored marker for borehole position
-                if showTilt:
-                    ax1.plot(
-                        [borehole.x, borehole.x + borehole.H*np.sin(borehole.tilt)*np.cos(borehole.orientation)],
-                        [borehole.y, borehole.y + borehole.H*np.sin(borehole.tilt)*np.sin(borehole.orientation)],
-                         linestyle='--',
-                         marker='None',
-                         color=color)
-                ax1.plot(borehole.x,
-                         borehole.y,
-                         linestyle='None',
-                         marker='o',
-                         color=color)
-            else:
-                # Draw black marker for borehole position
-                if showTilt:
-                    ax1.plot(
-                        [borehole.x, borehole.x + borehole.H*np.sin(borehole.tilt)*np.cos(borehole.orientation)],
-                        [borehole.y, borehole.y + borehole.H*np.sin(borehole.tilt)*np.sin(borehole.orientation)],
-                         linestyle='--',
-                         marker='None',
-                         color='k')
-                ax1.plot(borehole.x,
-                         borehole.y,
-                         linestyle='None',
-                         marker='o',
-                         color='k')
 
-        # Adjust figure to window
-        plt.tight_layout()
+
+        if self.solver.nMassFlow == 0:
+            # Configure figure and axes
+            fig = _initialize_figure()
+            ax1 = fig.add_subplot(121)
+            ax1.set_xlabel(r'$x$ [m]')
+            ax1.set_ylabel(r'$y$ [m]')
+            ax1.axis('equal')
+            _format_axes(ax1)
+            ax2 = fig.add_subplot(122)
+            ax2.set_xlabel(r'ln$(t/t_s)$')
+            ax2.set_ylabel(r'$\bar{T}_b$')
+            _format_axes(ax2)
+            # Plot curves for requested boreholes
+            for i, borehole in enumerate(self.solver.boreholes):
+                if i in iBoreholes:
+                    # Draw borehole wall temperature
+                    line = ax2.plot(lntts, T_b[iBoreholes.index(i)])
+                    color = line[-1]._color
+                    # Draw colored marker for borehole position
+                    if showTilt:
+                        ax1.plot(
+                            [borehole.x, borehole.x + borehole.H*np.sin(borehole.tilt)*np.cos(borehole.orientation)],
+                            [borehole.y, borehole.y + borehole.H*np.sin(borehole.tilt)*np.sin(borehole.orientation)],
+                             linestyle='--',
+                             marker='None',
+                             color=color)
+                    ax1.plot(borehole.x,
+                             borehole.y,
+                             linestyle='None',
+                             marker='o',
+                             color=color)
+                else:
+                    # Draw black marker for borehole position
+                    if showTilt:
+                        ax1.plot(
+                            [borehole.x, borehole.x + borehole.H*np.sin(borehole.tilt)*np.cos(borehole.orientation)],
+                            [borehole.y, borehole.y + borehole.H*np.sin(borehole.tilt)*np.sin(borehole.orientation)],
+                             linestyle='--',
+                             marker='None',
+                             color='k')
+                    ax1.plot(borehole.x,
+                             borehole.y,
+                             linestyle='None',
+                             marker='o',
+                             color='k')
+
+            # Adjust figure to window
+            plt.tight_layout()
+        else:
+            m_flow = self.solver.m_flow
+            if which is None:
+                which = [n for n in range(self.solver.nMassFlow)]
+            for n in which:
+                # Configure figure and axes
+                fig = _initialize_figure()
+                fig.suptitle(
+                    f'Borehole wall temperatures for m_flow={m_flow[n]} kg/s')
+                ax1 = fig.add_subplot(121)
+                ax1.set_xlabel(r'$x$ [m]')
+                ax1.set_ylabel(r'$y$ [m]')
+                ax1.axis('equal')
+                _format_axes(ax1)
+                ax2 = fig.add_subplot(122)
+                ax2.set_xlabel(r'ln$(t/t_s)$')
+                ax2.set_ylabel(r'$\bar{T}_b$')
+                _format_axes(ax2)
+                # Plot curves for requested boreholes
+                for i, borehole in enumerate(self.solver.boreholes):
+                    if i in iBoreholes:
+                        # Draw borehole wall temperature
+                        line = ax2.plot(lntts, T_b[iBoreholes.index(i)][n])
+                        color = line[-1]._color
+                        # Draw colored marker for borehole position
+                        if showTilt:
+                            ax1.plot(
+                                [borehole.x, borehole.x + borehole.H*np.sin(borehole.tilt)*np.cos(borehole.orientation)],
+                                [borehole.y, borehole.y + borehole.H*np.sin(borehole.tilt)*np.sin(borehole.orientation)],
+                                 linestyle='--',
+                                 marker='None',
+                                 color=color)
+                        ax1.plot(borehole.x,
+                                 borehole.y,
+                                 linestyle='None',
+                                 marker='o',
+                                 color=color)
+                    else:
+                        # Draw black marker for borehole position
+                        if showTilt:
+                            ax1.plot(
+                                [borehole.x, borehole.x + borehole.H*np.sin(borehole.tilt)*np.cos(borehole.orientation)],
+                                [borehole.y, borehole.y + borehole.H*np.sin(borehole.tilt)*np.sin(borehole.orientation)],
+                                 linestyle='--',
+                                 marker='None',
+                                 color='k')
+                        ax1.plot(borehole.x,
+                                 borehole.y,
+                                 linestyle='None',
+                                 marker='o',
+                                 color='k')
+
+                # Adjust figure to window
+                plt.tight_layout()
         return fig
 
     def visualize_temperature_profiles(
-            self, time=None, iBoreholes=None, showTilt=True):
+            self, time=None, iBoreholes=None, showTilt=True, which=None):
         """
         Plot the borehole wall temperature profiles at chosen time.
 
@@ -574,6 +825,11 @@ def visualize_temperature_profiles(
         showTilt : bool
             Set to True to show borehole inclination.
             Default is True
+        which : list of int, optional
+            Indices i of the diagonal variable mass flow rate g-functions for
+            which to plot borehole wall temperatures.
+            If None, all diagonal g-functions are plotted.
+            Default is None.
 
         Returns
         -------
@@ -587,55 +843,114 @@ def visualize_temperature_profiles(
         # Import temperature profiles
         z, T_b = self._temperature_profiles(time, iBoreholes)
 
-        # Configure figure and axes
-        fig = _initialize_figure()
-        ax1 = fig.add_subplot(121)
-        ax1.set_xlabel(r'$x$ [m]')
-        ax1.set_ylabel(r'$y$ [m]')
-        ax1.axis('equal')
-        _format_axes(ax1)
-        ax2 = fig.add_subplot(122)
-        ax2.set_xlabel(r'$T_b$')
-        ax2.set_ylabel(r'$z$ [m]')
-        ax2.invert_yaxis()
-        _format_axes(ax2)
-
-        # Plot curves for requested boreholes
-        for i, borehole in enumerate(self.solver.boreholes):
-            if i in iBoreholes:
-                # Draw heat extraction rate profile
-                line = ax2.plot(
-                    T_b[iBoreholes.index(i)], z[iBoreholes.index(i)])
-                color = line[-1]._color
-                # Draw colored marker for borehole position
-                if showTilt:
-                    ax1.plot(
-                        [borehole.x, borehole.x + borehole.H*np.sin(borehole.tilt)*np.cos(borehole.orientation)],
-                        [borehole.y, borehole.y + borehole.H*np.sin(borehole.tilt)*np.sin(borehole.orientation)],
-                         linestyle='--',
-                         marker='None',
-                         color=color)
-                ax1.plot(borehole.x,
-                         borehole.y,
-                         linestyle='None',
-                         marker='o',
-                         color=color)
-            else:
-                # Draw black marker for borehole position
-                if showTilt:
-                    ax1.plot(
-                        [borehole.x, borehole.x + borehole.H*np.sin(borehole.tilt)*np.cos(borehole.orientation)],
-                        [borehole.y, borehole.y + borehole.H*np.sin(borehole.tilt)*np.sin(borehole.orientation)],
-                         linestyle='--',
-                         marker='None',
-                         color='k')
-                ax1.plot(borehole.x,
-                         borehole.y,
-                         linestyle='None',
-                         marker='o',
-                         color='k')
-
-        plt.tight_layout()
+        if self.solver.nMassFlow == 0:
+            # Configure figure and axes
+            fig = _initialize_figure()
+            ax1 = fig.add_subplot(121)
+            ax1.set_xlabel(r'$x$ [m]')
+            ax1.set_ylabel(r'$y$ [m]')
+            ax1.axis('equal')
+            _format_axes(ax1)
+            ax2 = fig.add_subplot(122)
+            ax2.set_xlabel(r'$T_b$')
+            ax2.set_ylabel(r'$z$ [m]')
+            ax2.invert_yaxis()
+            _format_axes(ax2)
+
+            # Plot curves for requested boreholes
+            for i, borehole in enumerate(self.solver.boreholes):
+                if i in iBoreholes:
+                    # Draw borehole wall temperature profile
+                    line = ax2.plot(
+                        T_b[iBoreholes.index(i)],
+                        z[iBoreholes.index(i)])
+                    color = line[-1]._color
+                    # Draw colored marker for borehole position
+                    if showTilt:
+                        ax1.plot(
+                            [borehole.x, borehole.x + borehole.H * np.sin(borehole.tilt) * np.cos(borehole.orientation)],
+                            [borehole.y, borehole.y + borehole.H * np.sin(borehole.tilt) * np.sin(borehole.orientation)],
+                             linestyle='--',
+                             marker='None',
+                             color=color)
+                    ax1.plot(borehole.x,
+                             borehole.y,
+                             linestyle='None',
+                             marker='o',
+                             color=color)
+                else:
+                    # Draw black marker for borehole position
+                    if showTilt:
+                        ax1.plot(
+                            [borehole.x, borehole.x + borehole.H * np.sin(borehole.tilt) * np.cos(borehole.orientation)],
+                            [borehole.y, borehole.y + borehole.H * np.sin(borehole.tilt) * np.sin(borehole.orientation)],
+                             linestyle='--',
+                             marker='None',
+                             color='k')
+                    ax1.plot(borehole.x,
+                             borehole.y,
+                             linestyle='None',
+                             marker='o',
+                             color='k')
+
+            plt.tight_layout()
+        else:
+            m_flow = self.solver.m_flow
+            if which is None:
+                which = [n for n in range(self.solver.nMassFlow)]
+            for n in which:
+                # Configure figure and axes
+                fig = _initialize_figure()
+                fig.suptitle(
+                f'Borehole wall temperature profiles for m_flow={m_flow[n]} kg/s')
+                ax1 = fig.add_subplot(121)
+                ax1.set_xlabel(r'$x$ [m]')
+                ax1.set_ylabel(r'$y$ [m]')
+                ax1.axis('equal')
+                _format_axes(ax1)
+                ax2 = fig.add_subplot(122)
+                ax2.set_xlabel(r'$T_b$')
+                ax2.set_ylabel(r'$z$ [m]')
+                ax2.invert_yaxis()
+                _format_axes(ax2)
+
+                # Plot curves for requested boreholes
+                for i, borehole in enumerate(self.solver.boreholes):
+                    if i in iBoreholes:
+                        # Draw borehole wall temperature profile
+                        line = ax2.plot(
+                            T_b[iBoreholes.index(i)][n],
+                            z[iBoreholes.index(i)])
+                        color = line[-1]._color
+                        # Draw colored marker for borehole position
+                        if showTilt:
+                            ax1.plot(
+                                [borehole.x, borehole.x + borehole.H * np.sin(borehole.tilt) * np.cos(borehole.orientation)],
+                                [borehole.y, borehole.y + borehole.H * np.sin(borehole.tilt) * np.sin(borehole.orientation)],
+                                 linestyle='--',
+                                 marker='None',
+                                 color=color)
+                        ax1.plot(borehole.x,
+                                 borehole.y,
+                                 linestyle='None',
+                                 marker='o',
+                                 color=color)
+                    else:
+                        # Draw black marker for borehole position
+                        if showTilt:
+                            ax1.plot(
+                                [borehole.x, borehole.x + borehole.H * np.sin(borehole.tilt) * np.cos(borehole.orientation)],
+                                [borehole.y, borehole.y + borehole.H * np.sin(borehole.tilt) * np.sin(borehole.orientation)],
+                                 linestyle='--',
+                                 marker='None',
+                                 color='k')
+                        ax1.plot(borehole.x,
+                                 borehole.y,
+                                 linestyle='None',
+                                 marker='o',
+                                 color='k')
+
+                plt.tight_layout()
         return fig
 
     def _heat_extraction_rates(self, iBoreholes):
@@ -667,9 +982,18 @@ def _heat_extraction_rates(self, iBoreholes):
                 i0 = self.solver._i0Segments[i]
                 i1 = self.solver._i1Segments[i]
                 segment_ratios = self.solver.segment_ratios[i]
-                Q_t.append(
-                    np.sum(self.solver.Q_b[i0:i1,:]*segment_ratios[:,np.newaxis],
-                           axis=0))
+                if self.solver.nMassFlow == 0:
+                    Q_t.append(
+                        np.sum(
+                            self.solver.Q_b[i0:i1, :]
+                            * segment_ratios[:, np.newaxis],
+                            axis=0))
+                else:
+                    Q_t.append(
+                        np.sum(
+                            self.solver.Q_b[:, i0:i1, :]
+                            * segment_ratios[:, np.newaxis],
+                            axis=1))
         return Q_t
 
     def _heat_extraction_rate_profiles(self, time, iBoreholes):
@@ -712,13 +1036,26 @@ def _heat_extraction_rate_profiles(self, time, iBoreholes):
                 if time is None:
                     # If time is None, heat extraction rates are extracted at
                     # the last time step.
-                    Q_bi = self.solver.Q_b[i0:i1,-1].flatten()
+                    if self.solver.nMassFlow == 0:
+                        Q_bi = self.solver.Q_b[i0:i1, -1]
+                    else:
+                        Q_bi = self.solver.Q_b[:, i0:i1, -1]
                 else:
                     # Otherwise, heat extraction rates are interpolated.
-                    Q_bi = interp1d(self.time, self.solver.Q_b[i0:i1,:],
-                                    kind='linear',
-                                    copy=False,
-                                    axis=1)(time).flatten()
+                    if self.solver.nMassFlow == 0:
+                        Q_bi = interp1d(
+                            self.time,
+                            self.solver.Q_b[i0:i1,:],
+                            kind='linear',
+                            copy=False,
+                            axis=1)(time)
+                    else:
+                        Q_bi = interp1d(
+                            self.time,
+                            self.solver.Q_b[i0:i1,:],
+                            kind='linear',
+                            copy=False,
+                            axis=2)(time)
                 if self.solver.nBoreSegments[i] > 1:
                     # Borehole length ratio at the mid-depth of each segment
                     segment_ratios = self.solver.segment_ratios[i]
@@ -734,7 +1071,10 @@ def _heat_extraction_rate_profiles(self, time, iBoreholes):
                     z.append(
                         np.array([self.solver.boreholes[i].D,
                                   self.solver.boreholes[i].D + self.solver.boreholes[i].H]))
-                    Q_b.append(np.array(2*[np.asscalar(Q_bi)]))
+                    if self.solver.nMassFlow == 0:
+                        Q_b.append(np.repeat(Q_bi, 2, axis=0))
+                    else:
+                        Q_b.append(np.repeat(Q_bi, 2, axis=1))
         return z, Q_b
 
     def _temperatures(self, iBoreholes):
@@ -766,9 +1106,18 @@ def _temperatures(self, iBoreholes):
                 i0 = self.solver._i0Segments[i]
                 i1 = self.solver._i1Segments[i]
                 segment_ratios = self.solver.segment_ratios[i]
-                T_b.append(
-                    np.sum(self.solver.T_b[i0:i1,:]*segment_ratios[:,np.newaxis],
-                           axis=0))
+                if self.solver.nMassFlow == 0:
+                    T_b.append(
+                        np.sum(
+                            self.solver.T_b[i0:i1, :]
+                            * segment_ratios[:, np.newaxis],
+                            axis=0))
+                else:
+                    T_b.append(
+                        np.sum(
+                            self.solver.T_b[:, i0:i1, :]
+                            * segment_ratios[:, np.newaxis],
+                            axis=1))
         return T_b
 
     def _temperature_profiles(self, time, iBoreholes):
@@ -821,14 +1170,31 @@ def _temperature_profiles(self, time, iBoreholes):
                 if time is None:
                     # If time is None, temperatures are extracted at the last
                     # time step.
-                    T_bi = self.solver.T_b[i0:i1,-1].flatten()
+                    if self.solver.nMassFlow == 0:
+                        T_bi = self.solver.T_b[i0:i1, -1]
+                    else:
+                        T_bi = self.solver.T_b[:, i0:i1, -1]
                 else:
                     # Otherwise, temperatures are interpolated.
                     T_bi = interp1d(self.time,
-                                    self.solver.T_b[i0:i1,:],
+                                    self.solver.T_b[i0:i1, :],
                                     kind='linear',
                                     copy=False,
-                                    axis=1)(time).flatten()
+                                    axis=1)(time)
+                    if self.solver.nMassFlow == 0:
+                        T_bi = interp1d(
+                            self.time,
+                            self.solver.T_b[i0:i1,:],
+                            kind='linear',
+                            copy=False,
+                            axis=1)(time)
+                    else:
+                        T_bi = interp1d(
+                            self.time,
+                            self.solver.T_b[i0:i1,:],
+                            kind='linear',
+                            copy=False,
+                            axis=2)(time)
                 if self.solver.nBoreSegments[i] > 1:
                     # Borehole length ratio at the mid-depth of each segment
 
@@ -845,7 +1211,10 @@ def _temperature_profiles(self, time, iBoreholes):
                     z.append(
                         np.array([self.solver.boreholes[i].D,
                                   self.solver.boreholes[i].D + self.solver.boreholes[i].H]))
-                    T_b.append(np.array(2*[np.asscalar(T_bi)]))
+                    if self.solver.nMassFlow == 0:
+                        T_b.append(np.repeat(T_bi, 2, axis=0))
+                    else:
+                        T_b.append(np.repeat(T_bi, 2, axis=1))
         return z, T_b
 
     def _format_inputs(self, boreholes_or_network):
@@ -865,6 +1234,15 @@ def _format_inputs(self, boreholes_or_network):
             # use 'MIFT'
             if self.boundary_condition is None:
                 self.boundary_condition = 'MIFT'
+            # Extract mass flow rate from Network object if provided in the object
+            # and none of m_flow_borehole and m_flow_network are provided
+            if self.m_flow_borehole is None and self.m_flow_network is None:
+                if type(self.network.m_flow_network) is float:
+                     self.m_flow_network = self.network.m_flow_network
+                elif type(self.network.m_flow_network) is np.ndarray:
+                     self.m_flow_borehole = self.network.m_flow_network
+            if self.cp_f is None and type(self.network.cp_f) is float:
+                self.cp_f = self.network.cp_f
         else:
             self.network = None
             self.boreholes = boreholes_or_network
@@ -881,11 +1259,10 @@ def _format_inputs(self, boreholes_or_network):
                     "parallel-connected boreholes. Calculations will use the "
                     "'similarities' solver instead.")
                 self.method = 'similarities'
-            elif not (
-                    type(self.network.m_flow_network) is float or (
-                        type(self.network.m_flow_network) is np.ndarray and \
-                            np.allclose(self.network.m_flow_network,
-                                        self.network.m_flow_network[0]))):
+            elif not (self.m_flow_borehole is None
+                      or np.allclose(
+                          self.m_flow_borehole,
+                          self.m_flow_borehole[0])):
                 warnings.warn(
                     "\nSolver 'equivalent' is only valid for equal mass flow "
                     "rates into the boreholes. Calculations will use the "
@@ -918,9 +1295,19 @@ def _check_inputs(self):
             "There are duplicate boreholes in the borefield."
         assert (self.network is None and not self.boundary_condition=='MIFT') or isinstance(self.network, Network), \
             "The network is not a valid 'Network' object."
-        assert self.network is None or (self.network.m_flow_network is not None and self.network.cp_f is not None), \
-            "The mass flow rate 'm_flow_network' and heat capacity 'cp_f' must " \
-            "be provided at the instanciation of the 'Network' object."
+        if self.boundary_condition == 'MIFT':
+            assert not (self.m_flow_network is None and self.m_flow_borehole is None), \
+                "The mass flow rate 'm_flow_borehole' or 'm_flow_network' must " \
+                "be provided when using the 'MIFT' boundary condition."
+            assert not (self.m_flow_network is not None and self.m_flow_borehole is not None), \
+                "Only one of 'm_flow_borehole' or 'm_flow_network' can " \
+                "be provided when using the 'MIFT' boundary condition."
+            assert not self.cp_f is None, \
+                "The heat capacity 'cp_f' must " \
+                "be provided when using the 'MIFT' boundary condition."
+            assert not (type(self.m_flow_borehole) is np.ndarray and self.m_flow_borehole.ndim == 2 and not np.size(self.m_flow_borehole, axis=1)==self.network.nInlets), \
+                "The number of mass flow rates in 'm_flow_borehole' must " \
+                "correspond to the number of circuits in the network."
         assert type(self.time) is np.ndarray or isinstance(self.time, (np.floating, float)) or self.time is None, \
             "Time should be a float or an array."
         assert isinstance(self.alpha, (np.floating, float)), \
@@ -1406,6 +1793,27 @@ class _BaseSolver(object):
         (of type int) as an argument, or an array of size (nSegments[i],) when
         provided with an element of nSegments (of type list).
         Default is :func:`utilities.segment_ratios`.
+    m_flow_borehole : (nInlets,) array or (nMassFlow, nInlets,) array, optional
+        Fluid mass flow rate into each circuit of the network. If a
+        (nMassFlow, nInlets,) array is supplied, the
+        (nMassFlow, nMassFlow,) variable mass flow rate g-functions
+        will be evaluated using the method of Cimmino (2024)
+        [#gFunction-CimBer2024]_. Only required for the 'MIFT' boundary
+         condition. Only one of 'm_flow_borehole' and 'm_flow_network' can be
+         provided.
+        Default is None.
+    m_flow_network : float or (nMassFlow,) array, optional
+        Fluid mass flow rate into the network of boreholes. If an array
+        is supplied, the (nMassFlow, nMassFlow,) variable mass flow
+        rate g-functions will be evaluated using the method of Cimmino
+        (2024) [#gFunction-CimBer2024]_. Only required for the 'MIFT' boundary
+         condition. Only one of 'm_flow_borehole' and 'm_flow_network' can be
+         provided.
+        Default is None.
+    cp_f : float, optional
+        Fluid specific isobaric heat capacity (in J/kg.degC). Only required
+        for the 'MIFT' boundary condition.
+        Default is None.
     approximate_FLS : bool, optional
         Set to true to use the approximation of the FLS solution of Cimmino
         (2021). This approximation does not require the numerical evaluation of
@@ -1446,6 +1854,7 @@ class _BaseSolver(object):
 
     """
     def __init__(self, boreholes, network, time, boundary_condition,
+                 m_flow_borehole=None, m_flow_network=None, cp_f=None,
                  nSegments=8, segment_ratios=utilities.segment_ratios,
                  approximate_FLS=False, mQuad=11, nFLS=10,
                  linear_threshold=None, disp=False, profiles=False,
@@ -1485,6 +1894,20 @@ def __init__(self, boreholes, network, time, boundary_condition,
                             for i in range(nBoreholes)]
         self._i1Segments = [sum(self.nBoreSegments[0:(i + 1)])
                             for i in range(nBoreholes)]
+        self.nMassFlow = 0
+        self.m_flow_borehole = m_flow_borehole
+        if self.m_flow_borehole is not None:
+            if not self.m_flow_borehole.ndim == 1:
+                self.nMassFlow = np.size(self.m_flow_borehole, axis=0)
+            self.m_flow_borehole = np.atleast_2d(self.m_flow_borehole)
+            self.m_flow = self.m_flow_borehole
+        self.m_flow_network = m_flow_network
+        if self.m_flow_network is not None:
+            if not isinstance(self.m_flow_network, (np.floating, float)):
+                self.nMassFlow = len(self.m_flow_network)
+            self.m_flow_network = np.atleast_1d(self.m_flow_network)
+            self.m_flow = self.m_flow_network
+        self.cp_f = cp_f
         self.approximate_FLS = approximate_FLS
         self.mQuad = mQuad
         self.nFLS = nFLS
@@ -1555,17 +1978,6 @@ def solve(self, time, alpha):
         else:
             time_long = self.time
         nt_long = len(time_long)
-        # Initialize g-function
-        gFunc = np.zeros(nt)
-        # Initialize segment heat extraction rates
-        if self.boundary_condition == 'UHTR':
-            Q_b = 1
-        else:
-            Q_b = np.zeros((self.nSources, nt), dtype=self.dtype)
-        if self.boundary_condition == 'UBWT':
-            T_b = np.zeros(nt, dtype=self.dtype)
-        else:
-            T_b = np.zeros((self.nSources, nt), dtype=self.dtype)
         # Calculate segment to segment thermal response factors
         h_ij = self.thermal_response_factors(time_long, alpha, kind=self.kind)
         # Segment lengths
@@ -1578,10 +1990,17 @@ def solve(self, time, alpha):
         # Initialize chrono
         tic = perf_counter()
 
-        # Build and solve the system of equations at all times
-        p0 = max(0, p_long-1)
-        for p in range(nt_long):
-            if self.boundary_condition == 'UHTR':
+        if self.boundary_condition == 'UHTR':
+            # Initialize g-function
+            gFunc = np.zeros(nt)
+            # Initialize segment heat extraction rates
+            Q_b = 1
+            # Initialize borehole wall temperatures
+            T_b = np.zeros((self.nSources, nt), dtype=self.dtype)
+
+            # Build and solve the system of equations at all times
+            p0 = max(0, p_long-1)
+            for p in range(nt_long):
                 # Evaluate the g-function with uniform heat extraction along
                 # boreholes
 
@@ -1593,7 +2012,22 @@ def solve(self, time, alpha):
                 # The g-function is the average of all borehole wall
                 # temperatures
                 gFunc[p+p0] = np.sum(T_b[:,p+p0]*H_b)/H_tot
-            else:
+
+            # Linearize g-function for times under threshold
+            if p_long > 0:
+                gFunc[:p_long] = gFunc[p_long-1] * self.time[:p_long] / time_threshold
+                T_b[:,:p_long] = T_b[:,p_long-1:p_long] * self.time[:p_long] / time_threshold
+
+        elif self.boundary_condition == 'UBWT':
+            # Initialize g-function
+            gFunc = np.zeros(nt)
+            # Initialize segment heat extraction rates
+            Q_b = np.zeros((self.nSources, nt), dtype=self.dtype)
+            T_b = np.zeros(nt, dtype=self.dtype)
+
+            # Build and solve the system of equations at all times
+            p0 = max(0, p_long-1)
+            for p in range(nt_long):
                 # Current thermal response factor matrix
                 if p > 0:
                     dt = time_long[p] - time_long[p-1]
@@ -1609,32 +2043,78 @@ def solve(self, time, alpha):
                 T_b0 = self.temporal_superposition(
                     h_ij.y[:,:,1:], Q_reconstructed)
 
-                if self.boundary_condition == 'UBWT':
-                    # Evaluate the g-function with uniform borehole wall
-                    # temperature
-                    # ---------------------------------------------------------
-                    # Build a system of equation [A]*[X] = [B] for the
-                    # evaluation of the g-function. [A] is a coefficient
-                    # matrix, [X] = [Q_b,T_b] is a state space vector of the
-                    # borehole heat extraction rates and borehole wall
-                    # temperature (equal for all segments), [B] is a
-                    # coefficient vector.
-                    #
-                    # Spatial superposition: [T_b] = [T_b0] + [h_ij_dt]*[Q_b]
-                    # Energy conservation: sum([Q_b*Hb]) = sum([Hb])
-                    # ---------------------------------------------------------
-                    A = np.block([[h_dt, -np.ones((self.nSources, 1),
-                                                  dtype=self.dtype)],
-                                  [H_b, 0.]])
-                    B = np.hstack((-T_b0, H_tot))
-                    # Solve the system of equations
-                    X = np.linalg.solve(A, B)
-                    # Store calculated heat extraction rates
-                    Q_b[:,p+p0] = X[0:self.nSources]
-                    # The borehole wall temperatures are equal for all segments
-                    T_b[p+p0] = X[-1]
-                    gFunc[p+p0] = T_b[p+p0]
-                elif self.boundary_condition == 'MIFT':
+                # Evaluate the g-function with uniform borehole wall
+                # temperature
+                # ---------------------------------------------------------
+                # Build a system of equation [A]*[X] = [B] for the
+                # evaluation of the g-function. [A] is a coefficient
+                # matrix, [X] = [Q_b,T_b] is a state space vector of the
+                # borehole heat extraction rates and borehole wall
+                # temperature (equal for all segments), [B] is a
+                # coefficient vector.
+                #
+                # Spatial superposition: [T_b] = [T_b0] + [h_ij_dt]*[Q_b]
+                # Energy conservation: sum([Q_b*Hb]) = sum([Hb])
+                # ---------------------------------------------------------
+                A = np.block([[h_dt, -np.ones((self.nSources, 1),
+                                              dtype=self.dtype)],
+                              [H_b, 0.]])
+                B = np.hstack((-T_b0, H_tot))
+                # Solve the system of equations
+                X = np.linalg.solve(A, B)
+                # Store calculated heat extraction rates
+                Q_b[:,p+p0] = X[0:self.nSources]
+                # The borehole wall temperatures are equal for all segments
+                T_b[p+p0] = X[-1]
+                gFunc[p+p0] = T_b[p+p0]
+
+            # Linearize g-function for times under threshold
+            if p_long > 0:
+                gFunc[:p_long] = gFunc[p_long-1] * self.time[:p_long] / time_threshold
+                Q_b[:,:p_long] = 1 + (Q_b[:,p_long-1:p_long] - 1) * self.time[:p_long] / time_threshold
+                T_b[:p_long] = T_b[p_long-1] * self.time[:p_long] / time_threshold
+
+        elif self.boundary_condition == 'MIFT':
+            if self.nMassFlow == 0:
+                # Initialize g-function
+                gFunc = np.zeros((1, 1, nt))
+                # Initialize segment heat extraction rates
+                Q_b = np.zeros((1, self.nSources, nt), dtype=self.dtype)
+                T_b = np.zeros((1, self.nSources, nt), dtype=self.dtype)
+            else:
+                # Initialize g-function
+                gFunc = np.zeros((self.nMassFlow, self.nMassFlow, nt))
+                # Initialize segment heat extraction rates
+                Q_b = np.zeros(
+                    (self.nMassFlow, self.nSources, nt), dtype=self.dtype)
+                T_b = np.zeros(
+                    (self.nMassFlow, self.nSources, nt), dtype=self.dtype)
+
+            for j in range(np.maximum(self.nMassFlow, 1)):
+                # Build and solve the system of equations at all times
+                p0 = max(0, p_long-1)
+                a_in_j, a_b_j = self.network.coefficients_borehole_heat_extraction_rate(
+                        self.m_flow[j],
+                        self.cp_f,
+                        self.nBoreSegments,
+                        segment_ratios=self.segment_ratios)
+                k_s = self.network.p[0].k_s
+                for p in range(nt_long):
+                    # Current thermal response factor matrix
+                    if p > 0:
+                        dt = time_long[p] - time_long[p-1]
+                    else:
+                        dt = time_long[p]
+                    # Thermal response factors evaluated at t=dt
+                    h_dt = h_ij(dt)
+                    # Reconstructed load history
+                    Q_reconstructed = self.load_history_reconstruction(
+                        time_long[0:p+1], Q_b[j,:,p0:p+p0+1])
+                    # Borehole wall temperature for zero heat extraction at
+                    # current step
+                    T_b0 = self.temporal_superposition(
+                        h_ij.y[:,:,1:], Q_reconstructed)
+
                     # Evaluate the g-function with mixed inlet fluid
                     # temperatures
                     # ---------------------------------------------------------
@@ -1650,19 +2130,13 @@ def solve(self, time, alpha):
                     # [Q_{b,i}] = [a_in]*[T_{f,in}] + [a_{b,i}]*[T_{b,i}]
                     # Energy conservation: sum([Q_b*H_b]) = sum([H_b])
                     # ---------------------------------------------------------
-                    a_in, a_b = self.network.coefficients_borehole_heat_extraction_rate(
-                            self.network.m_flow_network,
-                            self.network.cp_f,
-                            self.nBoreSegments,
-                            segment_ratios=self.segment_ratios)
-                    k_s = self.network.p[0].k_s
                     A = np.block(
                         [[h_dt,
                           -np.eye(self.nSources, dtype=self.dtype),
                           np.zeros((self.nSources, 1), dtype=self.dtype)],
                          [np.eye(self.nSources, dtype=self.dtype),
-                          a_b/(2.0*pi*k_s*np.atleast_2d(Hb_individual).T),
-                          a_in/(2.0*pi*k_s*np.atleast_2d(Hb_individual).T)],
+                          a_b_j/(2.0*pi*k_s*np.atleast_2d(Hb_individual).T),
+                          a_in_j/(2.0*pi*k_s*np.atleast_2d(Hb_individual).T)],
                          [H_b, np.zeros(self.nSources + 1, dtype=self.dtype)]])
                     B = np.hstack(
                         (-T_b0,
@@ -1671,34 +2145,55 @@ def solve(self, time, alpha):
                     # Solve the system of equations
                     X = np.linalg.solve(A, B)
                     # Store calculated heat extraction rates
-                    Q_b[:,p+p0] = X[0:self.nSources]
-                    T_b[:,p+p0] = X[self.nSources:2*self.nSources]
+                    Q_b[j,:,p+p0] = X[0:self.nSources]
+                    T_b[j,:,p+p0] = X[self.nSources:2*self.nSources]
+                    # Inlet fluid temperature
                     T_f_in = X[-1]
                     # The gFunction is equal to the effective borehole wall
                     # temperature
                     # Outlet fluid temperature
-                    T_f_out = T_f_in - 2*pi*self.network.p[0].k_s*H_tot / (
-                        np.sum(
-                            np.abs(self.network.m_flow_network)
-                            * self.network.cp_f))
+                    T_f_out = T_f_in - 2 * pi * k_s * H_tot / (
+                        np.sum(np.abs(self.m_flow[j]) * self.cp_f))
                     # Average fluid temperature
                     T_f = 0.5*(T_f_in + T_f_out)
                     # Borefield thermal resistance
                     R_field = network_thermal_resistance(
-                        self.network, self.network.m_flow_network,
-                        self.network.cp_f)
+                        self.network, self.m_flow[j], self.cp_f)
                     # Effective borehole wall temperature
-                    T_b_eff = T_f - 2*pi*self.network.p[0].k_s*R_field
-                    gFunc[p+p0] = T_b_eff
-        # Linearize g-function for times under threshold
-        if p_long > 0:
-            gFunc[:p_long] = gFunc[p_long-1] * self.time[:p_long] / time_threshold
-            if not self.boundary_condition == 'UHTR':
-                Q_b[:,:p_long] = 1 + (Q_b[:,p_long-1:p_long] - 1) * self.time[:p_long] / time_threshold
-            if self.boundary_condition == 'UBWT':
-                T_b[:p_long] = T_b[p_long-1] * self.time[:p_long] / time_threshold
-            else:
-                T_b[:,:p_long] = T_b[:,p_long-1:p_long] * self.time[:p_long] / time_threshold
+                    T_b_eff = T_f - 2 * pi * k_s * R_field
+                    gFunc[j,j,p+p0] = T_b_eff
+
+            for i in range(np.maximum(self.nMassFlow, 1)):
+                for j in range(np.maximum(self.nMassFlow, 1)):
+                    if not i == j:
+                        # Inlet fluid temperature
+                        a_in, a_b = self.network.coefficients_network_heat_extraction_rate(
+                                self.m_flow[i],
+                                self.cp_f,
+                                self.nBoreSegments,
+                                segment_ratios=self.segment_ratios)
+                        T_f_in = (-2 * pi * k_s * H_tot - a_b @ T_b[j,:,p0:]) / a_in
+                        # The gFunction is equal to the effective borehole wall
+                        # temperature
+                        # Outlet fluid temperature
+                        T_f_out = T_f_in - 2 * pi * k_s * H_tot / np.sum(np.abs(self.m_flow[i]) * self.cp_f)
+                        # Borefield thermal resistance
+                        R_field = network_thermal_resistance(
+                            self.network, self.m_flow[i], self.cp_f)
+                        # Effective borehole wall temperature
+                        T_b_eff = 0.5 * (T_f_in + T_f_out) - 2 * pi * k_s * R_field
+                        gFunc[i,j,p0:] = T_b_eff
+
+            # Linearize g-function for times under threshold
+            if p_long > 0:
+                gFunc[:,:,:p_long] = gFunc[:,:,p_long-1] * self.time[:p_long] / time_threshold
+                Q_b[:,:,:p_long] = 1 + (Q_b[:,:,p_long-1:p_long] - 1) * self.time[:p_long] / time_threshold
+                T_b[:,:,:p_long] = T_b[:,:,p_long-1:p_long] * self.time[:p_long] / time_threshold
+            if self.nMassFlow == 0:
+                gFunc = gFunc[0,0,:]
+                Q_b = Q_b[0,:,:]
+                T_b = T_b[0,:,:]
+
         # Store temperature and heat extraction rate profiles
         if self.profiles:
             self.Q_b = Q_b
@@ -1834,9 +2329,19 @@ def _check_inputs(self):
             "objects."
         assert self.network is None or isinstance(self.network, Network), \
             "The network is not a valid 'Network' object."
-        assert self.network is None or (self.network.m_flow_network is not None and self.network.cp_f is not None), \
-            "The mass flow rate 'm_flow_network' and heat capacity 'cp_f' must be " \
-            "provided at the instanciation of the 'Network' object."
+        if self.boundary_condition == 'MIFT':
+            assert not (self.m_flow_network is None and self.m_flow_borehole is None), \
+                "The mass flow rate 'm_flow_borehole' or 'm_flow_network' must " \
+                "be provided when using the 'MIFT' boundary condition."
+            assert not (self.m_flow_network is not None and self.m_flow_borehole is not None), \
+                "Only one of 'm_flow_borehole' or 'm_flow_network' can " \
+                "be provided when using the 'MIFT' boundary condition."
+            assert not self.cp_f is None, \
+                "The heat capacity 'cp_f' must " \
+                "be provided when using the 'MIFT' boundary condition."
+            assert not (type(self.m_flow_borehole) is np.ndarray and not np.size(self.m_flow_borehole, axis=1)==self.network.nInlets), \
+                "The number of mass flow rates in 'm_flow_borehole' must " \
+                "correspond to the number of circuits in the network."
         assert type(self.time) is np.ndarray or isinstance(self.time, (float, np.floating)) or self.time is None, \
             "Time should be a float or an array."
         # self.nSegments can now be an int or list
@@ -1930,6 +2435,27 @@ class _Detailed(_BaseSolver):
         (of type int) as an argument, or an array of size (nSegments[i],) when
         provided with an element of nSegments (of type list).
         Default is :func:`utilities.segment_ratios`.
+    m_flow_borehole : (nInlets,) array or (nMassFlow, nInlets,) array, optional
+        Fluid mass flow rate into each circuit of the network. If a
+        (nMassFlow, nInlets,) array is supplied, the
+        (nMassFlow, nMassFlow,) variable mass flow rate g-functions
+        will be evaluated using the method of Cimmino (2024)
+        [#gFunction-CimBer2024]_. Only required for the 'MIFT' boundary
+         condition. Only one of 'm_flow_borehole' and 'm_flow_network' can be
+         provided.
+        Default is None.
+    m_flow_network : float or (nMassFlow,) array, optional
+        Fluid mass flow rate into the network of boreholes. If an array
+        is supplied, the (nMassFlow, nMassFlow,) variable mass flow
+        rate g-functions will be evaluated using the method of Cimmino
+        (2024) [#gFunction-CimBer2024]_. Only required for the 'MIFT' boundary
+         condition. Only one of 'm_flow_borehole' and 'm_flow_network' can be
+         provided.
+        Default is None.
+    cp_f : float, optional
+        Fluid specific isobaric heat capacity (in J/kg.degC). Only required
+        for the 'MIFT' boundary condition.
+        Default is None.
     approximate_FLS : bool, optional
         Set to true to use the approximation of the FLS solution of Cimmino
         (2021) [#Detailed-Cimmin2021]_. This approximation does not require the
@@ -2193,6 +2719,27 @@ class _Similarities(_BaseSolver):
         (of type int) as an argument, or an array of size (nSegments[i],) when
         provided with an element of nSegments (of type list).
         Default is :func:`utilities.segment_ratios`.
+    m_flow_borehole : (nInlets,) array or (nMassFlow, nInlets,) array, optional
+        Fluid mass flow rate into each circuit of the network. If a
+        (nMassFlow, nInlets,) array is supplied, the
+        (nMassFlow, nMassFlow,) variable mass flow rate g-functions
+        will be evaluated using the method of Cimmino (2024)
+        [#gFunction-CimBer2024]_. Only required for the 'MIFT' boundary
+         condition. Only one of 'm_flow_borehole' and 'm_flow_network' can be
+         provided.
+        Default is None.
+    m_flow_network : float or (nMassFlow,) array, optional
+        Fluid mass flow rate into the network of boreholes. If an array
+        is supplied, the (nMassFlow, nMassFlow,) variable mass flow
+        rate g-functions will be evaluated using the method of Cimmino
+        (2024) [#gFunction-CimBer2024]_. Only required for the 'MIFT' boundary
+         condition. Only one of 'm_flow_borehole' and 'm_flow_network' can be
+         provided.
+        Default is None.
+    cp_f : float, optional
+        Fluid specific isobaric heat capacity (in J/kg.degC). Only required
+        for the 'MIFT' boundary condition.
+        Default is None.
     approximate_FLS : bool, optional
         Set to true to use the approximation of the FLS solution of Cimmino
         (2021) [#Similarities-Cimmin2021]_. This approximation does not require
@@ -3494,6 +4041,27 @@ class _Equivalent(_BaseSolver):
         (of type int) as an argument, or an array of size (nSegments[i],) when
         provided with an element of nSegments (of type list).
         Default is :func:`utilities.segment_ratios`.
+    m_flow_borehole : (nInlets,) array or (nMassFlow, nInlets,) array, optional
+        Fluid mass flow rate into each circuit of the network. If a
+        (nMassFlow, nInlets,) array is supplied, the
+        (nMassFlow, nMassFlow,) variable mass flow rate g-functions
+        will be evaluated using the method of Cimmino (2024)
+        [#gFunction-CimBer2024]_. Only required for the 'MIFT' boundary
+         condition. Only one of 'm_flow_borehole' and 'm_flow_network' can be
+         provided.
+        Default is None.
+    m_flow_network : float or (nMassFlow,) array, optional
+        Fluid mass flow rate into the network of boreholes. If an array
+        is supplied, the (nMassFlow, nMassFlow,) variable mass flow
+        rate g-functions will be evaluated using the method of Cimmino
+        (2024) [#gFunction-CimBer2024]_. Only required for the 'MIFT' boundary
+         condition. Only one of 'm_flow_borehole' and 'm_flow_network' can be
+         provided.
+        Default is None.
+    cp_f : float, optional
+        Fluid specific isobaric heat capacity (in J/kg.degC). Only required
+        for the 'MIFT' boundary condition.
+        Default is None.
     approximate_FLS : bool, optional
         Set to true to use the approximation of the FLS solution of Cimmino
         (2021) [#Equivalent-Cimmin2021]_. This approximation does not require
@@ -3808,8 +4376,6 @@ def find_groups(self, tol=1e-6):
             self.network = _EquivalentNetwork(
                 self.boreholes,
                 pipes,
-                m_flow_network=self.network.m_flow_network,
-                cp_f=self.network.cp_f,
                 nSegments=self.nBoreSegments[0],
                 segment_ratios=self.segment_ratios[0])
 
@@ -4228,9 +4794,9 @@ def _check_solver_specific_inputs(self):
             assert np.all(np.array(self.network.c, dtype=int) == -1), \
                 "Solver 'equivalent' is only valid for parallel-connected " \
                 "boreholes."
-            assert type(self.network.m_flow_network) is float \
-                or (type(self.network.m_flow_network) is np.ndarray
-                    and np.allclose(self.network.m_flow_network, self.network.m_flow_network[0])), \
+            assert (self.m_flow_borehole is None
+                    or (self.m_flow_borehole.ndim==1 and np.allclose(self.m_flow_borehole, self.m_flow_borehole[0]))
+                    or (self.m_flow_borehole.ndim==2 and np.all([np.allclose(self.m_flow_borehole[:, i], self.m_flow_borehole[0, i]) for i in range(self.nBoreholes)]))), \
                 "Mass flow rates into the network must be equal for all " \
                 "boreholes."
             # Use the total network mass flow rate.
diff --git a/pygfunction/networks.py b/pygfunction/networks.py
index af37aec3..9fa5884c 100644
--- a/pygfunction/networks.py
+++ b/pygfunction/networks.py
@@ -1391,8 +1391,8 @@ def coefficients_network_heat_extraction_rate(
             # 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}]
-            a_in = np.reshape(np.sum(b_in*self.wBoreholes, axis=0), (1,-1))
-            a_b = np.reshape(np.sum(b_b*self.wBoreholes, axis=0), (1,-1))
+            a_in = np.reshape(self.wBoreholes.T @ b_in, (1,-1))
+            a_b = np.reshape(self.wBoreholes.T @ b_b, (1,-1))
 
             # Store coefficients
             self._set_stored_coefficients(


Test patch:
diff --git a/tests/conftest.py b/tests/conftest.py
index b14f4d4e..d531e3ca 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -524,8 +524,7 @@ def single_borehole_network(single_borehole):
         UTubes.append(UTube)
     # Initialize network
     network = gt.networks.Network(
-        boreField, UTubes, bore_connectivity=bore_connectivity,
-        m_flow_network=m_flow_network, cp_f=fluid.cp)
+        boreField, UTubes, bore_connectivity=bore_connectivity)
     return network
 
 
@@ -566,8 +565,7 @@ def single_borehole_network_short(single_borehole_short):
         UTubes.append(UTube)
     # Initialize network
     network = gt.networks.Network(
-        boreField, UTubes, bore_connectivity=bore_connectivity,
-        m_flow_network=m_flow_network, cp_f=fluid.cp)
+        boreField, UTubes, bore_connectivity=bore_connectivity)
     return network
 
 
@@ -608,8 +606,7 @@ def ten_boreholes_network_rectangular(ten_boreholes_rectangular):
         UTubes.append(UTube)
     # Initialize network
     network = gt.networks.Network(
-        boreField, UTubes, bore_connectivity=bore_connectivity,
-        m_flow_network=m_flow_network, cp_f=fluid.cp)
+        boreField, UTubes, bore_connectivity=bore_connectivity)
     return network
 
 
@@ -650,8 +647,7 @@ def ten_boreholes_network_rectangular_series(ten_boreholes_rectangular):
         UTubes.append(UTube)
     # Initialize network
     network = gt.networks.Network(
-        boreField, UTubes, bore_connectivity=bore_connectivity,
-        m_flow_network=m_flow_network, cp_f=fluid.cp)
+        boreField, UTubes, bore_connectivity=bore_connectivity)
     return network
 
 
@@ -692,6 +688,5 @@ def three_boreholes_network_series_unequal(three_boreholes_unequal):
         UTubes.append(UTube)
     # Initialize network
     network = gt.networks.Network(
-        boreField, UTubes, bore_connectivity=bore_connectivity,
-        m_flow_network=m_flow_network, cp_f=fluid.cp)
+        boreField, UTubes, bore_connectivity=bore_connectivity)
     return network
diff --git a/tests/gfunction_test.py b/tests/gfunction_test.py
index 151f1d6c..e60df9d9 100644
--- a/tests/gfunction_test.py
+++ b/tests/gfunction_test.py
@@ -155,69 +155,113 @@ def test_gfunctions_UHTR(field, method, opts, expected, request):
 # Test 'MIFT' g-functions for different bore fields using all solvers,
 # unequal/uniform segments, and with/without the FLS approximation
 # The 'equivalent' solver is not applied to series-connected boreholes.
-@pytest.mark.parametrize("field, method, opts, expected", [
+@pytest.mark.parametrize("field, method, m_flow_network, opts, expected", [
         #  'equivalent' solver - unequal segments
-        ('single_borehole_network', 'equivalent', 'unequal_segments', np.array([5.76597302, 6.51058473, 6.73746895])),
-        ('single_borehole_network_short', 'equivalent', 'unequal_segments', np.array([4.17105954, 5.00930075, 5.30832133])),
-        ('ten_boreholes_network_rectangular', 'equivalent', 'unequal_segments', np.array([12.66229998, 18.57852681, 20.33535907])),
+        ('single_borehole_network', 'equivalent', 0.05, 'unequal_segments', np.array([5.76597302, 6.51058473, 6.73746895])),
+        ('single_borehole_network_short', 'equivalent', 0.05, 'unequal_segments', np.array([4.17105954, 5.00930075, 5.30832133])),
+        ('ten_boreholes_network_rectangular', 'equivalent', 0.25, 'unequal_segments', np.array([12.66229998, 18.57852681, 20.33535907])),
         # 'equivalent' solver - uniform segments
-        ('single_borehole_network', 'equivalent', 'uniform_segments', np.array([5.78644676, 6.5311583, 6.75699875])),
-        ('single_borehole_network_short', 'equivalent', 'uniform_segments', np.array([4.17553236, 5.01476781, 5.31381287])),
-        ('ten_boreholes_network_rectangular', 'equivalent', 'uniform_segments', np.array([12.931553, 18.8892892, 20.63810364])),
+        ('single_borehole_network', 'equivalent', 0.05, 'uniform_segments', np.array([5.78644676, 6.5311583, 6.75699875])),
+        ('single_borehole_network_short', 'equivalent', 0.05, 'uniform_segments', np.array([4.17553236, 5.01476781, 5.31381287])),
+        ('ten_boreholes_network_rectangular', 'equivalent', 0.25, 'uniform_segments', np.array([12.931553, 18.8892892, 20.63810364])),
         #  'equivalent' solver - unequal segments, FLS approximation
-        ('single_borehole_network', 'equivalent', 'unequal_segments_approx', np.array([5.76596769, 6.51061169, 6.73731276])),
-        ('single_borehole_network_short', 'equivalent', 'unequal_segments_approx', np.array([4.17105984, 5.00931374, 5.30816983])),
-        ('ten_boreholes_network_rectangular', 'equivalent', 'unequal_segments_approx', np.array([12.66228879, 18.57863253, 20.33526092])),
+        ('single_borehole_network', 'equivalent', 0.05, 'unequal_segments_approx', np.array([5.76596769, 6.51061169, 6.73731276])),
+        ('single_borehole_network_short', 'equivalent', 0.05, 'unequal_segments_approx', np.array([4.17105984, 5.00931374, 5.30816983])),
+        ('ten_boreholes_network_rectangular', 'equivalent', 0.25, 'unequal_segments_approx', np.array([12.66228879, 18.57863253, 20.33526092])),
         #  'equivalent' solver - uniform segments, FLS approximation
-        ('single_borehole_network', 'equivalent', 'uniform_segments_approx', np.array([5.78644007, 6.53117706, 6.75684456])),
-        ('single_borehole_network_short', 'equivalent', 'uniform_segments_approx', np.array([4.17553118, 5.01478084, 5.31366109])),
-        ('ten_boreholes_network_rectangular', 'equivalent', 'uniform_segments_approx', np.array([12.93153466, 18.88937176, 20.63801163])),
+        ('single_borehole_network', 'equivalent', 0.05, 'uniform_segments_approx', np.array([5.78644007, 6.53117706, 6.75684456])),
+        ('single_borehole_network_short', 'equivalent', 0.05, 'uniform_segments_approx', np.array([4.17553118, 5.01478084, 5.31366109])),
+        ('ten_boreholes_network_rectangular', 'equivalent', 0.25, 'uniform_segments_approx', np.array([12.93153466, 18.88937176, 20.63801163])),
         # 'similarities' solver - unequal segments
-        ('single_borehole_network', 'similarities', 'unequal_segments', np.array([5.76597302, 6.51058473, 6.73746895])),
-        ('single_borehole_network_short', 'similarities', 'unequal_segments', np.array([4.17105954, 5.00930075, 5.30832133])),
-        ('ten_boreholes_network_rectangular', 'similarities', 'unequal_segments', np.array([12.66229998, 18.57852681, 20.33535907])),
-        ('ten_boreholes_network_rectangular_series', 'similarities', 'unequal_segments', np.array([3.19050169, 8.8595362, 10.84379419])),
+        ('single_borehole_network', 'similarities', 0.05, 'unequal_segments', np.array([5.76597302, 6.51058473, 6.73746895])),
+        ('single_borehole_network_short', 'similarities', 0.05, 'unequal_segments', np.array([4.17105954, 5.00930075, 5.30832133])),
+        ('ten_boreholes_network_rectangular', 'similarities', 0.25, 'unequal_segments', np.array([12.66229998, 18.57852681, 20.33535907])),
+        ('ten_boreholes_network_rectangular_series', 'similarities', 0.05, 'unequal_segments', np.array([3.19050169, 8.8595362, 10.84379419])),
         # 'similarities' solver - uniform segments
-        ('single_borehole_network', 'similarities', 'uniform_segments', np.array([5.78644676, 6.5311583, 6.75699875])),
-        ('single_borehole_network_short', 'similarities', 'uniform_segments', np.array([4.17553236, 5.01476781, 5.31381287])),
-        ('ten_boreholes_network_rectangular', 'similarities', 'uniform_segments', np.array([12.931553, 18.8892892, 20.63810364])),
-        ('ten_boreholes_network_rectangular_series', 'similarities', 'uniform_segments', np.array([3.22186337, 8.92628494, 10.91644607])),
+        ('single_borehole_network', 'similarities', 0.05, 'uniform_segments', np.array([5.78644676, 6.5311583, 6.75699875])),
+        ('single_borehole_network_short', 'similarities', 0.05, 'uniform_segments', np.array([4.17553236, 5.01476781, 5.31381287])),
+        ('ten_boreholes_network_rectangular', 'similarities', 0.25, 'uniform_segments', np.array([12.931553, 18.8892892, 20.63810364])),
+        ('ten_boreholes_network_rectangular_series', 'similarities', 0.05, 'uniform_segments', np.array([3.22186337, 8.92628494, 10.91644607])),
         # 'similarities' solver - unequal segments, FLS approximation
-        ('single_borehole_network', 'similarities', 'unequal_segments_approx', np.array([5.76596769, 6.51061169, 6.73731276])),
-        ('single_borehole_network_short', 'similarities', 'unequal_segments_approx', np.array([4.17105984, 5.00931374, 5.30816983])),
-        ('ten_boreholes_network_rectangular', 'similarities', 'unequal_segments_approx', np.array([12.66136057, 18.57792276, 20.33429034])),
-        ('ten_boreholes_network_rectangular_series', 'similarities', 'unequal_segments_approx', np.array([3.19002567, 8.85783554, 10.84222402])),
+        ('single_borehole_network', 'similarities', 0.05, 'unequal_segments_approx', np.array([5.76596769, 6.51061169, 6.73731276])),
+        ('single_borehole_network_short', 'similarities', 0.05, 'unequal_segments_approx', np.array([4.17105984, 5.00931374, 5.30816983])),
+        ('ten_boreholes_network_rectangular', 'similarities', 0.25, 'unequal_segments_approx', np.array([12.66136057, 18.57792276, 20.33429034])),
+        ('ten_boreholes_network_rectangular_series', 'similarities', 0.05, 'unequal_segments_approx', np.array([3.19002567, 8.85783554, 10.84222402])),
         # 'similarities' solver - uniform segments, FLS approximation
-        ('single_borehole_network', 'similarities', 'uniform_segments_approx', np.array([5.78644007, 6.53117706, 6.75684456])),
-        ('single_borehole_network_short', 'similarities', 'uniform_segments_approx', np.array([4.17553118, 5.01478084, 5.31366109])),
-        ('ten_boreholes_network_rectangular', 'similarities', 'uniform_segments_approx', np.array([12.93064329, 18.88844718, 20.63710493])),
-        ('ten_boreholes_network_rectangular_series', 'similarities', 'uniform_segments_approx', np.array([3.22139028, 8.92448715, 10.91485947])),
+        ('single_borehole_network', 'similarities', 0.05, 'uniform_segments_approx', np.array([5.78644007, 6.53117706, 6.75684456])),
+        ('single_borehole_network_short', 'similarities', 0.05, 'uniform_segments_approx', np.array([4.17553118, 5.01478084, 5.31366109])),
+        ('ten_boreholes_network_rectangular', 'similarities', 0.25, 'uniform_segments_approx', np.array([12.93064329, 18.88844718, 20.63710493])),
+        ('ten_boreholes_network_rectangular_series', 'similarities', 0.05, 'uniform_segments_approx', np.array([3.22139028, 8.92448715, 10.91485947])),
         # 'detailed' solver - unequal segments
-        ('single_borehole_network', 'detailed', 'unequal_segments', np.array([5.76597302, 6.51058473, 6.73746895])),
-        ('single_borehole_network_short', 'detailed', 'unequal_segments', np.array([4.17105954, 5.00930075, 5.30832133])),
-        ('ten_boreholes_network_rectangular', 'detailed', 'unequal_segments', np.array([12.66229998, 18.57852681, 20.33535907])),
-        ('ten_boreholes_network_rectangular_series', 'detailed', 'unequal_segments', np.array([3.19050169, 8.8595362, 10.84379419])),
+        ('single_borehole_network', 'detailed', 0.05, 'unequal_segments', np.array([5.76597302, 6.51058473, 6.73746895])),
+        ('single_borehole_network_short', 'detailed', 0.05, 'unequal_segments', np.array([4.17105954, 5.00930075, 5.30832133])),
+        ('ten_boreholes_network_rectangular', 'detailed', 0.25, 'unequal_segments', np.array([12.66229998, 18.57852681, 20.33535907])),
+        ('ten_boreholes_network_rectangular_series', 'detailed', 0.05, 'unequal_segments', np.array([3.19050169, 8.8595362, 10.84379419])),
         # 'detailed' solver - uniform segments
-        ('single_borehole_network', 'detailed', 'uniform_segments', np.array([5.78644676, 6.5311583, 6.75699875])),
-        ('single_borehole_network_short', 'detailed', 'uniform_segments', np.array([4.17553236, 5.01476781, 5.31381287])),
-        ('ten_boreholes_network_rectangular', 'detailed', 'uniform_segments', np.array([12.931553, 18.8892892, 20.63810364])),
-        ('ten_boreholes_network_rectangular_series', 'detailed', 'uniform_segments', np.array([3.22186337, 8.92628494, 10.91644607])),
+        ('single_borehole_network', 'detailed', 0.05, 'uniform_segments', np.array([5.78644676, 6.5311583, 6.75699875])),
+        ('single_borehole_network_short', 'detailed', 0.05, 'uniform_segments', np.array([4.17553236, 5.01476781, 5.31381287])),
+        ('ten_boreholes_network_rectangular', 'detailed', 0.25, 'uniform_segments', np.array([12.931553, 18.8892892, 20.63810364])),
+        ('ten_boreholes_network_rectangular_series', 'detailed', 0.05, 'uniform_segments', np.array([3.22186337, 8.92628494, 10.91644607])),
         # 'detailed' solver - unequal segments, FLS approximation
-        ('single_borehole_network', 'detailed', 'unequal_segments_approx', np.array([5.76596769, 6.51061169, 6.73731276])),
-        ('single_borehole_network_short', 'detailed', 'unequal_segments_approx', np.array([4.17105984, 5.00931374, 5.30816983])),
-        ('ten_boreholes_network_rectangular', 'detailed', 'unequal_segments_approx', np.array([12.66136057, 18.57792276, 20.33429034])),
-        ('ten_boreholes_network_rectangular_series', 'detailed', 'unequal_segments_approx', np.array([3.19002567, 8.85783554, 10.84222402])),
+        ('single_borehole_network', 'detailed', 0.05, 'unequal_segments_approx', np.array([5.76596769, 6.51061169, 6.73731276])),
+        ('single_borehole_network_short', 'detailed', 0.05, 'unequal_segments_approx', np.array([4.17105984, 5.00931374, 5.30816983])),
+        ('ten_boreholes_network_rectangular', 'detailed', 0.25, 'unequal_segments_approx', np.array([12.66136057, 18.57792276, 20.33429034])),
+        ('ten_boreholes_network_rectangular_series', 'detailed', 0.05, 'unequal_segments_approx', np.array([3.19002567, 8.85783554, 10.84222402])),
         # 'detailed' solver - uniform segments, FLS approximation
-        ('single_borehole_network', 'detailed', 'uniform_segments_approx', np.array([5.78644007, 6.53117706, 6.75684456])),
-        ('single_borehole_network_short', 'detailed', 'uniform_segments_approx', np.array([4.17553118, 5.01478084, 5.31366109])),
-        ('ten_boreholes_network_rectangular', 'detailed', 'uniform_segments_approx', np.array([12.93064329, 18.88844718, 20.63710493])),
-        ('ten_boreholes_network_rectangular_series', 'detailed', 'uniform_segments_approx', np.array([3.22139028, 8.92448715, 10.91485947])),
+        ('single_borehole_network', 'detailed', 0.05, 'uniform_segments_approx', np.array([5.78644007, 6.53117706, 6.75684456])),
+        ('single_borehole_network_short', 'detailed', 0.05, 'uniform_segments_approx', np.array([4.17553118, 5.01478084, 5.31366109])),
+        ('ten_boreholes_network_rectangular', 'detailed', 0.25, 'uniform_segments_approx', np.array([12.93064329, 18.88844718, 20.63710493])),
+        ('ten_boreholes_network_rectangular_series', 'detailed', 0.05, 'uniform_segments_approx', np.array([3.22139028, 8.92448715, 10.91485947])),
+    ])
+def test_gfunctions_MIFT(
+        field, method, m_flow_network, opts, expected, request):
+    # Extract the bore field from the fixture
+    network = request.getfixturevalue(field)
+    # Extract the g-function options from the fixture
+    options = request.getfixturevalue(opts)
+    # Fluid is propylene-glycol (20 %) at 20 degC
+    fluid = gt.media.Fluid('MPG', 20.)
+    # Mean borehole length [m]
+    H_mean = np.mean([b.H for b in network.b])
+    alpha = 1e-6    # Ground thermal diffusivity [m2/s]
+    # Bore field characteristic time [s]
+    ts = H_mean**2 / (9 * alpha)
+    # Times for the g-function [s]
+    time = np.array([0.1, 1., 10.]) * ts
+    # g-Function
+    gFunc = gt.gfunction.gFunction(
+        network, alpha, time=time, m_flow_network=m_flow_network,
+        cp_f=fluid.cp, method=method, options=options, boundary_condition='MIFT')
+    assert np.allclose(gFunc.gFunc, expected)
+
+
+# Test 'MIFT' g-functions for different bore fields using all solvers, unequal
+# segments, and with the FLS approximation for variable mass flow rate
+# g-functions
+@pytest.mark.parametrize("field, method, m_flow_network, opts, expected", [
+        # 'detailed' solver - unequal segments, FLS approximation
+        ('single_borehole_network', 'detailed', np.array([-0.25, 0.05]), 'unequal_segments_approx', np.array([[[5.60566489, 6.38071943, 6.62706326], [5.603355, 6.3615959, 6.5979712]], [[5.60284597, 6.36233529, 6.5996188], [5.76596769, 6.51061169, 6.73731276]]])),
+        ('single_borehole_network_short', 'detailed', np.array([-0.25, 0.05]), 'unequal_segments_approx', np.array([[[4.17130904, 5.0108666, 5.31097968], [4.1708935, 5.00979555, 5.30928557]], [[4.1708873, 5.0097772, 5.30925837], [4.17105984, 5.00931374, 5.30816983]]])),
+        ('ten_boreholes_network_rectangular', 'detailed', np.array([-0.5, 0.25]), 'unequal_segments_approx', np.array([[[11.19891898, 17.50417429, 19.51221751], [11.61378246, 17.69400368, 19.56564714]], [[11.44417287, 17.5662909, 19.47928482], [12.66136057, 18.57792276, 20.33429034]]])),
+        ('ten_boreholes_network_rectangular_series', 'detailed', np.array([-0.25, 0.05]), 'unequal_segments_approx', np.array([[[7.57467211, 13.97746499, 16.08569783], [25.31574066, 31.796452, 33.86566759]], [[17.37750733, 23.69326533, 25.75080157], [3.19002567, 8.85783554, 10.84222402]]])),
+        # 'similarities' solver - unequal segments, FLS approximation
+        ('single_borehole_network', 'similarities', np.array([-0.25, 0.05]), 'unequal_segments_approx', np.array([[[5.60566489, 6.38071943, 6.62706326], [5.603355, 6.3615959, 6.5979712]], [[5.60284597, 6.36233529, 6.5996188], [5.76596769, 6.51061169, 6.73731276]]])),
+        ('single_borehole_network_short', 'similarities', np.array([-0.25, 0.05]), 'unequal_segments_approx', np.array([[[4.17130904, 5.0108666, 5.31097968], [4.1708935, 5.00979555, 5.30928557]], [[4.1708873, 5.0097772, 5.30925837], [4.17105984, 5.00931374, 5.30816983]]])),
+        ('ten_boreholes_network_rectangular', 'similarities', np.array([-0.5, 0.25]), 'unequal_segments_approx', np.array([[[11.19891898, 17.50417429, 19.51221751], [11.61378246, 17.69400368, 19.56564714]], [[11.44417287, 17.5662909, 19.47928482], [12.66136057, 18.57792276, 20.33429034]]])),
+        ('ten_boreholes_network_rectangular_series', 'similarities', np.array([-0.25, 0.05]), 'unequal_segments_approx', np.array([[[7.57467211, 13.97746499, 16.08569783], [25.31574066, 31.796452, 33.86566759]], [[17.37750733, 23.69326533, 25.75080157], [3.19002567, 8.85783554, 10.84222402]]])),
+        # 'equivalent' solver - unequal segments, FLS approximation
+        ('single_borehole_network', 'equivalent', np.array([-0.25, 0.05]), 'unequal_segments_approx', np.array([[[5.60566489, 6.38071943, 6.62706326], [5.603355, 6.3615959, 6.5979712]], [[5.60284597, 6.36233529, 6.5996188], [5.76596769, 6.51061169, 6.73731276]]])),
+        ('single_borehole_network_short', 'equivalent', np.array([-0.25, 0.05]), 'unequal_segments_approx', np.array([[[4.17130904, 5.0108666, 5.31097968], [4.1708935, 5.00979555, 5.30928557]], [[4.1708873, 5.0097772, 5.30925837], [4.17105984, 5.00931374, 5.30816983]]])),
+        ('ten_boreholes_network_rectangular', 'equivalent', np.array([-0.5, 0.25]), 'unequal_segments_approx', np.array([[[11.19974879, 17.50502384, 19.51351783], [11.61465222, 17.69478305, 19.56679743]], [[11.44502591, 17.56709008, 19.48049408], [12.66228879, 18.57863253, 20.33526092]]])),
     ])
-def test_gfunctions_MIFT(field, method, opts, expected, request):
+def test_gfunctions_MIFT_variable_mass_flow_rate(
+        field, method, m_flow_network, opts, expected, request):
     # Extract the bore field from the fixture
     network = request.getfixturevalue(field)
     # Extract the g-function options from the fixture
     options = request.getfixturevalue(opts)
+    # Fluid is propylene-glycol (20 %) at 20 degC
+    fluid = gt.media.Fluid('MPG', 20.)
     # Mean borehole length [m]
     H_mean = np.mean([b.H for b in network.b])
     alpha = 1e-6    # Ground thermal diffusivity [m2/s]
@@ -227,8 +271,8 @@ def test_gfunctions_MIFT(field, method, opts, expected, request):
     time = np.array([0.1, 1., 10.]) * ts
     # g-Function
     gFunc = gt.gfunction.gFunction(
-        network, alpha, time=time, method=method, options=options,
-        boundary_condition='MIFT')
+        network, alpha, time=time, m_flow_network=m_flow_network,
+        cp_f=fluid.cp, method=method, options=options, boundary_condition='MIFT')
     assert np.allclose(gFunc.gFunc, expected)