Optimization approaches for the design and operation of openloop shallow geothermal systems
Fabian Böttcher
Kai Zosseder
Thomas Hamacher
The optimization of openloop shallow geothermal systems, which includes both design and operational aspects, is an important research area aimed at improving their efficiency and sustainability and the effective management of groundwater as a shallow geothermal resource. This paper investigates various approaches to address optimization problems arising from these research and implementation questions about GWHP systems. The identified optimization approaches are thoroughly analyzed based on criteria such as computational cost and applicability. Moreover, a novel classification scheme is introduced that categorizes the approaches according to the types of groundwater simulation model and the optimization algorithm used. Simulation models are divided into two types: numerical and simplified (analytical or datadriven) models, while optimization algorithms are divided into gradientbased and derivativefree algorithms. Finally, a comprehensive review of existing approaches in the literature is provided, highlighting their strengths and limitations and offering recommendations for both the use of existing approaches and the development of new, improved ones in this field.
 Article
(666 KB)  Fulltext XML
 BibTeX
 EndNote
Openloop shallow geothermal systems, also known as groundwater heat pumps (GWHPs), have emerged as a promising solution for decarbonizing the residential heating and cooling sector (Russo et al., 2012). The performance of GWHPs is primarily influenced by groundwater temperature (Kim and Nam, 2016), which remains relatively stable throughout the year and is elevated in urban areas due to the subsurface urban heat island effect (Menberg et al., 2013; Epting and Huggenberger, 2013; Böttcher and Zosseder, 2021). These systems harness the thermal energy of the aquifer by extracting groundwater from one or more extraction wells and returning it to the same aquifer via injection wells after heat exchange in a heat pump (Florides and Kalogirou, 2007; Stauffer et al., 2014; García Gil et al., 2022). Since the temperature of the reinjected water is different from that of the extracted (lower in the heating and higher in the cooling case), this results in thermal plumes in the aquifer that propagate downstream along the groundwater flow direction. If these plumes reach the extraction wells of neighboring downstream GWHPs, this can result in either negative or positive thermal interference (Perego et al., 2022). Figure 1 provides an overview of potential thermal interferences that can occur between neighboring systems, depicting scenarios where the operation of downstream systems can be either degraded (negative interference) or enhanced (positive interference).
It is also important to recognize that GWHPs have a thermal impact on groundwater, which serves as a vital source of drinking water in many places (e.g. Blum et al., 2021). To mitigate the aforementioned negative interactions and improve the efficiency and sustainability of thermal groundwater use, resource management strategies need to be implemented (Epting et al., 2020). This includes optimizing the design, particularly well placement, and operation of GWHP systems, since the propagation of thermal plumes is affected by injection well locations, system operation (pumping rates), and aquifer characteristics. For example, optimal well placement can minimize negative thermal interference between neighboring GWHPs (Halilovic et al., 2022a) or even maximize their positive interference (GarcíaGil et al., 2020). Hence, optimization of multiple neighboring systems plays an important role in urban planning strategies aimed at enhancing sustainability. In addition, optimization of GWHP systems is crucial for managing groundwater resources and maintaining the current state of groundwater, i.e. preventing adverse changes in its physical, chemical, and biological properties. Furthermore, it is essential to ensure adequate spacing between wells within the same GWHP to prevent hydraulic and thermal breakthroughs (Böttcher et al., 2019). Thus, optimization of individual systems is also important to maximize their efficiency and sustainability. Optimization of individual GWHP systems and concurrent optimization of multiple neighboring systems are challenging due to the complexity of the resulting optimization problems and the necessity for novel and efficient optimization approaches to solve them.
This paper presents a comprehensive overview of optimization approaches for the design and operation of GWHP systems. The approaches are critically evaluated and compared based on several criteria, and a novel classification scheme is introduced to effectively categorize these approaches. Furthermore, the current status of the approaches found in the literature is presented and possible future research directions are discussed.
GWHP systems affect the groundwater body both hydraulically and thermally (García Gil et al., 2022), which can also affect its chemical and biological conditions to a minor degree (Blum et al., 2021). The hydraulic head increases around injection wells and decreases around extraction wells, which also changes the hydraulic gradient and groundwater flow patterns. Thermal impacts are present due to the previously described thermal plumes. To analyze these impacts of GWHP systems on groundwater conditions, simulation models are commonly used. For a particular system design and operation, a simulation model can quantify its impacts on groundwater and, based on that, analyze the performance of the system. Therefore, simulation models play a crucial role for the computation of the resulting groundwater temperature field and serve as an integral component within the optimization procedures.
These simulation models generally fall into three categories: numerical, analytical, and datadriven. Numerical models use partial differential equations (PDEs) to describe the underlying physical phenomena, i.e. groundwater flow and heat transport in aquifers. The resulting system of PDEs can be solved with general PDE solving software or computational fluid dynamics (CFD) software, but there are also several software packages that include specialized domains of numerical simulation for shallow geothermal resources, such as: FEFLOW (Diersch, 2014) – based on the finite element method (FEM) or PFLOTRAN (Hammond et al., 2012) – based on the finite volume method (FVM). Numerical models can incorporate various complex subsurface conditions, including spatially heterogeneous groundwater parameters (e.g. hydraulic conductivity) and conditions (e.g. velocity, temperature, hydraulic head), complex boundary conditions, coupled physical processes, multiple subsurface layers, etc., while simultaneously simulating thermal and hydraulic effects of GWHP systems on the groundwater body. Therefore, they are closest to reality given sufficient quality of input data, but are generally computationally expensive.
The second category of models uses analytical formulas to approximate numerical solutions and is commonly applied to estimate thermal plumes associated with smaller GWHPs, whose energy consumption is less than 45 000 kWh per year (Ohmer et al., 2022). Due to their analytical nature, these models offer significant computational advantages over numerical models. In Pophillat et al. (2020) three prominent analytical models for estimating GWHP thermal plumes were analyzed and compared. These models include the radial heat transport model (RHM) (Guimerà et al., 2007), the linear advective heat transport model (LAHM) (Kinzelbach, 1992), and the planar advective heat transport model (PAHM) (Hähnlein et al., 2010). The authors concluded that although analytical solutions are less accurate compared to numerical models, they still have value for evaluating the thermal impact of GWHPs. Analytical solutions are particularly useful for performing initial assessments of potential negative interference between neighboring GWHPs (Pophillat et al., 2020).
Finally, datadriven models are gaining popularity in this area of research, primarily due to the emergence of machine learning. A common example is the use of neural networks (NNs) to predict thermal plumes (Russo et al., 2014; Leiteritz et al., 2022; Davis et al., 2023). Datadriven models, such as NNs, offer the advantage of fast evaluation, but rely on extensive training data and require additional time for the training process. Acquiring this training data is often challenging due to the limited measurement and monitoring of hydrogeological data. One possible solution is the use of physicsinformed neural networks (PINNs) that integrate physical laws driven by PDEs, mitigating the need for extensive training data (Raissi et al., 2017).
This section provides a comprehensive analysis of two key aspects related to the optimization of GWHP systems. First, in Sect. 3.1, the underlying optimization problems are discussed. Second, in Sect. 3.2, a detailed overview of the approaches for solving these optimization problems is provided. In the following section, we present a generalized problem related to the optimization of GWHP systems, which prepares the way for further analysis in subsequent sections.
3.1 Optimization problems
The highlevel optimization problem concerning GWHP systems can be formulated as follows:
where x_{d} is the vector of optimization variables related to the design of GWHP system(s), x_{o} is the vector of optimization variables related to the operation of GWHP system(s), f_{obj} is the objective function to be minimized, f_{cost} is the function describing technical costs, f_{env} is the function describing negative environmental impacts, α_{1} and α_{2} are the weighting factors, ℱ_{sim} is the simulation model in a residual form, g are the inequality constraints, h are the equality constraints.
In this generalized problem, we differentiate between two types of optimization variables: design variables x_{d} and operational variables x_{o}. An example of design variables are the number and spatial layout of GWHP wells, while an example of operational variables are the pumping rates of each well. The design variables are constant in time, whereas the operational variables are usually timedependent.
The objective function f_{obj} contains two parts: f_{cost}, accounting for the technical costs of GWHP systems, and f_{env}, accounting for the negative environmental impacts. The term f_{cost} can represent various costs associated with the installation and operation of GWHPs, which can be reduced through different means, such as proper sizing of systems (reduced investment costs) or optimal operation of systems (increased efficiency and lifetime). On the other hand, the term f_{env} covers different environmental categories, such as negative impacts on groundwater or CO_{2} emissions indirectly caused by the operation of GWHP systems. It should be noted that environmental considerations are usually incorporated into the problem through constraints and not directly within the objective function.
The simulation model ℱ_{sim} (see Sect. 2) that describes the subsurface phenomena is incorporated into the optimization as a single or multiple equality constraints. This model can be of any type discussed previously: numerical, analytical, or datadriven, and it can also have an explicit form, such as PDEs or algebraic equations, or an implicit ”blackbox” form, such as numerical simulation tools or NNs.
In addition to the simulation model, other inequality g or equality h constraints may be present in the optimization problem. These can be technical constraints, such as upper and lower limits on pumping rates, regulatory constraints, such as the maximum allowed change in groundwater temperatures, or any other additional constraints.
Depending on how certain elements are specified in the generalized problem (Eq. 1), the resulting optimization problems can be classified according to different criteria:

Optimization variables: if the only optimization variables are the design parameters x_{d} and the operation of the system(s) is predefined, i.e. x_{o} is fixed, the problem (Eq. 1) becomes a design optimization problem. On the other hand, the problem becomes an optimal control problem when the system design is specified and the operating parameters x_{o} are the optimization variables. Finally, simultaneous optimization of the design and operation of GWHP systems is possible, i.e. considering both x_{d} and x_{o} as optimization variables. This generally leads to improved optimal solutions since there are more degrees of freedom to be optimized, but the resulting problems are usually more difficult to solve due to increased problem complexity (e.g. from linear to nonlinear).

Objective function: in problem (Eq. 1), the objective function f_{obj} is a weighted sum of technical costs f_{cost} and quantified negative environmental impacts f_{env}. Setting one of the weights α to 0, the problem becomes either a purely economic or an ecological optimization problem, i.e. a singleobjective optimization problem. If both weighting factors are kept positive, then both the cost and the environmental impact are minimized simultaneously and the problem becomes a type of multiobjective optimization. This means that by changing the values of α_{1} and α_{2}, different Paretooptimal solutions are obtained (Marler and Arora, 2010).

Application: the application types can be divided into two main groups: optimization of a single standalone system and optimization of multiple neighboring GWHP systems. The application type directly changes the format of the optimization variables and the objective function. In addition, different applications may involve different optimization constraints, such as the threshold for negative interference between neighboring systems in the case of optimizing multiple systems.

Simulation model: in the mathematical sense, the choice of simulation model ℱ_{sim} fundamentally changes the type of the optimization problem. These optimization problem types belong to different branches of optimization and therefore require corresponding optimization approaches to be solved efficiently. For this reason, the entire following section is dedicated to the optimization approaches for the problem (Eq. 1).
3.2 Optimization approaches
In this study, the term “optimization approach” is considered to encompass not only the specific methodology used to solve a given optimization problem, such as the choice of an algorithm, but also the way in which the problem is formulated, which includes the selection of a groundwater simulation model. The classification of optimization approaches is shown in Fig. 2, where four different classes are identified. The categorization is based on the simulation model used, whether it is a PDE model or a simplified model, and the optimization algorithm employed, either gradientbased or derivativefree algorithms. In the following, each of these four classes is explained and references to relevant literature sources are provided.
Class I comprises optimization approaches where the simulation model is a numerical PDE model, and the optimization is performed using gradientbased algorithms. These approaches are referred to as PDEconstrained optimization (PDECO) problems, which are recognized as the most mathematically complex problems of the four classes considered. The complexity arises due to the multidisciplinary nature of these problems, necessitating expertise in several areas, including computational optimization, functional analysis, and numerical analysis. For example, stateoftheart groundwater simulation tools usually lack the automatic provision of gradient information, requiring users to estimate gradients manually. There are two main methods to solve this problem: automatic differentiation of existing simulation tools (Naumann, 2011) or the development of custom numerical simulators within frameworks such as Firedrake (Rathgeber et al., 2017) and FEniCS (Logg et al., 2012), which can automatically provide the required gradients. Therefore, a comprehensive understanding of PDE solving is essential in the initial stages of developing a Class I approach. Various strategies exist for solving PDECO problems, including fullspace and reducedspace methods, as well as discretizethenoptimize and optimizethendiscretize approaches (Hinze et al., 2008). For more information on PDECO, the reader is referred to the books by Hinze et al. (2008) and Tröltzsch (2010).
Class II includes approaches that utilize simplified models for groundwater simulation and employ gradientbased optimization algorithms to solve the underlying optimization problems. In this context, the simplified models primarily take the form of analytical models (see Sect. 2). These models are expressed through analytical formulas, allowing for direct integration into the optimization problem. The conceptualization of the overall optimization problem determines the resulting problems, which typically correspond to wellestablished types of mathematical programming (optimization) problems, such as linear programming (LP), mixedinteger linear programming (MILP), quadratic programming (QP), and similar types. These problems are extensively studied in the optimization community, and consequently, efficient algorithms and solvers (software implementations of the algorithms) are readily available. Comprehensive information on these types of optimization problems can be found in numerous literature sources, including the books by Nocedal and Wright (1999), Schrijver (1998) and Bonnans et al. (2006).
Class III encompasses approaches that combine PDE simulation models with derivativefree optimization (DFO) algorithms. This form of optimization is commonly referred to as simulationbased optimization, where the simulation model is treated as a blackbox, meaning that only the inputs and outputs of the simulator are observed and used by the optimization algorithm to guide the optimization process. As a result, the term blackbox optimization (BBO) can also be used, which refers to optimization problems where either the objective function or some constraints are treated as blackboxes. However, it is important to note that the blackbox in BBO is not limited to numerical simulation models. For instance, Class IV also falls under the umbrella of BBO. Furthermore, the terms BBO and DFO are closely interconnected and can be used interchangeably in certain contexts. These distinct terminologies have emerged over time, highlighting different aspects: the conceptual characteristics of BBO, and the algorithmic features of DFO. DFO algorithms, as the name suggests, do not rely on derivative information during optimization iterations to determine optimal solutions; instead, they use only the values of the objective function and constraints. These algorithms include: heuristic methods, e.g. genetic algorithms, particle swarm optimization, simulated annealing, etc. (BozorgHaddad et al., 2017); direct search methods, e.g. MADS algorithm (Le Digabel, 2011); and modelbased methods, e.g. modelbased trustregion (Conn et al., 2000). For more detailed information on DFO and BBO, interested readers are referred to Audet and Hare (2017) and Conn et al. (2009).
Class IV comprises optimization approaches that involve the coupling of simplified groundwater simulation models, such as analytical models or NNs, with DFO algorithms. Similar to Class III, Class IV falls into the DFO and BBO branches of optimization. Consequently, the same or similar optimization algorithms can be applied to both classes. The main difference between the two classes lies in the fidelity of the simulation models used. Class III uses highfidelity models, while Class IV relies on lowfidelity models. Thus, the computational cost associated with evaluating potential solution candidates during optimization iterations is significantly lower in Class IV. It is important to note that Class IV has similarities to situations where modelbased DFO algorithms are applied in Class III. The difference is that the approximate models in Class III are constructed dynamically during the optimization iterations, based on evaluations of the PDE model. In contrast, in Class IV, the simplified models are predefined and remain constant throughout the optimization process.
Finally, it should be mentioned that the four introduced classes do not encompass all conceivable approaches, since combined approaches also exist. For example, the solution obtained from Class II can serve as an initialization for the optimization process in Class I. However, within the context of this study, the division into four distinct classes seems both logical and practical, since there are substantial differences between these classes. The following section reviews previous research studies on the optimization of GWHP systems.
3.2.1 Current status of the approaches used for GWHP systems
Despite the increasing importance of GWHP optimization as a research area, the number of existing studies on this topic remains limited. Park et al. (2020) propose a simulationbased optimization approach to optimize pumping rates for a single GWHP system. The approach couples a numerical groundwater simulation model with a genetic optimization algorithm. Furthermore, the same approach was extended in Park et al. (2021) to optimize both well locations and pumping rates within a single system. Since the approach in these two studies uses a PDE simulation model along with a DFO algorithm, it falls under Class III of the proposed classification.
To date, only one research study has been identified that applies the approach of Class I, i.e. the PDECO framework. This study (Halilovic et al., 2022a) introduces a novel approach for concurrently optimizing the well locations of multiple neighboring systems. The approach was illustrated using a case study with ten systems, where the optimization objective is to minimize negative interactions between systems and maximize the overall efficiency of all systems. The proposed approach uses the adjoint method to efficiently compute gradients from the numerical simulation model, which are required by the optimization algorithm.
There is also only one research study that implements the approach belonging to Class II. In Halilovic et al. (2023), the authors introduce an approach that integrates an analytical groundwater simulation model directly into the optimization problem. Specifically, the analytical model used to calculate thermal plumes is the LAHM model (Kinzelbach, 1992) and the resulting optimization problem is formulated as an MILP problem. The study applies this approach to optimize the locations of systems and their associated wells within an urban area comprising 56 potential systems. The objective is to satisfy relevant regulations while maximizing heat extraction from the aquifer. An opensource implementation of the proposed approach can be accessed at Halilovic and Böttcher (2022).
No research studies have been identified that apply the approach of Class IV, which involves the combination of simplified models with DFO algorithms. It should be noted that other studies on GWHP optimization exist, focusing on aspects such as optimizing the components of a heat pump or determining optimal control strategies. However, these studies do not consider underground processes and are therefore outside the scope of this work. Furthermore, there are other research studies (e.g. Zhou and Zhou, 2009; Lo Russo and Civita, 2009; Gao et al., 2013) that address the optimal design or operation of GWHPs using methods such as scenario comparison or sensitivity analysis. However, these methods are not optimization methods, and as such, they are not further discussed in the present study. In the subsequent section, a comparative analysis is conducted between optimization approaches, i.e. the identified four classes.
3.2.2 Comparison of the optimization approaches
The primary factors for comparing optimization approaches, i.e. their respective classes, are the computational cost and applicability criteria. Figure 3 shows a qualitative comparison of these classes, considering two dimensions. The vertical axis represents the computational cost required to solve optimization problems with the approaches of the respective class. The computational cost of an approach is of major importance, since in practical planning procedures a relatively fast solution is required. Moreover, the computational cost increases proportionally with the number of optimization parameters (variables) and the size of the simulation domain. Consequently, approaches with high computational costs are limited to scenarios with a small number of optimization parameters and small domains. As the number of optimization parameters increases, inefficient approaches quickly become computationally impractical, even when using highperformance computers.
The horizontal axis represents the complexity (fidelity) of the groundwater simulation model used in these approaches. In the context of this study, the complexity of a simulation model refers to the level of detail in representing physical phenomena in the subsurface, such as the propagation of thermal plumes, that are relevant to the optimization problem under consideration. Assuming that the required input data, such as groundwater parameters, are available in sufficient quality, more complex simulation models are more accurate, i.e. closer to reality. However, it is important to recognize that data on groundwater parameters and conditions are often limited, which limits the use of complex models. Model complexity is essentially limited by the available data, making the use of highly complex models impractical in the absence of the necessary data. Nevertheless, simpler PDE simulation models, such as a 2D model with uniform groundwater conditions, are applicable even with restricted data availability and generally offer higher accuracy than analytical models with identical input data. Since simulation models are an integral component of optimization approaches, their complexity directly affects the applicability of the obtained optimization results. For instance, the results of an approach that uses a complex groundwater simulation model provided with highquality data can be applied in practice with greater confidence than the results of an approach based on less accurate models.
In the context of computational costs, two key aspects deserve attention: the convergence rate and the computational cost associated with the evaluation of each candidate solution (a unique combination of optimization variables). The former quantifies the number of optimization iterations required to reach the optimal solution, while the latter describes the runtime required for each model simulation used to evaluate the current candidate solution within the optimization iterations. In general, gradientbased algorithms significantly outperform derivativefree algorithms in terms of convergence rate and therefore it is recommended to use gradientbased algorithms when gradient information is readily available and can be obtained at a reasonable cost (Audet and Hare, 2017; Conn et al., 2009). As a result, Class II will almost always outperform Class IV, and Class I will outperform Class III, due to the use of gradientbased algorithms in the former classes (I and II) and derivativefree algorithms in the latter classes (III and IV). Another disadvantage of Classes III and IV is that derivativefree algorithms generally only find nearoptimal solutions and do not guarantee optimality (Audet and Hare, 2017). Furthermore, classes that use simplified simulation models (II and IV) commonly have lower computational costs than classes that use PDE models (I and III). This is a direct consequence of the computational costs associated with evaluating the simulation model during optimization iterations. Considering all of the above, a hierarchy of classes based on overall computational cost can be established. Class II entails the least computational cost, while Class III is the most computationally demanding. Classes I and IV fall somewhere in between, with Class IV usually outperforming Class I, although the specific problem characteristics (number of optimization variables, size of the simulation domain, etc.) can also influence this comparative performance. Despite the computationally demanding nature of Class III, this class is frequently used in the simulation community owing to its userfriendly nature. By coupling standard standalone simulation software with an existing implementation of a DFO algorithm, typically an evolutionary algorithm, users can develop and apply such approaches relatively quickly.
In terms of the complexity/fidelity of the simulation model used in optimization approaches, it is evident that the classes employing PDE models (I and III) outperform those employing simplified models (II and IV). The complexity of the simulation model directly influences the validity of the optimization results, thereby affecting the applicability of the corresponding classes. Consequently, it is reasonable to use approaches from different classes for different application scenarios. For instance, the classes with more complex models (I and III) are suitable for detailed planning of large GWHP systems, while the other two classes (II and IV) can be applied for initial assessments of potential negative interactions between neighboring systems or estimations of geothermal potential on a larger scale.
While there are a limited number of research studies (see Sect. 3.2.1) that address the optimization of GWHP systems, the research area remains insufficiently explored, which provides an opportunity to pose new research questions and develop novel optimization approaches. The existing approaches in this field have certain limitations and do not cover all relevant applications. For example, the approach of Class III presented in Park et al. (2020, 2021) is limited by the number of optimization variables it can efficiently handle, since the computational cost increases exponentially as the number of variables increases. Similarly, the only study (Halilovic et al., 2022a) using the approach of Class I does not cover all relevant aspects of GWHP optimization, such as optimizing the number of wells in a large GWHP system, optimizing pumping rates, or simultaneously optimizing pumping rates and well locations.
Class I (PDECO) seems to be the most promising among the four classes because it uses PDE simulation models and has lower computational costs compared to Class III. Here, the complexity level of the PDE model can be selected based on data availability, as discussed in Sect. 3.2.2. However, this class presents significant challenges due to its mathematical complexity and multidisciplinary nature. To overcome the challenges associated with developing new approaches within Class I and facilitate their further advancement, collaboration within multidisciplinary teams will be required in the future.
The limitations of the classes that use simplified simulation models (Class II and IV) are directly related to the limitations of the simulation models employed. Consequently, improving the simplified simulation models directly enhances the approaches within these classes. The main goal is to maintain the simulation models as fast and simple to evaluate while enhancing their closeness to reality. By further improving the accuracy of these simplified models, their scope can be extended to new applications, such as detailed design of large GWHP systems comprising multiple extraction and injection wells. Moreover, the simplified models are wellsuited for integration into energy system optimization models (ESOMs), where GWHP systems play an important role (Halilovic et al., 2022c). This is because the coupling of a numerical groundwater simulation with an ESOM is impractical and computationally demanding (Halilovic et al., 2022b). By using simplified models, the computational cost can be significantly reduced while still capturing the essential aspects of GWHP systems in the broader context of energy system optimization. This further enables the development of automated urban planning tools, which will increase sustainability.
Another important consideration in GWHP optimization is the inherent uncertainty associated with subsurface parameters and conditions. The complex nature of aquifers and the limited availability of measurement and monitoring data contribute to the presence of uncertainties (Gelhar et al., 1992). Incorporating these uncertainties into optimization approaches leads to stochastic programming problems (Birge and Louveaux, 2011), which constitute a separate field of optimization. The inherent stochastic nature of these problems significantly increases the complexity and computational cost compared to deterministic problems. Stochastic problems are often solved with modified deterministic optimization approaches or by using a deterministic equivalent of the stochastic problem (Hannah, 2015; Li and Grossmann, 2021). By minimizing the computational cost of deterministic optimization approaches, researchers can better address the challenges of stochastic problems and develop efficient approaches to such problems. Therefore, it is crucial for the deterministic optimization approaches discussed in this study to reduce their computational cost. Several strategies can be employed to reduce the computational cost of the existing approaches. For instance, the Class III approach (Park et al., 2020) may improve its efficiency by using modelbased algorithms instead of a genetic algorithm. Similarly, the Class I approach (Halilovic et al., 2022a) can improve its efficiency by using suitable (re)meshing techniques or by finetuning optimization algorithm parameters.
It is important to note that the classification and comparison of optimization approaches presented in Sect. 3.2 is not limited to the optimization of GWHP systems, but can be applied to any optimization problem where the underlying physical phenomena are described by PDEs. Moreover, the approaches developed for GWHP systems (see Sect. 3.2.1) and their future advancements or new approaches in this area can be extended to other applications. First, they can be extended to applications that share the same underlying physics, such as optimization of aquifer thermal energy storage (ATES) systems, calibration of numerical hydrothermal groundwater simulation models, and optimization of observation well placement. Second, these approaches can be extended to other areas of shallow geothermal energy, including optimization of vertical and horizontal closedloop shallow geothermal systems and optimization of borehole thermal energy storage (BTES) systems. Lastly, these approaches can be further extended to areas involving different physical phenomena, such as the optimization of wind farms or tidal power plants. Nevertheless, it is important to note that advances in these other areas, particularly in the area of shallow geothermal energy, can reciprocally contribute to the improvement of optimization approaches for GWHP systems. Namely, optimization approaches and principles used in other areas have the potential to be adapted and applied to GWHP systems.
This paper presents a comprehensive analysis and overview of approaches for optimizing the design and operation of GWHP systems. First, the optimization problems arising from this research and practice question were investigated, using a generalized problem as a basis. Then, optimization approaches were identified and compared, and a novel classification of the approaches is proposed. The identified approaches were divided into four distinct classes based on the type of groundwater simulation model used (PDEbased or simplified models) and the optimization algorithm applied (gradientbased or derivativefree). Finally, the paper includes a thorough review of the existing approaches in the literature, highlighting their limitations and outlining opportunities for future improvements.
Based on the analysis performed, several conclusions can be drawn:

Optimization approaches that rely on gradientbased optimization algorithms are preferable, as they consistently outperform derivativefree algorithms.

The choice of a simulation model used in an optimization approach has a significant impact on its applicability. For example, approaches using PDE models are more suitable for detailed design of largescale GWHPs, while simplified models offer practical advantages for assessing the geothermal potential of large areas. However, it is important to note that the degree of model complexity is limited by the availability of hydrogeological data.

The existing research on GWHP optimization is limited, with only a few studies addressing this topic.

Existing approaches have certain limitations and do not cover all relevant applications and research questions in GWHP optimization. One of the main limitations is the high computational cost, which limits the number of optimization parameters and the size of the simulation domain that can be effectively considered. In addition, some approaches are limited in applicability due to the use of simplified groundwater simulation models. Moreover, applications such as optimizing the number and placement of wells in large GWHP systems or simultaneous optimization of pumping rates and well placements remain unexplored. Consequently, there is an ongoing need to develop new and improve existing approaches to address these limitations and fill the research gaps.

The efficient optimization approaches developed for GWHP systems have the potential to be extended to other shallow geothermal applications as well as to other optimization problems where the underlying physical phenomena are described by PDEs. At the same time, approaches from other areas can be adapted and used for GWHP optimization in the future.
Finally, this study can provide a valuable foundation for researchers and practitioners involved in the management and optimization of shallow geothermal energy systems. In particular, it provides valuable insights and recommendations for the application and development of optimization approaches for GWHP systems. Despite its challenging nature, optimization of GWHP systems is of utmost importance to improve thermal management of groundwater and to unlock the full potential and attractiveness of GWHP technology.
No data sets were used in this article.
SH: Conceptualization, Writing – original draft, Methodology, Investigation, Visualization. FB: Writing – reviewing and editing. KZ: Writing – reviewing and editing. TH: Writing – reviewing and editing, Supervision.
The contact author has declared that none of the authors has any competing interests.
Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.
This article is part of the special issue “European Geosciences Union General Assembly 2023, EGU Division Energy, Resources & Environment (ERE)”. It is a result of the EGU General Assembly 2023, Vienna, Austria, 23–28 April 2023.
No external funding was received for conducting this study. We would like to thank Jannis Epting for reviewing the paper and providing valuable comments that improved its quality.
This work was supported by the Technical University of Munich (TUM) in the framework of the Open Access Publishing Program.
This paper was edited by Michael Kühn and Giorgia Stasi, and reviewed by Jannis Epting and one anonymous referee.
Audet, C. and Hare, W.: Derivativefree and blackbox optimization, Springer, ISBN 9783319689135, https://doi.org/10.1007/9783319689135, 2017. a, b, c
Birge, J. R. and Louveaux, F.: Introduction to stochastic programming, Springer Science & Business Media, ISBN 9781461402374, https://doi.org/10.1007/9781461402374, 2011. a
Blum, P., Menberg, K., Koch, F., Benz, S. A., Tissen, C., Hemmerle, H., and Bayer, P.: Is thermal use of groundwater a pollution?, J. Contam. Hydrol., 239, 103791, https://doi.org/10.1016/j.jconhyd.2021.103791, 2021. a, b
Bonnans, J.F., Gilbert, J. C., Lemaréchal, C., and Sagastizábal, C. A.: Numerical optimization: theoretical and practical aspects, Springer Science & Business Media, ISBN 9783540354475, https://doi.org/10.1007/9783540354475, 2006. a
Böttcher, F. and Zosseder, K.: Thermal influences on groundwater in urban environments – A multivariate statistical analysis of the subsurface heat island effect in Munich, Sci. Total Environ., 810, 152193, https://doi.org/10.1016/j.scitotenv.2021.152193, 2021. a
Böttcher, F., Casasso, A., Götzl, G., and Zosseder, K.: TAP – Thermal aquifer Potential: A quantitative method to assess the spatial potential for the thermal use of groundwater, Renew. Energy, 142, 85–95, https://doi.org/10.1016/j.renene.2019.04.086, 2019. a
BozorgHaddad, O., Solgi, M., and Loáiciga, H. A.: Metaheuristic and evolutionary algorithms for engineering optimization, John Wiley & Sons, ISBN 9781119387053, 2017. a
Conn, A. R., Gould, N. I., and Toint, P. L.: Trust region methods, SIAM, ISBN 9780898714609, 2000. a
Conn, A. R., Scheinberg, K., and Vicente, L. N.: Introduction to derivativefree optimization, SIAM, ISBN 9780898716689, 2009. a, b
Davis, K., Leiteritz, R., Pflüger, D., and Schulte, M.: Deep learning based surrogate modeling for thermal plume prediction of groundwater heat pumps, arXiv [preprint], arXiv:2302.08199, https://doi.org/10.48550/arXiv.2302.08199, 2023. a
Diersch, H.J. G.: FEFLOW, Springer, Berlin, Heidelberg, ISBN 9783642387388, https://doi.org/10.1007/9783642387395, 2014. a
Epting, J. and Huggenberger, P.: Unraveling the heat island effect observed in urban groundwater bodies – Definition of a potential natural state, J. Hydrol., 501, 193–204, https://doi.org/10.1016/j.jhydrol.2013.08.002, 2013. a
Epting, J., Böttcher, F., Mueller, M. H., GarcíaGil, A., Zosseder, K., and Huggenberger, P.: Cityscale solutions for the energy use of shallow urban subsurface resources – Bridging the gap between theoretical and technical potentials, Renew. Energy, 147, 751–763, https://doi.org/10.1016/j.renene.2019.09.021, 2020. a
Florides, G. and Kalogirou, S.: Ground heat exchangers – A review of systems,models and applications, Renew. Energy, 32, 2461–2478, https://doi.org/10.1016/j.renene.2006.12.014, 2007. a
Gao, Q., Zhou, X.Z., Jiang, Y., Chen, X.L., and Yan, Y.Y.: Numerical simulation of the thermal interaction between pumping and injecting well groups, Appl. Therm. Eng., 51, 10–19, https://doi.org/10.1016/j.applthermaleng.2012.09.017, 2013. a
GarcíaGil, A., Mejías Moreno, M., Garrido Schneider, E., Marazuela, M. Á., Abesser, C., Mateo Lázaro, J., and Sánchez Navarro, J. Á.: Nested Shallow Geothermal Systems, Sustainability, 12, 5152, https://doi.org/10.3390/su12125152, 2020. a
García Gil, A., Garrido Schneider, E. A., Mejías Moreno, M., and Santamarta Cerezal, J. C.: Shallow Geothermal Energy, Springer, https://doi.org/10.1007/9783030922580, 2022. a, b
Gelhar, L. W., Welty, C., and Rehfeldt, K. R.: A Critical Review of Data on FieldScale Dispersion in Aquifers, Water Resour. Res., 28, 1955–1974, https://doi.org/10.1029/92WR00607, 1992. a
Guimerà, J., Ortuño, F., Ruiz, E., Delos, A., and PérezParicio, A.: Influence of groundsource heat pumps on groundwater, in: Conference Proceedings: European Geothermal Congress, 30 May–1 June 2007, Unterhaching, Germany, https://www.geothermalenergy.org/pdf/IGAstandard/EGC/2007/250.pdf (last access: 18 December 2023), 2007. a
Hähnlein, S., MolinaGiraldo, N., Blum, P., Bayer, P., and Grathwohl, P.: Ausbreitung von Kältefahnen im Grundwasser bei Erdwärmesonden, Grundwasser, 15, 123–133, https://doi.org/10.1007/s007670090125x, 2010. a
Halilovic, S. and Böttcher, F.: Optimization of GWHP well layouts using analytic models, Zenodo [code], https://doi.org/10.5281/zenodo.7230875, 2022. a
Halilovic, S., Böttcher, F., Kramer, S. C., Piggott, M. D., Zosseder, K., and Hamacher, T.: Well layout optimization for groundwater heat pump systems using the adjoint approach, Energ. Convers. and Manag., 268, 116033, https://doi.org/10.1016/j.enconman.2022.116033, 2022a. a, b, c, d
Halilovic, S., Odersky, L., Böttcher, F., Davis, K., Schulte, M., Zosseder, K., and Hamacher, T.: Optimization of an Energy System Model Coupled with a Numerical Hydrothermal Groundwater Simulation, in: Mapping the Energy Future – Voyage in Uncharted Territory, 43rd IAEE International Conference, 31 July–3 August 2022, International Association for Energy Economics, http://www.iaee.org/proceedings/article/17725 (last access: 18 December 2023), 2022b. a
Halilovic, S., Odersky, L., and Hamacher, T.: Integration of groundwater heat pumps into energy system optimization models, Energy, 238, 121607, https://doi.org/10.1016/j.energy.2021.121607, 2022c. a
Halilovic, S., Böttcher, F., Zosseder, K., and Hamacher, T.: Optimizing the spatial arrangement of groundwater heat pumps and their well locations, Renew. Energy, 217, 119148, https://doi.org/10.1016/j.renene.2023.119148, 2023. a
Hammond, G., Lichtner, P., Lu, C., and Mills, R. T.: PFLOTRAN: Reactive flow & transport code for use on laptops to leadershipclass supercomputers, Groundwater React. Trans. Models, 5, 141–159, 2012. a
Hannah, L. A.: Stochastic optimization, Int. Encycloped. Social Behav. Sci., 2, 473–481, 2015. a
Hinze, M., Pinnau, R., Ulbrich, M., and Ulbrich, S.: Optimization with PDE constraints, in: vol. 23, Springer Science & Business Media, ISBN 9781402088391, https://doi.org/10.1007/9781402088391, 2008. a, b
Kim, J. and Nam, Y.: A Numerical Study on System Performance of Groundwater Heat Pumps, Energies, 9, 4, https://doi.org/10.3390/en9010004, 2016. a
Kinzelbach, W.: Numerische Methoden zur Modellierung des Transports von Schadstoffen im Grundwasser, Oldenbourg, ISBN 9783486263473, 1992. a, b
Le Digabel, S.: Algorithm 909: NOMAD: Nonlinear optimization with the MADS algorithm, ACM Trans. Math. Softw., 37, 1–15, https://doi.org/10.1145/1916461.1916468, 2011. a
Leiteritz, R., Davis, K., Schulte, M., and Pflüger, D.: A Deep Learning Approach for Thermal Plume Prediction of Groundwater Heat Pumps, arXiv [preprint], arXiv:2203.14961, https://doi.org/10.48550/arXiv.2203.14961, 2022. a
Li, C. and Grossmann, I. E.: A review of stochastic programming methods for optimization of process systems under uncertainty, Front. Chem. Eng., 2, 34, https://doi.org/10.3389/fceng.2020.622241, 2021. a
Logg, A., Mardal, K., and Wells, G. N.: Automated Solution of Differential Equations by the Finite Element Method, Springer, https://doi.org/10.1007/9783642230998, 2012. a
Lo Russo, S. and Civita, M. V.: Openloop groundwater heat pumps development for large buildings: A case study, Geothermics, 38, 335–345, https://doi.org/10.1016/j.geothermics.2008.12.009, 2009. a
Marler, R. T. and Arora, J. S.: The weighted sum method for multiobjective optimization: new insights, Struct. Multidiscip. Optimiz., 41, 853–862, https://doi.org/10.1007/s0015800904607, 2010. a
Menberg, K., Bayer, P., Zosseder, K., Rumohr, S., and Blum, P.: Subsurface urban heat islands in German cities, Sci. Total Environ., 442, 123–133, https://doi.org/10.1016/j.scitotenv.2012.10.043, 2013. a
Naumann, U.: The Art of Differentiating Computer Programs, Society for Industrial and Applied Mathematics, https://doi.org/10.1137/1.9781611972078, 2011. a
Nocedal, J. and Wright, S. J.: Numerical optimization, Springer, ISBN 9780387227429, https://doi.org/10.1007/b98874, 1999. a
Ohmer, M., Klester, A., Kissinger, A., Mirbach, S., Class, H., Schneider, M., Lindenlaub, M., Bauer, M., Liesch, T., Menberg, K., and Blum, P.: Berechnung von Temperaturfahnen im Grundwasser mit analytischen und numerischen Modellen, Grundwasser, 27, 113–129, https://doi.org/10.1007/s00767022005092, 2022. a
Park, D., Lee, E., Kaown, D., Lee, S.S., and Lee, K.K.: Determination of optimal well locations and pumping/injection rates for groundwater heat pump system, Geothermics, 92, 102050, https://doi.org/10.1016/j.geothermics.2021.102050, 2021. a, b
Park, D. K., Kaown, D., and Lee, K.K.: Development of a simulationoptimization model for sustainable operation of groundwater heat pump system, Renew. Energy, 145, 585–595, https://doi.org/10.1016/j.renene.2019.06.039, 2020. a, b, c
Perego, R., Dalla Santa, G., Galgaro, A., and Pera, S.: Intensive thermal exploitation from closed and open shallow geothermal systems at urban scale: unmanaged conflicts and potential synergies, Geothermics, 103, 102417, https://doi.org/10.1016/j.geothermics.2022.102417, 2022. a
Pophillat, W., Attard, G., Bayer, P., HechtMéndez, J., and Blum, P.: Analytical solutions for predicting thermal plumes of groundwater heat pump systems, Renew. Energy, 147, 2696–2707, https://doi.org/10.1016/j.renene.2018.07.148, 2020. a, b
Raissi, M., Perdikaris, P., and Karniadakis, G. E.: Physics Informed Deep Learning (Part I): Datadriven Solutions of Nonlinear Partial Differential Equations, arXiv [preprint], arXiv:1711.10561, https://doi.org/10.48550/arXiv.1711.10561, 2017. a
Rathgeber, F., Ham, D. A., Mitchell, L., Lange, M., Luporini, F., Mcrae, A. T. T., Bercea, G.T., Markall, G. R., and Kelly, P. H. J.: Firedrake, ACM Trans. Math. Softw., 43, 1–27, https://doi.org/10.1145/2998441, 2017. a
Russo, S. L., Taddia, G., and Verda, V.: Development of the thermally affected zone (TAZ) around a groundwater heat pump (GWHP) system: A sensitivity analysis, Geothermics, 43, 66–74, https://doi.org/10.1016/j.geothermics.2012.02.001, 2012. a
Russo, S. L., Taddia, G., Gnavi, L., and Verda, V.: Neural network approach to prediction of temperatures around groundwater heat pump systems, Hydrogeol. J., 22, 205–216, https://doi.org/10.1007/s1004001310722, 2014. a
Schrijver, A.: Theory of linear and integer programming, John Wiley & Sons, ISBN 9780471982326, 1998. a
Stauffer, F., Bayer, P., Blum, P., Giraldo, N. M., and Kinzelbach, W.: Thermal use of shallow groundwater, CRC Press, Boca Raton, Florida, ISBN 9781466560192, 2014. a
Tröltzsch, F.: Optimal control of partial differential equations: theory, methods, and applications, in: Graduate Studies in Mathematics, Band 112, American Mathematical Society, ISBN 9780821849040, 2010. a
Zhou, Y.z. and Zhou, Z.f.: Simulation of Thermal Transport in Aquifer: A GWHP System in Chengdu, China, J. Hydrodynam., 21, 647–657, https://doi.org/10.1016/S10016058(08)601961, 2009. a