TAP - Thermal aquifer Potential: A quantitative method to assess the spatial potential for the thermal use of groundwater

This paper proposes a method to assess the potential for thermal use of groundwater and its integration in spatial energy planning. The procedure can be adapted to local regulatory and operational limits, thus estimating legally and technically achievable ﬂ ow rates and subsequently, the thermal power that can be exchanged with the aquifer through a well doublet. The constraints applied to ﬂ ow rates are i Þ a drawdown threshold in the extraction well, ii Þ a limit for the groundwater rise in the injection well and iii Þ a threshold to avoid the hydraulic breakthrough between the two wells. For the spatial assessment, the hydraulic in ﬂ uence on neighbouring well doublets is simulated with the maximum ﬂ ow rates before the hydraulic breakthrough occurs. The Thermal Aquifer Potential (TAP) method combines mathematical relations derived through non-linear regression analysis using results from numerical parameter studies. A demonstration of the TAP method is provided with the potential assessment in Munich, Germany. The results are compared with monitoring data from existing open-loop systems, which prove that conservative peak extraction estimates are achieved.


Introduction
The European Union (EU) has set ambitious goals to reduce greenhouse gas emissions, increase the share of renewable energy (to 32% in 2030) and improve its energy efficiency [1]. Given that heating and cooling accounts for 50% of the EU's energy needs, any efficiency improvement in this sector will contribute significantly to reaching those targets. In particular, the use of shallow geothermal energy (SGE) mitigates greenhouse gas emissions by reducing the primary energy consumption and hence increasing energy efficiency [2e5]. Thus, geothermal heat pumps can provide a strategic contribution in the de-carbonization of the heating and cooling sector [6e8].
SGE can be used through open loop systems, which exchange heat with groundwater by means of a heat pump, commonly known as groundwater heat pumps (GWHPs). Separate wells are installed as a doublet to extract water from a shallow aquifer and re-inject it after thermal exchange [9]. In consequence, the availability of groundwater in a sufficient quantity, quality, temperature and depth is the requirement for a sustainable operation of open-loop systems [10e12]. This makes the thermal use of groundwater a renewable, but spatially limited resource with a greatly varying availability [13,14]. Detailed knowledge about the local hydrogeology and the resulting technical potential is crucial, for the installation of efficient open loop systems [15]. If not available, cost-intensive exploration, like the drilling of observation wells, the conduction of well tests and time-intensive data inquiry can hamper the diffusion of GWHP systems [16]. Therefore, spatial potential evaluations for the thermal use of groundwater are key to identifying suitable areas and to fostering its use [17].
The legislative procedures for licencing SGE systems commonly follow the "first come, first served" principle [18,19]. However, a development of procedures towards active resource management would support authorities in the implementation of energy strategies [20,21]. Spatial potential evaluations serve as a data basis on which those management concepts can be built and on which approval decisions can be made [22].
The already intensive thermal use of groundwater in suitable areas, like Munich, shows how SGE resources have become relevant and should be taken into account in any energy development scenario [23]. In detail, assessment concepts have to provide quantitative extraction values at the smallest scale of urban energy planning, which assure a sustainable operation within legal and technological boundaries [24]. Previous studies, listed in Table 1, already identified the most relevant issues that constrain the potential of GWHPs. However, no study included all of them, nor did they use a uniform numerical or analytic method to perform a potential assessment per plot of land. Further, the defined constraints have to properly address the relevant issues on the specific scale of assessment. For example, on a single plot multiple well doublets would rather be aligned perpendicular to the groundwater flow to avoid thermal interference. Thus, the influence of doublets to the side will be more relevant than the extent of thermal anomalies downstream. Concepts for closed loop systems are already capable of detailed quantitative potential estimation with an application in energy planning [25e30], whereas the open loop potential has not yet been assessed with all identified requirements.
This paper presents a quantitative method to assess the technical potential of thermal groundwater use for an application in spatial energy planning. The method integrates the relevant regulatory and operational constraints for an estimation of technologically achievable flow rates in well doublets and hence of the thermal power, which can be extracted at a specific site. The following section 2 presents the assessment concept and the methodology. Section 3 addresses an exemplary application in Munich with a comparison between the TAP results and monitoring data in section 3.3. Subsequently, the derived method as well as its assumptions and limitations are discussed in section 4 with reference to the studies introduced in Table 1.

The TAP methodology e a legislative-operational potential
In the following, the TAP method is presented. First, the general approach is introduced to define the scope of the potential assessment (section 2.1). Second, the numerical approach (section 2.2) is explained to provide the basis for the mathematical formulation of the technical flow rate equations (section 2.3) and finally, the hydraulic footprint of a well doublet is studied to deliver values on appropriate well spacing (section 2.4).

General approach of the potential assessment
The major challenge in geothermal potential evaluation is the deduction from a theoretical (natural) potential to a more applicable estimation, which includes technological, legal and economical barriers [39,40]. In addition, a potential estimation should have a high general validity to make it practically useful and sufficiently conservative. Since costs and subsidies are temporally and locally more variable than e.g. regulatory elements [41], an economic potential assessment is not in the focus of this study.
The approach integrates constraints within several steps. First, the legislative framework is reviewed to obtain the spatial restrictions of water protection areas and drilling limitations, as well as system related restrictions, like temperature, abstraction, drawdown and injection limits. Second, operational limits are defined with safety margins according to common practice, e.g. minimum injection of 4 C-cold water. All identified issues limiting the technical potential, as summarised in Table 2, are merged to evaluate a sustainable (technical) flow rate potential. Subsequently, it can be used to calculate the extractable thermal power with approvable temperature differences as well as absolute minimum and maximum injection temperature limits, i.e. 4 C and 20 C in the case of Munich. With the assessment results, the potential of thermal groundwater use can be integrated into energy planning at the building scale. Since the initial potential estimates follow the previously mentioned "first come, first served" principle, possible negative influences on neighbours are not yet considered. Instead, a maximum extraction per plot or raster cell is provided within legal boundaries (cf. section 3.2). This strategy is also proposed for the first stage of licence application in recent SGE management concepts [42]. The following section 2.2 describes the modelling approach to derive a functional dependence for the respective constraints.

Numerical modelling concept
For a comprehensive consideration of influential parameters, a numerical approach under the same simplifying assumptions is used. Therefore, idealised finite-element models are built, to determine the relevant influences on achievable groundwater flow rates in parameter studies. The simulations that vary the significant parameters against each other are conducted in idealised boxmodels with the following characteristics:   1a shows the design of the numerical 2D-box-model and its boundary conditions. The model represents all material and geometric properties in an isotropic and constant way. The extraction and injection wells are located in the centre of the model, parallel to the groundwater flow at varying distances. The variation ranges of significant parameters are set according to their distribution in the study area, as shown for Munich in Table 3.
For an application in other areas, where certain parameters might not be spatially known, ranges can also be derived from Table 1 Studies on the quantitative potential assessment of thermal groundwater use with applied constraints. *(share of saturated aquifer thickness).  [33] 25%* and 100 L/s (Thiem) e e e G€ otzl et al. [34] Mean and low mean (Thiem) e e e Javandel & Tsang [35] 100%* (Theis) e e complete capture zone Muñoz et al. [36] 100%* (Logan) e e e PDGN [37] e e e thermal anomaly (GED) Urich et al. [38] e 5% inter-flow & 0.5 K feedback (HST3D) e thermal anomaly (HST3D) representative local measurements or literature values. Scenarios for the variations are given by the Cartesian product, which is composed of the respective set of significant parameters. Subsequently, the simulation results serve as a basis for a non-linear regression analysis to fit the mathematical formulations of the TAP method, as presented in the following sections [44]. With an estimation of regression functions, an integration of the assessment procedure in existing GIS-workflows is straightforward, e.g. in energy planning (cf. section 3.2).

Mathematical formulation of the technical flow rate
The quantitative assessment of technical flow rates copes with the main design constraints of well doublets. The addressed issues are iÞ the drawdown in the extraction well, which can lead to the depletion of the aquifer (section 2.3.1); iiÞ the possible return of reinjected water to the extraction well, with the consequence of thermal recycling and decreasing efficiency (section 2.3.2); iiiÞ the level increase at the reinjection well, which can result in groundwater flooding (section 2.3.3). The flow rate estimation is based on the characteristics of one entire well doublet including their mutual interference.

Maximum drawdown constraint
The estimation of an extractable well yield is fundamental in the quantitative evaluation of groundwater resources [45]. The flow rate is generally defined based on a maximum acceptable hydraulic head decline (drawdown) in the well. Considering fully developed wells, a commonly agreed drawdown limit is one third of the saturated aquifer thickness as a technical threshold for sustainable operation [31]. However, it can also be adjusted in the simulation procedure.
The parameter study was performed with the significant parameters, hydraulic conductivity and aquifer thickness within the ranges shown in Table 3. Simulation records from 36 (6,6, cf. Table 3) scenarios were used to fit the constant regression coefficient to 0.195 in equation (1). where: : saturated aquifer thickness Equation (1) calculates the maximum flow rate at the drawdown threshold, in dependence of saturated aquifer thickness and hydraulic conductivity. The maximum extractable flow rate is therefore proportional to the transmissivity and the saturated thickness, which is in line with common approaches such as proposed by Misstear et al. [46]. The scatter-plot in Fig. 2a displays the   fit of the regression estimates compared with the simulation results. In this setting, the maximum flow rate at the drawdown constraint is determined with an influence of the injection well at the limit conditions before a hydraulic breakthrough occurs. Fig. 1b displays a cross section of the hydraulic head through the two wells, where this influence becomes visible. Considering only a single extraction well would result in a larger groundwater decline and therefore in lower pumping rates, as a product of a missing hydraulic support from the injection well [32]. To quantify this influence, the parameter study was repeated with only one extraction well. The results showed the same functional dependency, but consistently 18% lower pumping rates at this drawdown threshold.

Hydraulic breakthrough constraint
Apart from a maximum abstraction estimate, additional measures have to be applied to guarantee the sustainable operation of a well doublet. The injection well, located downstream of the extraction well, returns the thermally altered water to the aquifer. Therefore, extensive pumping in the extraction well can result in a hydraulic breakthrough, where the injected water re-enters the extraction well. As a consequence, thermal recycling will eventually occur and the efficiency of the system diminishes [47]. To prevent thermal recycling also for continuous annual loads, the maximum pumping rate without a hydraulic breakthrough is simulated. This leads to conservative flow rate values, where a thermal interaction between the wells is excluded [48e50].
The numerical parameter study to determine the hydraulic breakthrough constraint is also conducted with the model, shown in Fig. 1a. The maximised flow rate of a well doublet with no hydraulic feedback is determined by the Darcy velocity vector as it approaches zero in the stagnation point. In the simulation concept, the flow rate of the well doublet changes iteratively, until the Darcy velocity in the middle of the wells is between ± 0.001 [m/d]. The resulting hydraulic head is displayed in Fig. 1b and only pumping rates causing this limit state are considered in the regression analysis.
In the box-model, the significant parameters, hydraulic conductivity, hydraulic gradient and aquifer thickness are varied within the ranges shown in Table 3. In addition, the parameter study is extended by a set of five well distances in the range from 10 m to the distance where also a drawdown of 1 = 3 aquifer thickness is exceeded. where: The resulting record of 1080 (6 3 ,5, cf. Table 3) simulations is used to fit the regression coefficient (1.96) in equation (6). Fig. 2 b displays the simulation results against the flow rate predictions. (1) and (2) can be solved for x w . As a result, equation (3) calculates the maximum well distance for a well doublet operating simultaneously at the drawdown and the breakthrough threshold. Thus, equation (3) also calculates the maximum well distance used for the five x w scenarios.
: saturated aquifer thickness x w ½m : extraction to injection well distance where: b ½m : saturated aquifer thickness x w ½m : extraction to injection well distance i ½ À : hydraulic gradient

Injection constraint
In areas where the depth to the groundwater surface is low, the re-injection of water should not cause a rise in the groundwater level, which might result in surface or basement flooding. Therefore, the injection capacity of the aquifer is an important factor in sensitive areas and a sufficient safety margin should be applied [32,51]. To prevent flooding problems with a standard well design, an injection limit is included in the assessment of technical flow rate potentials.
For the simulations, the same model is used (cf. Fig. 1a). The parameter study of 216 (6 3 , cf. Table 3) scenarios includes the variation of the hydraulic conductivity, hydraulic gradient and aquifer thickness. The regression coefficients of equation (4) (0.789 and 29.9) are fitted, with the record of the induced groundwater level rise at the injection well and the related flow rates. The achieved accuracy of the regression function's predictions is presented in Fig. 2c. Equation (4) calculates the flow rate, which can be injected until a rise from the natural to the maximum groundwater level to the maximum allowed level (z max ) is reached. The scenarios are also simulated with only one infiltration well in order to quantify the hydraulic influence of the extraction well. Because the hydraulic influence of the extraction well diminishes the infiltration cone, the well doublet is able to supply higher pumping rates of around 10% at the same level of groundwater rise.
Finally, the three flow rates ( _ V drawdown , _ V injection , _ V breakthrough ) are compared. The minimum flow rate is the constraining factor and sets the technical flow rate for a well doublet ( _ V tech ).

Spatial flow rate potential
Based on the "internal" distance of the well doublet, a suitable "external" distance to hypothetically neighbouring doublets can be derived with the following results. This step transfers the technical into a spatial potential, because the hydraulic footprint of a well doublet is considered and the water budget in each footprint area is balanced.
The relation between external to internal well distance and the percentage of cycled water in a well doublet (inter-flow) is determined in numerical simulations. For this purpose, a slight modification of the model presented in Fig. 1a is used. The area is enlarged to 4000 m edge length to host multiple well doublets, which are lined up perpendicular to the groundwater flow direction in the model centre for 2000 m (cf. Fig. 3a). This model setup was constructed with four "internal" well distances (10, 25, 50, 100 m).
In a sensitivity analysis, external to internal well distance ratios from 1.0 to 3.0 are simulated in 0.2 steps (cf. Fig. 3b). The sensitivity is also analysed for different "internal" well distances and for the influential parameters shown in Table 3, summing up to 242 (ð6,3 þ 4Þ,11) cases. The flow rate ( _ V breakthrough ) for each scenario is calculated from equation (2) and it is kept constant for a variation of the external well distances.
The closer the well doublets are spaced, the more water is circulated between injection and extraction well, because the width of upstream capture zones and downstream release fronts is narrowed. Fig. 3b shows the resulting inter-flow percentages for the ratio range of the sensitivity analysis. Since _ V breakthrough is proportional to the respective parametrisation, it can be observed that the remaining influence on the inter-flow is caused only by the ratio of well distances. In consequence, a ratio can be defined, depending on the required conservativeness of spatial estimates. The ratio can further be used to calculate the cell size of the doublet's hydraulic footprint from a certain internal well distance. Fig. 3a displays the resulting flow pattern for the hydraulic footprint at an internal well distance of 10 m at a 20 m resolution. This corresponds to a well distance ratio of 2.0 and would limit the interflow to moderate values from 4.2 to 4.9% at maximum spacing density.

The case study: Munich
The proposed TAP method is applied in Munich and enables the city's urban energy planners and decisional authorities to consider the potential of thermal groundwater use.

Hydrogeological setting
The city of Munich is located on the so-called Munich gravel plain. Since the Pleistocene period, fluvio-glacial and fluvial sediments have accumulated on a Tertiary palaeo-surface and built a widespread alluvial plain. A large unconfined groundwater body lies in the Quaternary and Holocene sandy gravels, where the average hydraulic conductivity is 5,10 À3 m=s. The underlying Upper Freshwater Molasse as topmost Tertiary formation consists of heterogeneously changing layers of predominantly fine sand, silt and clay. Therefore, the Tertiary sediments have a hydraulic conductivity mostly below 1,10 À6 m=s and the Tertiary palaeo-surface normally acts as the shallow aquifer bottom [52]. Due to channel structures carved into the Tertiary deposits by the former drainage system, the saturated groundwater thickness varies largely in the city area [53]. In consequence, a robust spatial interpolation of the aquifer thickness requires a very detailed spatial knowledge of the Quaternary/Tertiary boundary and the groundwater table.
To acquire this information, over 48,000 measurement points from borehole logs and outcrop observations have been reviewed and checked for plausibility in the GEothermal POtential (GEPO) project from 2012 to 2015 [54]. With this data basis the Quaternary/ Tertiary boundary was interpolated throughout the hole Munich gravel plain. In addition, an elaborate measurement campaign containing over 6000 groundwater wells was carried out in April 2014. Fig. 4 (top) displays the saturated thickness of Quaternary groundwater. Two major channel structures in North-South orientation can be observed, from Sendling-Westpark (21) to Moosach (23) and from Pasing-Obermenzing (07) to Allach-Untermenzing (10). Here the aquifer thickness is constantly above 10 m, whereas areas with no groundwater indicate that the Tertiary surface is elevated above the surrounding groundwater level. This is the case along the Isar and at some Tertiary hills, such as the "Aubinger Lohe" (22), which are elevated above ground as inliers. In consequence, the aquifer thickness varies largely and the potential for thermal use of groundwater is very heterogeneous. The potential is also influenced by the distribution of depth to water, which drives the economic accessibility of the resource or may hinder a sustainable groundwater re-injection. The depth to water increases towards the south and is very low in the northern districts of Bogenhausen (13), Feldmoching-Hasenbergl (24) and Aubing-Lochhausen-Langwied (22) (cf. Fig. 4 bottom). Furthermore, the aquifer is already intensively used from over 2600 users in 2017. This situation highlights the need for a spatial potential assessment, before an integration of the resource in energy planning can be elaborated.

Assessment of the thermal groundwater potential with the TAP method
As presented in section 2.1, the first assessment step is the consideration of legislative and operational constraints valid for Munich. Initially, the maximum approvable drilling depth, i.e. the Quaternary/Tertiary layer boundary, and the drinking water protection areas are reviewed. In the next step, the maximum flow rate values area calculated according to the constrains for drawdown, injection and hydraulic breakthrough (cf. section 2.3). Fig. 5 shows the resulting technical flow rates for well distances of 100 m (top) and 10 m (bottom). The resolution of the potential maps is set to 200 m and 20 m respectively, to consider the hydraulic footprint of the well doublets. As presented in Fig. 3b, this limits the inter-flow, i.e. the share of circulated water in a doublet, to below 5%, if systems in all neighbouring raster cells would be realised. Two different kinds of hatching indicate the areas where the constraint on groundwater drawdown (up to one third of the saturated aquifer thickness) or rise (up to 0.5 m below ground surface) applies. However, for both well distances the hydraulic breakthrough is the prevailing limiting factor in most of the area. Only in Schwabing-Freimann (area no. 12 in Fig. 4 top) and around the river Isar, the drawdown constraint is stronger than the one on hydraulic breakthrough, due to the low aquifer thickness. The injection limits the potential flow rate only in the north-western part such as Aubing-Lochhausen-Langwied (area no. 22 in Fig. 4 top), where the depth to water table can already be smaller than the safety threshold of 0.5 m. Increasing the well distance to 100 m leads to a ten-fold increase of the related flow rate ( _ V breakthrough ), thus amplifying the effect of the other constrains. The drawdown becomes the limiting factor in larger areas of low saturated aquifer thickness west of the Isar and in areas of a low hydraulic gradient (cf. Fig. 4 top). With an elevated _ V breakthrough limit, also the injection constraint is active in a larger area of low depth's to water, especially in Aubing-Lochhausen-Langwied (22). However, on the properties of these two houses well distances larger than 10 m could be accommodated, although they are way below 100 m, i.e. the other hypothetical well distance considered in Fig. 5 (top). The flexibility of TAP also allows to calculate the openloop geothermal potential of every single plot of land, considering the best practice of installing extraction and injection wells aligned with the groundwater flow. Fig. 6 b shows an example of this plotwise assessment, with the available space for a well doublet. Buffer distances provided by regulations are included, such as the 3 m minimum distance from neighbouring plots or buildings prescribed in Bavaria [55]. Plots are marked in red if the minimum prescribed well distance of 10 m is not feasible. With the plot-wise approach, the TAP method meets all requirements for an application in spatial energy planning on the building scale. Since only one plot with one well doublet is considered at a time, the flow rate results cannot be aggregated as mutual hydraulic and thermal interferences are not taken into account. If the required well distance cannot be installed, a comparative evaluation of the three constrained flow rates can indicate the feasibility of a specific thermal short circuit analysis.

Comparison of TAP results with existing installations
In Bavaria, shallow geothermal applications with an approved annual extraction of more than 100,000 m 3 have to report the monthly extracted volumes to the approving authority. Hence, the maximum monthly flow rates in 2015 are compared to the technical flow rate estimates with a validation data set of 97 well doublets. Flow rates are calculated for the drawdown limit from equation (1) and the hydraulic breakthrough limit from equation (2). The values of relevant hydrogeological parameters are averaged within a circle around the wells. The radius of the circle is calculated with respect to the reported pumping rates after Todd & Mays [56]. In addition, the well doublets are considered to be parallel to the regional groundwater flow and horizontal filter wells are excluded. Fig. 7 shows the cumulative distribution of the over-and undershooting systems compared to the estimates. The abscissa shows the ratio of maximum monthly means and estimated flow  rates, thus systems with a value higher than one extract more than expected. 10% of the reported extractions exceed the estimated flow rate at the drawdown limit and 26% overshoot the rate at the hydraulic breakthrough limit (cf. dashed lines in Fig. 7).
Two characteristics of the monitoring data should be highlighted. Although only large systems have been evaluated, they do not necessarily use the maximum potential at their site. Additionally, the full-load operation time of the peak month is not known. Apart from that, thermal recycling will need a certain duration of excessive pumping to develop. Especially for the examined doublets with a mean distance of 140 m. In consequence, mean values are more appropriate to evaluate a breakthrough risk than short full-load extractions.
From Fig. 7, we can observe that the majority of systems do not exceed the established drawdown and hydraulic breakthrough limits. No severe overshooting is present and at five well doublets, which exceed the estimated hydraulic breakthrough threshold by 1.8e2.6 times, negative influences of thermal recycling have already been reported. However, the reported issues are not severe enough to endanger a sustainable operation, which indicates the conservative character of the hydraulic breakthrough constraint. This conservative character is also visible in the comparison of the constraints, where the breakthrough limit derives lower flow rates than the drawdown limit. Some of the systems could also be deliberately designed with a hydraulic support to cover the demand and thermal recycling is tolerated to a minor extent. In summary, it can be stated that the technical flow rates have proven their suitability as a conservative peak extraction estimate.

Discussion
The presented method offers a workflow to assess a quantitative potential for the thermal use of groundwater. It integrates the relevant regulatory and operational constraints for a sustainable operation of well doublets with a consideration of their hydraulic impact. Thus, spatial inquiries of technologically achievable flow rates are possible. The necessary functional relations are derived through non-linear regression analyses of numerical simulation results, as reported in section 2.3. The defined hydraulic constraints assume a continuous plant operation, which provides a conservative potential estimation.
The derived approach can be applied in spatial energy planning at the building scale and provides a basis to compare the resource with other renewable options. Further, the assessed potentials can be used for thermal groundwater management, enabling authorities to adapt approval procedures and rise the awareness for this heating and cooling option early in the planning process.
Specific constraints for a sustainable flow rate estimation in open-loop systems have been introduced in previous potential assessment studies (cf. Table 1). Drawdown limits dependent on aquifer thickness are generally used to estimate reliable well yields, with fractions from 25% to 70% of the saturated thickness [33,45]. In this paper, we provide formulae valid for a drawdown threshold of 1 = 3 of the saturated thickness, as prescribed in Munich and recommended in France [31]. However, considering only a linear correlation between flow rate and drawdown does not include all technologically constraining factors. Casasso and Sethi [32] included an additional re-injection limit to avoid flooding with a  safety distance to the surface of 3 m by using the Cooper-Jacob equation [57]. With the presented equation (4), this threshold is adaptable and can be set according to local requirements, as done in Munich with a value of 0.5 m. In addition, the TAP method considers the mutual hydraulic effect of abstraction and reinjection in numerical simulations, which provides specific adaptations for well doublets.
Further, the flow rate of a well doublet is constraint to prevent a hydraulic breakthrough. Only the breakthrough constraint allows a utilization of the potential values for spatial energy planning, because it includes the internal well distance as a spatial influence on the hydraulically affected zone. Since the propagation of thermal anomalies is retarded by a factor of 2e3, the hydraulic breakthrough occurs before the thermal breakthrough [58]. Therefore, it serves as a feasible and conservative threshold to eliminate the risk of thermal recycling [50,59,60]. For the calculation of a hydraulic breakthrough, authors commonly refer to an analytical solution by Lippmann & Tsang [48] and Clyde & Madabhushi [49], which is derived from the same simplifying assumptions and leads to comparable results. M€ oderl et al. derived well spacing formulae for doublets in horizontal and vertical direction [61]. The study also fitted results of parameter variations in numerical simulations to an empirical equation by solving a not further documented minimisation problem. However, the aquifer thickness was not treated as independent parameter and a combination of 5% injected water which return to the abstraction well as inter-flow and 0.5 K of thermal recycling was chosen as acceptable breakthrough criteria for flow rates up to 1 L/s. Although the results are not directly comparable, the well distances by M€ oderl et al. are significantly lower then the results from equation (2), especially for higher flow rates, which confirms the conservativeness of the approach.
Analytical solutions for an optimised external well spacing perpendicular to the groundwater flow are proposed by Javandel & Tsang for groundwater treatment [35]. The derived formulae calculate well spacings that prevent flow to pass between the pumping wells, neglecting an interference between injection wells. Javandel & Tsang did not integrate an influence of the wells on the drawdown. Also the study of Clyde & Madabhushi gives only qualitative information on an external well spacing [49] The objective of our study was to derive the footprint area of hydraulically sustainable well doublets with a balanced water budget. Therefore, the use of numerical models provided a specific solution, which allows spatial queries.

Assumptions and limitations
The presented method assumes several commonly used simplifications [32,33,37]. Since idealised 2D box-models are used for the simulations, the methodology is only suitable for porous aquifers with a certain hydrogeological homogeneity and the method is designed for unconfined conditions, which are more frequent in shallow aquifers. Further, the TAP-method should only be applied within the varied parameter ranges, which however are very wide as shown in Table 3, and in the absence of complex boundary conditions, e.g. recharge from surface water bodies. In the procedure, the hydrogeological parameters are averaged in the assessed area, i.e. the hydraulic footprint or the plot extent. Thus, the error introduced through averaging may increase for larger well distances in heterogeneous conditions. In addition, all flow rate relations are studied in steady-state simulations with fully penetrating wells. This is intended to gain results with general validity, because the estimates represent a constant operation of the openloops system. However, domestic heating and cooling loads, as well as injection temperatures are typically very variable. In consequence, the potential assessment serves as a conservative initial measure, but can not substitute demand specific designs for abstraction or injection, like the use of horizontal filter wells or infiltration ditches. This is also the case for installations where the injection well can not be built directly down-gradient or the extraction occurs at a different elevation from the injection.
In a future perspective, a consideration of thermal anomalies can lead to a spatial assessment of thermal potential, without negative interaction of installations. Additionally, an inclusion of dynamic hydraulic and thermal processes would reveal synergetic effects between installations or enable users to test the influence of different demand scenarios. At a certain point however, a method designed for application in a GIS-workflow can not cope with a sitespecific numerical model, which comprehensively considers the effect of significant influences in transient simulations.

Conclusion
This paper proposes a method to assess the spatial extraction potential for thermal use of groundwater within the relevant operational and regulatory limits. The technical flow rate potential between extraction and injection well is estimated by three hydraulic constraints. The constraints ensure that the drawdown in the extraction well does not exceed a specific limit, that the groundwater table in the injection well does not rise above a specific limit, and that the flow rate does not lead to a hydraulic breakthrough in the well doublet. Additionally, the mutual hydraulic interference of well doublets is calculated to derive a feasible spatial density without inducing a significant cycling of water between the wells.
Finally, the potential assessment was performed in Munich. Since areas with a low aquifer thickness or depth to water are present, the technical flow rate assessment displays the importance of each constraint. A comparison of technical flow rate estimates with monthly groundwater extractions from large open-loop systems showed the suitability of the estimates as conservative peak extraction values. Specifically, the TAP-method delivers quantitative potential values, which are based on a highly transferable and adaptable numerical approach. The method is thus well-suited to serve as a tool for the integration of thermal groundwater use in future energy strategies not only in Munich but in similar environments worldwide.