67. The Barends Problem
Underground thermal energy storage (UTES) systems, sometimes referred to as aquifer thermal energy storage (ATES) systems, can help bolster sustainable energy supply portfolios. In UTES, energy can be injected into an aquifer at one time to be withdrawn as needed at a later time. However, “thermal bleeding,” or in other words conductive heat exchange with over and underlying geological formations may degrade the efficacy of UTES systems as injected heat diffuses away from the storage site. To better understand the potential of UTES at various sites, numerical models are often used to verify field tests and further explore the suitability of a site for UTES (Barends, F. B. J., 2010). This demonstration of MODFLOW 6 is patterned after the “Barends” problem.
67.1. Example description
This problem uses a 2D vertical cross-section comprised of two zones - a 100 \(m\) thick saturated zone that is overlain by a 100 \(m\) thick overburden (Figure 67.1). Within the saturated zone, groundwater flow is simulated as a 1-dimensional (1D) flow-field moving from left to right. The bottom boundary of the groundwater reservoir is considered sealed for both flow and heat. However, the upper boundary of the groundwater reservoir - the interface with the overburden - is sealed for flow but not for heat transfer.

Figure 67.1 A simplified depiction of the DIS grid (Langevin et al., 2017) showing basic model components that corresponds to the (Barends, F. B. J., 2010) problem. For this setup, 100 layers are used to represent the groundwater reservoir; however, thickness-weighted injection/extraction amounts within each layer of the groundwater reservoir result in 1D groundwater flow below the overburden. There is no groundwater flow in the overburden, only conductive heat exchange between these two geologic units
For simulating only conductive heat exchange at the
overburden-groundwater reservoir interface, the vertical hydraulic
conductivities are set extremely low throughout the entire model domain
(tab-ex-gwe-barends-01
). The
horizontal hydraulic conductivities are set such that the horizontal
groundwater velocity is \(1.2649 \times 10^{-8} \frac {m}{s}\)
within the groundwater reservoir. Flow is injected in the left-most
groundwater reservoir cells and extracted from the right-most cells at
the same rate. Within the overburden, there is no convection (as already
mentioned) and thermal conduction is 2-dimensional (2D), a significant
departure from the Barends problem setup
(Figure 67.1).
The model domain is initialized with a uniform temperature of \(80^{\circ}C\) at the beginning of the simulation (\(T_0\)). \(30^{\circ}C\) water enters on the left side of the groundwater reservoir at the start of the simulation (\(T_0 > 0\)). As cold water flows into the groundwater reservoir, a thermal gradient between the overburden and saturated zone is established that results in the aforementioned thermal bleeding (conduction) of heat from the overburden into the groundwater reservoir. Geometric grid refinement is added in the vertical direction on both sides of the overburden-groundwater reservoir interface (Figure 67.2)

Figure 67.2 A zoomed-in view of the numerical grid showing enhanced vertical refinement above and below the overburden-groundwater reservoir interface (y=100 \(m\)).
The rate of thermal bleeding at each overburden-groundwater reservoir cell interface changes as the pulse of cold water injected into the groundwater reservoir migrates inward.
67.2. Example Results
After 200 years of simulation time, lower temperatures show an advance into the groundwater reservoir (Figure 67.3). As cooler water is injected into the aquifer, conduction from the overburden into the groundwater reservoir retards the advance of colder temperatures in the upper portion of the groundwater reservoir compared to the lower portion as shown by the curved temperature contours within the groundwater reservoir.

Figure 67.3 Temperatures predicted by the MODFLOW 6 GWE model after 200 years of simulation time.
67.3. References Cited
Barends, F. B. J. (2010). Complete solution for transient heat transport in porous media, following lauwerier’s concept. In Society of petroleum engineers annual technical conference and exhibition, september 19–22, 2010, florence, italy. Deltares; Technical University Delft; Society of Petroleum Engineers. https://doi.org/doi.org/10.2118/134670-MS
Langevin, C. D., Hughes, J. D., Provost, A. M., Banta, E. R., Niswonger, R. G., & Panday, S. (2017). Documentation for the MODFLOW 6 Groundwater Flow (GWF) Model. https://doi.org/10.3133/tm6A55
67.4. Jupyter Notebook
The Jupyter notebook used to create the MODFLOW 6 input files for this example and post-process the results is: