Model for dimensioning borehole heat exchanger applied to mixed-integer-linear-problem (MILP) energy system optimization

This paper introduces three novel approaches to size geothermal energy piles in a MILP, offering fresh perspectives and potential solutions. The research overlooks MILP models that incorporate the sizing of a geothermal borefield. Therefore, this paper presents a new model utilizing a g-function model to regulate the power limits. Geothermal energy is an essential renewable source, particularly for heating and cooling. Complex energy systems, with their diverse sources of heating and cooling and intricate interactions, are crucial for a climate-neutral energy system. This work significantly contributes to the integration of geothermal energy as a vital energy source into the modelling of such complex systems. Borehole heat exchangers help generate heat in low-temperature energy systems. However, optimizing these exchangers using mixed-integer-linear programming (MILP), which only allows for linear equations, is complex. The current research only uses R-C, reservoir, or g-function models for pre-sized borefields. As a result, borehole heat exchangers are often represented by linear factors such as 50 W/m for extraction or injection limits. A breakthrough in the accuracy of bore-hole heat exchanger sizing has been achieved with the development of a new model, which has been rigorously compared to two simpler models. The geothermal system was configured for three energy systems with varying ground and bore field parameters. The results were then compared with existing geothermal system tools. The new model provides more accurate depth sizing with an error of less than 5 % compared to simpler models with an error higher than 50 %, although it requires more calculation time. The new model can lead to more accurate borefield sizing in MILP applications to optimize energy systems. This new model is especially beneficial for large-scale projects that are highly dependent on borefield size.


Introduction
Designing an energy system with the aim of a minimal economic or ecological impact is a very complex task.Possible applications span from commercial and residential building design, district heating and cooling systems, and industrial process heating and cooling to agricultural applications.The increasing complexity and possibilities in energy sources, time and space lead to an even more difficult problem.These problems are challenging to solve analytically.Instead, mathematical programs can identify optimal energy system designs (Baños et al. (2011)).These programs can use Mixed-Integer-Linear-Programming (MILP) solvers to determine the optimal solution concerning their assumptions and model accuracy during a reasonable calculation time compared to Non-Linear-Programming solvers (Moretti et al. 2021;Epelle and Gerogiorgis 2020).MILP-equations contain only linear equations and real and integer variables (Floudas and Pardalos 2009).
These MILPs are often used to design energy systems (Wirtz et al. 2021a;Gabrielli et al. 2020;Kümpel et al. 2022;Robineau et al. 2016;Sigurdardottir et al. 2015;Wolisz et al. 2017;Gabrielli et al. 2018;Hoffmann et al. 2021;Kotzur et al. 2021b;van der Heijde et al. 2019;Blanke 2022).Geothermal systems and their integration into a multi-source system are essential, especially for districts with a low-temperature district heating and cooling network (Wirtz et al. 2021b;Buffa et al. 2019;Wirtz et al. 2020Wirtz et al. , 2023;;Volkova et al. 2022).Simple linear models assume a limiting power factor like 50 W/m multiplied by the length (VDI 2019).These simple models result in notable inaccuracies, mainly when dealing with large-scale systems.Using R-C, reservoir, or g-function models is more accurate but requires a longer calculation time (Cimmino and Cook 2022;Gabrielli et al. 2020;Kümpel et al. 2022;Claesson and Johansson 1980;Zhang et al. 2019;Hellström 1991;Sigurdardottir et al. 2015).Gabrielli et al. (2020) use a simple g-function model to approximate the temperatures of a borefield in their optimization.Their model is only applicable to optimize operation and not to size borefields.Kümpel et al. (2022) compare a simple g-function and an R-C model.Their results show that the g-function model leads to smaller deviations from the measured data.Both models only apply to optimize operation and not to size borefields.Sigurdardottir et al. (2015) use a reservoir model to optimize operation.Since all these papers only focus on optimizing operations, there is a lack of complex models used to size borefields.
This paper aims to close this gap by providing three methods.Besides, it aims to enable MILP-Energy system optimization to size borefields.These methods are a simple linear model, a mean load model, and a complex g-function model.The simple model limits the power to a limiting factor like 50 W/m to the total borefield depth (VDI 2019).This model is chosen because of its simplicity in calculation and implementation.The mean load model extends the simple model by a factor limiting the average load over a decent number of time steps.This model has been selected since it is a straightforward extension of the simple model and is still easy to calculate and implement.The complex g-function model is based on the approach of Peere et al. (2021) and relies on a g-function to approximate the thermal response of the borefield.This model was chosen since it may yield the most accurate results but is the most complex.An R-C model has not been implemented since Kümpel et al. (2022) favoured a g-function model over it, and it is challenging to parameterize.A reservoir model has not been chosen because it is difficult to parameterize.Furthermore, it is hard to linearize it for sizing in an MILP.The aim is to check which models lead to the most accurate results within a reasonable calculation time.So, they will be compared in terms of accuracy and calculation time.EED and GHEtool are used to compare the MILP results to their accuracy (BLOCON 2023;Peere and Blanke 2022).

Model explanation
The modelling of borehole heat exchangers is complex (Zhang et al. (2019)), and many different models exist (Stober and Bucher (2021); Cimmino and Cook (2022); Blank et al. (2021); Caulk et al. (2016); Claesson and Johansson (1980); Hellström (1991); Lund and Östman (1985)).The considered models are a simple model, a mean load model and a complex model.The simple model limits the extractable or injectable power to a fixed value multiplied by the borehole length.The mean load model extends the simple model by a limiting power factor for a running average load over a certain number of time steps.The complex model uses a g-function and monthly loads to limit the power.

Simple model
The simple model is limiting the maximal extracted ( P + g ) or injected heat ( P − g ) to a fixed value ( f + g , f − g , i.e.50 W/m) times the borehole depth ( L ) and number of boreholes ( n): Furthermore, the length can either be used for extracting or injecting heat.The following equation ensures that the geothermal system is not simultaneously used over its capacity: The simple model has the fewest and most simple constraints.Therefore, it does not consider any degradation of the borefield over time due to the previously injected or extracted heat.The model will likely calculate fast but with a significant error.The advantage is the easy implementation and a potentially short calculation time.However, its disadvantage is that it may cause short-and long-term inaccuracies.Its scope of application lies in small-scale projects, whose results do not depend heavily on the borefield size or simple manual calculations.

Mean load model
For the mean load model the simple model equations are extended to ensure, that the peak load of heat, which will be extracted from or injected to the underground, is below a certain limit for an average over a certain number of time steps l ( F + g , F − g ): (1) (2) (3) The mean load model has a higher number of and more complex equations than the simple model.It considers short-term degradation of the borefield over time but misses long-term effects.The model will most likely calculate slower than the simple model, potentially leading to a lower error.The advantage is the consideration of short-term degradation, easy implementation, and potentially short calculation time.However, its disadvantage is that it may cause an inaccuracy in the long-term effects.

Complex g-function model
The g-function model is considering a g-function for a fixed depth and borefield configuration.It uses the monthly average load to determine the temperature gradient over the years.
The monthly average loads of month m use an equally yearly split of 730 h as the g-func- tion is calculated for this time step: The g-function for the peak ( g p ) and the monthly load ( g * ) is calculated using GHEtool (Peere and Blanke (2022); Peere et al. (2021)).
For all time steps t the temperature constraints for the extracted and injected heat have to be met based on Peere et al. (2021): These equations can be simplified with maximal and minimal temperature differences to the ground temperature ( T max , T min ) to: (4) (9) The depth of the borehole ( L ) is a variable and can therefore be extracted from the mid- dle part to the boundaries: This constraint has to hold for every time step t within the first and last year of the geo- thermal system.
To simplify the problem, the convolution term has been separated ( P M,con ) and a monthly maximal ( P max ) and minimal value ( P min ) are implemented: However, there exists a challenge associated with the given equation, namely, the temperature difference tends to be dependent on depth.Consequently, the product of T and L leads to a quadratic term.This paper will elucidate two potential solutions to address this issue.The initial solution involves computing an average value over the (11) (15) P M,con = P M * g * , (16) ) ∀t anticipated depth range.For instance, when considering a depth range from 50 m to 150 m, the average ground temperature can be determined.In this scenario, the average depth ( L ) is set at 100 m.The second option involves introducing a piecewise linear term to approximate this behaviour.
The first option can be implemented using a geothermal heat flux ( α [ W/m 2 ]) and a start temperature difference ( T 0 ): As an alternative, the temperature differences can be determined using a piecewise linear approach by introducing new variables ( T L,i ): While they are linearly ( m i • x + b i ) dependent on the bore hole depth, binary ( B i ) and trash variables ( T L,tr,i ) are needed.The trash variables store the values if the bore field length exceeds their range.m i is the current piece's gradient over borehole depth, and b i is the axis intercept: The further equations to ensure that just one new variable is selected are explained in the Appendix.
The complex model has the highest number of variables and most complex equations.It considers short-term and long-term degradation of the borefield over time.The model will most likely calculate slower than the other models, leading to the lowest error.The piecewise linear temperature difference should lead to a higher calculation time and more accurate results than the constant temperature difference.This model has the advantage of potentially high accuracy for short-and long-term effects.However, it also has the disadvantages of potentially high computational time and complex integration into MILPs and practical applications.Its scope of application lies in large-scale projects whose results are highly dependent on the borefield size.

Methodology
Three different cases will be evaluated to determine the influence of the different models on real borefield dimensioning tools, such as EED or GHEtool.Figure 1 shows the energy systems for these cases.The first system, called "Case New" (representing a new building, with actual energetic standards), has to cover the heat and cooling load of a residential building.Therefore, it can use a heat storage, a geothermal heat pump, and an electrical heater to meet the heating demand.The cooling demand has to be covered by the geothermal system, the cold storage and the geothermal heat pump.The second system is called "Case Old" (representing an old, existing building).It just has to cover the heating demand of the same kind of building.So, it does not use cooling.The third system, called "Case District" has to cover the heating and cooling demand of a district (representing a district with a couple of (19 new buildings).Therefore, it is using the same technologies as Case New.Besides, it can use an electrical chiller to cover the cooling demand.These three cases are the most common applications of borehole probes in building contexts.Case Old represents a refurbishing case of an existing building.Case New is a newly constructed building that can also regenerate the probe.Since the 5th generation district heating grid is a typical new heating grid, the Case District is considered here.
Figure 2 shows the schematic of this study for the different cases considered.The following steps are applied to every borefield configuration.First, the energy system will be optimized considering the case's loads and options and the borefield configuration.The energy system results are reformatted to be usable in EED and GHEtool.Afterwards, the depth considering the injected and extracted heat is determined using EED and GHEtool.The depth and temperature error is then calculated.This process is repeated for all borefield configurations.A detailed explanation of each step will be provided in the upcoming section.
A MILP problem is formulated with the primary goal of efficiently minimizing the investment and operative costs, as described below.It uses the abovementioned geothermal model constraints and the cost constraints stated below.Furthermore, the energy model equations are explained in the Appendix since they are not the focus of this study.The model has been solved using GUROBI 10.0.0 with the maximal presolve value and a heuristic value of 0.01 (Gurobi (2023)).The rest of the settings are kept as default values.
The MILP model equations are based on Blanke et al. (2022); Welder et al. (2018) and the prices are based on Kotzur et al. (2018b).The equations for the energy system components are explained in section 4. The price and technical assumptions are explicitly not the focus of this study.The aim of the optimization is to minimize investment ( C investment ) and operation costs ( C operation ): The investment costs are a combination of the geothermal heat pump ( C HP ), electrical heater ( C EH ), heat storage ( C HT ), cold storage ( C CT ), geothermal system ( C Geo ) and the electrical chiller ( C EC ): (22) min(C investment + C operation ).The operation costs are the summation over the 8760 h of the year and the bought electricity ( P buy ) using a tax-free electricity price of 0.25 €/(kW h) (BDEW 2019): The analysis encompasses a wide range of parameters, including varying thermal conductivities of the underground (ranging from 1 to 4 W/(m K) in a 1 W/(m K) step), different borehole spacings (3, 5, 7, 9 m for Case New and Case Old, and 4, 6, 8, 10 m for Case District), and various borefield configurations.For Case New and Case Old, the borefield configurations include 1x1 pipe, 2x1 pipe, 3x1 pipes, and 2x2 pipes setups.For Case District, the configurations consist of 25x25 pipes, 60x12 pipes, 39x13 pipes, and 30x30 pipes.Case District has different spacing due to the configurations with an high amount of borehole pipes.This results in a total of 80 sensitivity cases.
The thermal capacity is set at 2160 kJ/(m 3 K) , while the borehole's thermal resist- ance is 0.05 m K/W.This thermal resistance value corresponds to a coaxial steel pile with an outer steel pile (diameter (outer/inner): 150 mm/ 140 mm) and an inner PE pipe (diameter (outer/inner): 90 mm/84 mm) with a mass flow rate of 2 kg/s (Blanke   (2021b, 2021a)).The heat transfer fluid used is a 30% glycol-water mixture, and the burial depth is 5 m.
All load profiles are based on the test reference year for Düren TRY 2045 (DWD (2017)) and are provided in Blanke (2023).The yearly heating load for the Case New adds up to 8172 kWh and the cooling load adds up to 874 kWh.The yearly heating load for the Case Oldadds up to 18 040 kWh.The yearly heating load for the Case District adds up to 3056 MWh and the cooling load adds up to 475 MWh.For Case Old and Case New, the heating and cooling demand is multiplied by the number of boreholes.In these cases, scaling is imperative as the borehole configurations result in insufficient depth for the models applied.In the district case, this can be compensated by the electrical chiller.
For the MILP, the following assumptions regarding the geothermal system have been made.For the simple model, the gain is calculated using equation 25.The g-function is determined with GHEtool using the bore field configuration and ground properties: For the mean load model, the same approach is used to determine F + g and F − g but with a peak duration of six hours instead of one.For the complex model, the g-function is provided also from GHEtool using a start depth ( L start ).The temperature difference with a piecewise linear temperature difference (Complex G) and a constant temperature difference (Complex N) are considered for the complex model.The fluid temperature needs to stay 10 K below and 7 K above the ground temperature.The start depth results for all cases from GHEtool sizing of the energy system's heating and cooling load fully covered by the heat pump and the geothermal system.The calculation with EED were done using the same input data (ground properties, borehole and geothermal probe, heat carrier fluid, heat/cold loads) as for the GHEtool.EED does not include the option to set a certain value for the buried depth..
Besides the full year optimization, a typical period's optimization is performed.It is done for 10, 20, 30, 40, 50 and 60 typical periods with 24 h's length.The clustering is done by the Time Series Aggregation Module (TSAM) using a k-mediods clustering algorithm (Kotzur et al. 2021a(Kotzur et al. , 2018a;;Hoffmann et al. 2021).Further details on the clustering algorithms can be found in Kotzur et al. (2018a).
Besides the MILP-optimization, the results determined by the MILP are compared with results from GHEtool and EED.Three different scenarios will be compared.The temperature error with a ground temperature increase of 3 K/100 m and the depth to cover the load profile with temperature increase with depth will be calculated.The results will be compared with the depth and temperatures calculated by the MILP.Furthermore, the calculation time needed to solve the problem on an Intel Xeon W-2155 CPU will be compared. (25)

Results
Tables 1, 2 and 3 in the appendix list the temperature and depth results for the different configurations and models.The following figures illustrates this results.There are no tables for the typical periods results since there are too many results.
Figure 3 shows the temperature error between the MILP model and the GHEtool calculation for the MILP-determined depth for the different cases and models.The results are shown as a box-plot diagram with the corresponding average value ( + ).A positive error means temperature exceeds the limits of 0 • C or 17 • C and a negative value means that it stays below the limit.In this case, it is not using the complete potential.The simple model leads to the largest error for all three cases with its maximal error above 125 K.The second worst results for all three cases is achieved by the average model.The complex models achieve the best results with an error below 1 K or negative values for the first two cases.For the Case District, the error is the highest for all models.The complex models achieve the lowest error with a maximal error below 4 K.The complex model with no gradient (Complex N) achieves the lower errors than the complex model with the piecewise linear temperature gradient (Complex G).The complex model error increases with the difference of the depth for which the g-function is calculated.This trend is especially dominant in the Case District.Especially for Case District, the error is higher because the MILP-sizing is reaching the depth boundaries of 50 m or 150 m in more sensitivity cases.Figure 4 shows the depth error.A value below zero means that the MILP-depth is higher than the one calculated by the software.The simple model determines a depth which is too small to cover the demand.The maximal depth difference laying outside the graph is − 690 m which is the case for the simple model.The average model leads to an error in nearly the same range as the simple model, while the median and average have an lower error than the simple model.The complex models lead to the smallest error.The error is especially for the first two cases negligibly small.The complex model with no gradient (Complex N) achieves lower errors than the complex model with the piecewise linear temperature gradient (Complex G).
Figure 5 shows the calculation time for the MILP and the different models.The simple and mean load models are computed for all cases in under 13 s.The complex model, on the other hand, exhibits notably longer calculation times, sometimes exceeding 400 s at its maximum.This time increase depends on the case.In the simple Case Old, the calculation time is around 150 s and around 70 s to 150 s for Case New.In the case of Case District, this time increases to approximately 260 s, and it also shows a wider distribution across the sensitivity cases.The complex model with no gradient (Complex N) achieves a better calculation speed than the complex model with the piecewise linear temperature gradient (Complex G).The exception is the Case New.Here, the calculation time of the Complex G is around the same as for Complex N.
Figure 6 shows the depth error if typical periods are used.The upper plots (a-c) represent the depth error for a full year simulation and the lower plots represent the depth error compared to GHEtool sizing.The error compared to the full year simulation is low especially for the first two cases and decreases over the number of typical periods.The most significant error occurs in the Case District, with depth inaccuracies below 100 ms for a number of typical periods of 20 and beyond.The error compared to the GHEtool sizing is nearly the same as in the full year simulation.The calculation time has been reduced using typical periods by to factor of 2 to 10, depending on the number of typical periods.
Figure 7 shows the depth error between the optimized depth to GHEtool and EED.The figure shows these results for different thermal conductivities and the three considered cases.The GHEtool error decreases with increasing thermal conductivity with the exception of Case District.The EED results are close to the GHEtool results.The EED and GHEtool error is for the first two cases within 7 m.For the Case District, this depth error is in most cases below 20 m.

Discussion
The results show that the presented complex models lead to the lowest errors.This low error is reasonable since the simple and mean load model does not consider long-term effects.The error is small for cases close to the g-function depth.This minor error is reasonable given the relatively small error of the g-function in that region.This trend, however, is a bigger problem for large bore field geometries because a slight depth difference in borehole depth leads to a big difference in the total length.However, this error is still significantly lower than for the simpler models.The complex model with a piecewise linear temperature gradient leads to more minor errors than the one using a constant temperature gradient.The gradient error and the g-function error compensate each other to some extent.The complex models come with the trade-off of a significant calculation time increase.
The application of the method using typical periods can solve this issue.The usage of typical periods reduces the calculation time significantly.For more than 10 typical periods, the calculation error is low.Typical periods are a suitable option to reduce the calculation time.
The results were cross-verified with EED in addition to GHEtool, yielding nearly identical outcomes.This similarity confirms the model's applicability for depth sizing.Discrepancies between GHEtool and EED primarily stem from varying program assumptions.
One key factor contributing to these differences is the option to set the buried depth in GHEtool, which is not available in EED and can influence the obtained results.
This new complex model has multifaceted applications spanning commercial and residential building designs, district heating and cooling systems, industrial processes for heating and cooling, and agricultural applications.Architects and engineers can optimize energy-efficient structures in commercial and residential building design by sizing geothermal borefields according to heating and cooling demands.Municipalities and district heating companies can optimize their systems by sizing geothermal bore fields appropriately, leading to more efficient district heating and cooling networks.Industrial facilities can optimize energy usage and reduce environmental impact using accurately sized geothermal systems for process heating and cooling.Similarly, agricultural operations can benefit from optimized geothermal installations for climate control, ensuring ideal growing conditions while minimizing energy costs.

Conclusion
This paper presented three approaches to size geothermal energy piles in a MILP.The research overlooks MILP models that incorporate the sizing of a geothermal borefield.Therefore, this paper presents a new model utilizing a g-function model to regulate the power limits.
The MILP models are compared for three energy systems, different bore field configurations and ground properties.These results show that the new g-function-based model leads to significantly lower errors than simpler models and can approximate the geothermal behaviour.Besides, the simple models lead to significant temperature and depth sizing errors far above 100 %.This high error emphasizes the importance of this research and the new model.However, this comes with the trade-off of an increase in calculation time.Also, the model is capable of handling typical periods.Application of typical periods can reduce the calculation time without increasing the error.
Therefore, adopting the new g-function model is strongly advisable, particularly in scenarios where the outcome is exceptionally sensitive to the geothermal energy system's size.This sensitivity can be the case for the energy system optimization of commercial or residential buildings.Furthermore, the model presented here can effectively restrict geothermal power when sizing is unnecessary.Consequently, the complex model is also recommended for optimization related to operational efficiency.This operational efficiency can be essential for optimizing the operation of district heating or cooling grids.
Future research can focus on further validation with a broader and more diverse dataset and other Software tools to strengthen the findings.Furthermore, the model can be performance optimized to reduce the calculation time by reformulating or simplifying the constraints.Also, approaches are possible to calculate the number of boreholes with their binary selectable g-function instead of sizing the length.This approach would lead to more binary variables but also more accurate results.Future studies could explore further refinement of the model.The model's applicability in different geographical contexts can be explored.Further economic analyses can enhance the practical value and impact of the model.

Appendix
See Tables 1, 2 and 3.         Using the maximal cooling power of the load profile ( max( Qcool ) ) the reference electri- cal chiller power is limited: The electrical power of the electrical chiller ( P EC,el ) is limited to the reference power divided by the energy efficiency ratio (EER) of 2.82 based on a reversible heat pump of the hplib (W5A35) (Hoops et al. 2022): The cooling power ( QEC,cool ) is the electrical power times the COP:

Heat storage
The heatings storage costs are linearly dependent on the reference size ( Q HT,ref ): The storage state ( Q HT ) is limited to the reference size: The storage state is dependent on the charging ( QHT,ch ) and discharging power ( QHT,dis ) with its corresponding efficiencies taken from Kotzur et al. (2018b): The storage state at the start ( t = 1 ) and end ( t = 8760 ) are linked: The typical period's optimization uses the model of Blanke et al. (2022).

Cold storage
The cold storage costs are linearly dependent on the reference size ( Q CT,ref ): The storage state ( Q CT ) is limited to the reference size:    (43) The storage state is dependent on the charging ( QCT,ch ) and discharging power ( QCT,dis ) with its corresponding efficiencies taken from Kotzur et al. (2018b): The storage state at the start ( t = 1 ) and end ( t = 8760 ) are linked: The typical period's optimization uses the model of Blanke et al. (2022).

Geothermal system
The geothermal system costs are linearly dependent on the reference size ( L ) and on the binary variable for the constant price ( b Geo ): Using the maximal depth ( L max ) the length is limited.The maximal depth is the start depth multiplied with the width and length of the geothermal system times ten: The rest of the equations for the geothermal system are explained above.

Balances
The summation of all heat fluxes of the components has to cover the heat demand ( Qheat ): The summation of all cooling fluxes of the components has to cover the cooling demand ( Qcool ): The bought electricity is the sum of electricity for the geothermal heat pump, electrical heater and electrical chiller:

Piecewise linear geothermal equations
This section explains the further equations to create a piecewise linear temperature gradient.
The new variables and trash variables are linked to the binary variables so just one of them can be active:     (50) 0 = P buy (t) − P HP,el (t) − P EH,el (t) − P EC,el (t). (51

Fig. 2
Fig. 2 Schematic of the optimization study

Fig. 3 Fig. 4
Fig. 3 Temperature error between the temperature limits and the predicted temperatures calculated with the g-function for the optimization determined depth for the different models and cases [Case New (a), Case Old (b), Case District (c)]

Fig. 5 Fig. 6
Fig. 5 Calculation time for the different models and cases Case New (a), Case Old (b), Case District (c)

Fig. 7
Fig. 7 Depth error for the different models using EED and GHEtool over different thermal conductivities and cases [Case New (a), Case Old (b), Case District (c)] 1812 C /k W • P EC,ref + 4729 C • b EC .(35) P EC,ref ≤ b EC • max( Qcool ).

Table 1
Results of case New (n.c.= not calculated)

Table 2
Results of case Old (n.c.= not calculated)

Table 3
Results of case District (n.c.= not calculated)