Cloud-Based Optimal Control of Individual Borehole Heat Exchangers in a Geothermal Field

Integrating renewable energy sources is a crucial component in reducing CO2 emissions in the building sector. In particular, shallow geothermal energy is expected to play a significant role in the regenerative energy supply of buildings. An effective control strategy for the geothermal field is crucial to reduce the overall energy consumption. This paper analyzes the benefits of controlling an existing field’s individual borehole heat exchangers (BHE) using nonlinear model predictive control (NMPC) and moving horizon estimation. The considered geothermal field consists of 41 BHEs and is used for heating and cooling. Each BHE is equipped with temperature sensors for in- and outflow and has individually controllable valves, while a central hydraulic pump feeds all BHEs. The sensor measurements are accessed through a cloud platform, enabling also set point writing for the pump speed and the valve positions. To control the BHEs individually, we propose a two-stage process. In the calibration stage, a moving horizon estimator estimates the actual borehole and ground temperatures for each BHE. In the second stage, first, a nonlinear model predictive controller optimizes the number of active BHEs necessary to meet the buildings’ energy demand. With the estimated ground temperatures as a basis, it is determined which BHEs shall be (de)-activated. The active BHEs are fed with a fixed volume flow of 24 L/min to ensure turbulent heat transfer. To reduce the power usage of the pumps, an optimal control problem based on a simple hydraulic model of the geothermal field is used. The methodology is analyzed through simulations first and then validated experimentally. The results show that half or more of the BHEs could be deactivated most of the time, leading to 67% savings in electricity consumption by the hydraulic pump. The experimental validation confirms the high energy saving potential of the proposed methodology, reducing the consumption of electrical energy by 71%. Additionally, the deactivated BHEs regenerate faster and improve the field’s long-term behavior. In conclusion, the proposed strategy improves the short and long-term performance of the geothermal field.


Introduction
Integrating renewable energy sources in buildings becomes more and more critical to tackle climate change and the depletion of fossil energy resources. At this point, shallow geothermal fields provide a stable and efficient heat source for heat pumps, which can also be used for cooling applications. For these reasons, the use of shallow geothermal energy increased over the last years [1]. The geothermal energy is harnessed in geothermal fields by pumping heat transfer fluid through borehole heat exchangers (BHE). On the one hand, the field's efficiency is primarily influenced by the power consumption of the hydraulic pumps [2].
On the other hand, the temperature of the ground influences the heat pump's coefficient of performance (COP) and, therefore, the power consumption of the heat pump. If the field is imbalanced, for example, primarily used for heating, the ground temperature changes, decreasing the long-term efficiency of the field and the heat pump [3][4][5]. A well-tuned operation strategy is necessary to ensure an efficient operation of the field in the short-term and sustainable use in the long term.
Due to its predictive behavior and its ability to deal with multiple-input multiple-output (MIMO) systems as well as constraints, model predictive control (MPC) is well suited to be such a strategy. MPC for building control application gained increased interest in recent years and is also applied to building energy systems powered by geothermal energy [6][7][8][9].
In this work, we present a model predictive controller for a real-life geothermal field equipped with extensive monitoring, which offers the unique capability to control the single boreholes of the geothermal field individually. To evaluate the controller's performance, the methodology is tested on a simulation model of the geothermal field first. Afterward, we show how a state-of-the-art cloud platform is used to implement the algorithm in the real field and present the first results of the real-life demonstration.
In the following, we discuss related work on modeling and model predictive control of geothermal fields and present our real-life example. Afterward, the methodology is explained, and its performance is evaluated.

Related work
Developing a model predictive controller requires an accurate system model to control. The challenge of modeling geothermal fields is to account for short-and long-term effects like borehole interactions. A comprehensive review of modeling geothermal fields for control applications is given in Ref. [10]. Thermal response factor-based g-function models often model the long-term effects of a geothermal field. These models employ a temporal and spatial superimposition of the temperature distribution around each BHE induced by the heat extraction rates in discrete periods. A detailed description of g-function models and a toolbox to model arbitrarily positioned and sized boreholes are presented in Ref. [11]. To account for the short-term dynamics of boreholes, resistance-capacitance (RC) network modeling approaches are often used [12][13][14]. Furthermore, a combination of both techniques can be used for detailed modeling.
Several examples of model predictive controller implementations for geothermal fields can be found in the literature. In Ref. [4], a model predictive control strategy based on a g-function model for an undersized borehole in combination with an electric heating element is deployed, reducing the peak energy demand of the borehole by 58%. Weeratunge et al. [15] present a model predictive controller for a solar-assisted ground source heat pump system based on a g-function model combined with mixed-integer optimization. Verhelst [7] developed an MPC scheme to control a ground-coupled heat pump in combination with concrete core activation. The author concludes that an RC model for the geothermal field with less than six states shows sufficient control quality. Optimizing the geothermal field's load together with the heat pump operation results in energy savings of 20% to 40% compared to a baseline strategy for different design cases. Picard [8] optimizes the volume flow rate of a geothermal field, assuming a constant wall temperature and linear heat pump characteristics. In Ref. [9], a shadow cost function combined with active regeneration is introduced to assure the long-term efficiency of a geothermal field. Compared with a baseline MPC controller, this approach reduces the total operating costs over ten years by 9.7%. The authors of Ref. [16] implement a model predictive controller in a real-life building powered by a ground source heat pump. Compared to a standard rule-based controller, the MPC saves 53.5% energy consumption while increasing comfort by 36.9% in transient seasons. Atam et al. [17] compare a model predictive controller to dynamic programming for a ground source heat pump system. In this work, the authors use a reduced order model of a borehole with six states in a nonlinear controller to predict the temperature dynamics accurately.
To initialize an MPC, knowledge of the current systems states is crucial. An estimator is required since not all states can be measured, like the underground temperatures in an RC model. In Ref. [13], the performance of a Moving Horizon Estimator (MHE) and a Time-Varying Kalman Filter (TVKF) to estimate the load history of a geothermal field is evaluated. It is shown that the MHE performs better than the TVKF with the drawback of higher computational costs.
In the literature presented previously, the geothermal field operation is either manipulated by optimizing the overall mass flow or the overall heat extraction and insertion, respectively. An approach to optimize the load of individual boreholes using linear programming is presented in Ref. [18]. With the implemented strategy, the temperature change after 30 years in the ground can be decreased by 18%, advocating the usage of individually controllable boreholes. This methodology is extended by Hecht-Méndez et al. [19] to include groundwater flows. Based on this work, Bayer et al. show that the number of boreholes can significantly be reduced if the heat extraction of individual boreholes is optimized [20,21]. The latter can significantly reduce the construction costs of geothermal fields.
A simplified approach is proposed by Luo et al. [2], where a geothermal field consisting of 596 boreholes is grouped into twenty units. Based on the timely energy demand of the building, the number of active groups is assigned. This strategy leads to savings of 3000 EUR in operation costs.
An example of a real-life borehole field with individual controllable field sections is the borehole field in Emmaboda, Sweden [22]. However, the focus of this work is the utilization of waste heat, using rule-based control strategies for sections of a geothermal field and not individual boreholes [23].
Therefore, this work's novelty is a methodology for how a borehole field's short and long-term performance can be improved by individually controlling single boreholes with a model predictive controller and moving horizon estimation. Furthermore, the model predictive controllers for geothermal fields presented in the literature are solely based on simulation studies. In contrast, we present a real-life demonstration of a model predictive controlled geothermal field.

Case study-the main building of the E.ON Energy Research Center
We develop the methodology for the geothermal field of the E.ON Energy Research Center's (ERC) main building in Aachen, Germany. The building consists of office rooms and laboratories on 7222 m 2 with a complex heat and cooling demand. The energy system comprises different heat and cold suppliers to meet this demand. The central heat and cold supply component is a reversible heat pump with a nominal thermal power of 250 kW. In addition to that, a combined heat and power plant and two gas boilers are used to cover peak loads in heat demand. More than 12 000 data points provided by different sensors in the building are monitored. The building itself is used as a demonstrator in several research projects. A detailed description of the energy system and its current mode-based control strategy can be found in Refs. [24] and [25]. The methodology developed in this paper is applied to the ERC geothermal field with 41 BHEs, which is shown in Fig. 1. The BHEs are organized in three shafts (east, west, and south). The southern shaft is positioned under a meadow. The eastern and western shafts are predominantly located under an asphalt surface. The field is used for either free cooling or heating in combination with the heat pump. The current mode-based control strategy of the energy system uses the return temperature of the field to decide whether to use it. When using the field, a constant volume flow of 30 L/min is applied to each BHE. For this purpose, two pumps (Wilo ChronoLine-IL-E 65) with a maximum electric power of 6 kW are running alternating to prevent wear. In 2020, the pumps had a power consumption of 35.39 MWh.
The 41 BHEs have individually controllable valves (Belimo R413 control valve with electric actuation) and are each equipped with temperature sensors for in-and outlet flow and ultrasonic volume flow sensors. Since 2014, the sensor measurements have been stored in a cloud platform with a sampling rate of about 30 s. The cloud platform can also send control signals to the 41 valves and the two pumps.
The available monitoring data reveal several weaknesses of the implemented energy management. Since the field is unbalanced and mainly used for cooling, a temperature increase in several parts of the field was detected. Especially the eastern shaft shows poorer cooling performance due to denser BHE placement compared to the other shafts. On the other hand, the field's significant, constant volume flow results in a slight temperature difference between in-and outlet during part-load operation. Thus, deactivating selected BHEs while keeping the nominal mass flow in the active ones could lead to a more efficient field operation by increasing the temperature difference and reducing the necessary pumping power. Additionally, the deactivated BHEs can regenerate and thus improve long-term performance.
In the following, a methodology to determine the optimal number of active BHEs and select the activated ones is presented and tested in simulation and experiment.

Methodology
We developed a control scheme based on nonlinear model predictive control, moving horizon estimation, and optimal control to improve the field efficiency. In this section, we first present the overall control strategy before giving a detailed description of the individual components. Finally, we explain the cloud-based implementation of the proposed method.

Proposed control scheme
An overview of the proposed methodology is given in Fig. 2. Starting with measurements from the geothermal field, an MHE is used to estimate the initial temperatures T model,i,o , and parameters of an RC model for each BHE i. With new measurements from the field, the RC models then continuously simulate the ground temperatures T m o d e l , i . Based on the mean simulated ground temperatures of the active BHEs, a nonlinear model predictive control (NMPC) scheme is initialized. The NMPC uses a heat exchanger model in combination with a model of the geothermal field. Here, the geothermal field is represented by a single BHE RC model with the mean temperatures T model,mean scaled by the number of active BHEs n BHE,active . The heat exchanger connects the geothermal field's glycol circuit to the building's water circuit. The NMPC then optimizes the number of active BHEs n BHE,active to meet the thermal demand of the building at the heat exchanger. In the next step, the algorithm decides which BHEs to use. At this point, the BHEs that need recalibration due to a large modeling error (evaluated as the root-mean-squared error (RMSE model,i )) are activated first. A recalibration with the MHE is requested when the modeled outlet temperature of active BHEs deviates too much from the measured outlet temperature or when a BHE has been inactive for more than 30 days. The other BHEs are sorted according to their ground temperature depending on the load case. The n coldest BHEs are selected for cooling, and the n Fig. 2 Overview of the developed control algorithm warmest BHEs for heating. In the final step, an optimal control problem (OCP) is solved to determine the pump speed N pump and the valve positions u valve,i to obtain a nominal volume flow of 24 L/min in the activated BHEs while minimizing pumping power. The nominal volume flow of 24 L/min ensures efficient, turbulent heat transfer in the boreholes. In the following sections, a detailed description of the individual components of the presented strategy is given, starting with the RC borehole model.

Short-term modeling of geothermal fields
In this section, we present and validate the RC BHE model used for the simulator, the MHE, and in the NMPC. The short-term field dynamics, neglecting the interaction of the BHEs, are captured by a simple RC approach. Here, lumped capacities with temperatures T j in discrete shells j with diameters d j around the borehole are used to model the thermal behavior of the soil, as shown in Fig. 3.
With the in-and outlet temperatures of the heat carrier fluid T in and T out and the thermal power of the probe fluid Q  . In this model, all relevant parameters like the capacities C j , heat transfer coefficients k, and the areas A are derived by material properties, geometry, and standard procedures like a thermal response test. Since this model is developed and used for short-term performance prediction, the temperature T soil adjoining the last shell is kept constant and must be estimated by the moving horizon estimator. Therefore, the remaining model tuning parameters are the number and diameters of the soil shells.
When choosing the shell diameters, the minimum shell size is defined by the radial penetration min min min , 5 , with the thermal diffusivity of the borehole α and the shortest time interval t min for a change in load condition [26]. Assuming the shortest time interval to be 24 h, the minimum shell size is calculated to r min =0.41 m. Thus, the borehole diameter, 0.3 m, is chosen as the minimum shell size. Eskilson [26] recommends the maximum shell size based on the maximum simulation time max , which is the maximum radial penetration depth for the given time. In the given example, this depth is 4.9 m for a maximum simulation length of one month. Since later on, the MHE would need an equal time horizon to provide an educated guess of the initial values, we allow a smaller shell size than the recommended maximum for model calibration.
Since all lumped shell temperatures cannot be measured and need to be estimated by solving an optimization problem, the number of shells should be kept small. Therefore, different model formulations are calibrated and tested on historical data to evaluate the necessary model complexity. For model validation, measurement data of four representative months, including a heating, cooling, and transition phase, were used. Fig. 4 shows exemplarily the temperatures of a three-shell model for the fifteenth BHE. Here, the first day of the month is used to estimate the initial temperatures as well as the adjoining soil temperature as a constant parameter. Apparently, the model is well suited to predict the outlet temperature T out,mod for a given inlet temperature and mass flow over one month. The root-mean-squared error (RMSE) for the outlet temperature prediction is referred to evaluate the model quality. The RMSE of models with two, three, and four shells is presented in Fig. 5 for each week of the predicted month. Overall, it can be seen that all variants achieve very good accuracy but become less accurate with increasing prediction horizon. The latter results from the assumption of a constant adjacent temperature T soil that is not valid with increasing time. It is further shown that the three-shell model is best suited for forecasts up to three weeks and is only slightly worse than the four-shell model with a larger maximum shell size at a horizon of four weeks.
For this reason, the three-shell model is used for further development. The entire geothermal field is simulated in a computationally efficient way by 41 single BHE models that receive the same inlet temperature, and individual volume flows. The thermal power of the field is then the sum of the thermal power of each BHE.
In the following, we show how the MHE uses the presented model.

Moving horizon estimation of ground temperatures
To predict the performance of the individual BHEs with the analytical BHE model, the current ground temperatures of the RC model must be known as initial values as well as the adjacent temperature T soil . These can be determined with a Kalman filter or a Moving Horizon Estimator [13]. An MHE determines the states and parameters of a process model by minimizing the deviation of the model from the measured values for a measured period M. Analogous to model predictive control, this is done with a moving horizon. This means that the last M measured values are taken into account to determine the current variables. The simplified scheme of an MHE is represented by Eqs. (5)-(7) [27].
At this point, y k denotes the modeled outlet temperature andˆk y denotes the measured outlet temperature. Possible unknown measurement noise is considered in the term w k . In order to consider past measured values that lie outside the horizon M, the arrival cost term  [28]. Since the optimization problem is linear with respect to the variables to be optimized, the linear programming optimizer Gurobi [29] is used.
The main design parameter of the estimator is the length of the horizon M and the time step size dt to discretize the differential set of equations. Fig. 6 shows the accuracy of a monthly simulation initialized with estimates from different estimation horizons. However, a short estimation horizon leads to acceptable results that differ strongly from BHE to BHE. Thus, a short estimation horizon leads to unreliable results. On the other hand, a long horizon does not necessarily improve the overall result but significantly reduces accuracy variations. However, since the respective BHE must be active for parameter estimation, a long horizon is not practical. For this reason, 24 h with a time step size of five minutes is chosen as the horizon. Thus, each BHE must be operated at least for one day for parameter estimation.

Nonlinear model predictive control (NMPC) for part-load operation
During the analysis of the historical data of the field operation, it became clear that the geothermal field is only operated with a minimal difference between flow and return temperature. The latter results from the mode-based operation strategy, according to which all BHEs are operated, although the building only requires little power from the geothermal field [25]. A targeted deactivation of individual BHEs with simultaneous control of the active BHEs to a constant nominal mass flow leads to a lower field-side mass flow and, therefore, an increase in the temperature difference in the heat exchanger. At the same time, pumping power can be saved, and the deactivated BHEs can regenerate. In order to determine the number of active BHEs necessary to meet the building's energy demand, a nonlinear model predictive controller was developed, which takes into account the behavior of the field as well as the behavior of the heat exchanger. Fig. 7 shows the system boundaries and the relevant variables. The NMPC uses the RC model presented in section 2.1, scaled by the number of active BHEs n BHE,active . The RC model is initialized with the average temperatures of the active BHEs from the last time step. When selecting a suitable heat exchanger model, the influence of the two mass flows on the transferable heat must be captured with sufficient accuracy. At the same time, the feasibility of the optimization problem has to be ensured. For this reason, the cross-convection model proposed in Ref. [30] was used. The following equations summarize the model: In this approach, the outlet temperatures on each side of the heat exchanger (T out,HX,B for the building side and T in,gtf for the geothermal field side) are related to inlet temperatures (T in,HX,B for the building side and T out,gtf for the geothermal filed side)and mass flows gtf m  and B m  .
Here, the parameters p 1 , p 2 , d, and c were calibrated on measured data. The behavior of the heat exchanger is assumed to be stationary so that the time derivatives with the time constants τ 1 and τ 2 can be neglected. A minimum temperature difference of 0.4 K is introduced for a more conservative estimation, which corresponds well to the measurements. When used in the NMPC, this allows the optimization of the transferred heat flux k Q  and the number of active BHEs. In the NMPCs cost function (12) With the variables: At this point, the design parameters of the control are the prediction N=24 h, horizon, time step size dt=1 h, and the weights Q=5 and R=0.1. A collocation scheme with three collocation points is chosen to discretize the nonlinear differential equations. To avoid integer variables in the optimization, the number of active BHEs is treated as a continuous quantity, which is then rounded up for a conservative approach. The controller is implemented in Pyomo [31] using the interior point solver Ipopt [32].
To maintain a minimum flow rate of the pump, the minimum number of active BHEs is set to five. Based on the optimized number of BHEs, the individual BHEs to be activated have to be identified. For that reason, we first select the BHEs that need recalibration, either because they have been inactive for more than 30 days or because the model error becomes too large. The other BHEs are sorted according to their temperature in the second model shell T 2 . When the field is in cooling mode, the coldest BHEs are selected, and in heating mode, the warmest. For the BHEs to be activated, a nominal volume flow of 24 L/min is requested. This nominal volume flow is chosen based on Ref. [33] to ensure an efficient, turbulent heat transfer in the borehole heat exchangers.

Optimal control of valve positions and pump speed
In order to reduce the electric pumping power while providing the requested nominal volume flow for the selected BHEs, we need to solve an optimal control problem. Therefore, a model of the hydraulic system is required. According to DIN EN 60535-2-1 [34], the mass flow i m  of BHE i can be expressed as a function of the conductance of the valves C V and the pressure p valve,i that occurs at the valve. The volume flow can easily be derived from the mass flow by considering the density. We assume that the valve's pressure can be represented by the pump's pressure corrected by an individual parameter h i . This parameter accounts, for example, for the individual piping length of each BHE and other pressure losses. Furthermore, the conductance of the valve is a function of the valve opening ϕ, and the pump's pressure is a function of the pumps' speed and the total volume flow of the field tot V  . Thus, we can write the mass flow of BHE i as: The volume flow of the field is calculated as the sum of the volume flows of the individual BHEs: Additionally, the pump's electric power can also be expressed as a function of pump speed and volume flow.
The model was implemented in CasADi using look-up tables for the pressure and power of the pump and the valve conductance. These look-up tables were generated by datasheets from the manufacturers. To calibrate the model, the correction factors h i need to be parameterized. For this purpose, experiments in the real field were conducted. Here, the pump speed was varied between 30% and 100% with open valves. The measurements are then used in a least-squares optimization problem, which minimizes the model error to identify the 41 parameters h i . Fig. 8 shows the measured and modeled volume flows plotted against the piping length of each BHE for a pump speed of 100% and 40%, respectively.
In this figure, the piping length shown considers the length from the shafts to the borehole heads. This explains why the western shaft's volume flows, with the shortest piping between the shaft and the pump, are consistently higher than those of the others. Overall, the hydraulic model shows a good accuracy with a maximum deviation of 1 L/min, which is assumed to be sufficient.
To calculate the necessary pump speed and the valve openings to attain the requested mass flows req,i m  , the hydraulic model is used in the following nonlinear optimal control problem: The optimal control problem is implemented in CasADi [28] and solved with Ipopt [32] in less than 2 seconds.

Cloud-based implementation of the proposed control scheme
A commercial cloud platform is used to control the real geothermal field of the E.ON ERC main building with the proposed algorithm (see Fig. 9). The core of this cloud platform is a time-series database, which can be accessed via MQTT or HTTPS using a REST-API. To interact with the field, an edge device collects data from the building's automation layer, streams measurements into the cloud, and receives setpoints via the MQTT protocol. The pump speed and valve openings setpoints are written into the building automation system by the edge device using the BACnet (Building Automation and Control Network) protocol. Furthermore, the cloud platform provides real-time data visualization with Grafana 1 dashboards.
The developed control algorithm is split into two parts, which can run separately on two virtual machines. The first part contains the moving horizon estimator and the simulator and thus acts as a "digital twin" for the individual BHEs. The digital twin continuously downloads the latest measurements and streams simulation results back into the database with a sample rate of five minutes. Furthermore, the digital twin module checks the model accuracy of the active BHE and recalibrates these using the MHE if necessary. If no recalibration data are available, either because the BHE is inactive or operating with discontinuities, an operation request is submitted. The model predictive controller and the optimal control scheme are implemented in the second part. The controller part loads the heat exchanger measurements and the data provided by the digital twin first. Since no load forecast is currently implemented, the Fig. 9 Real-life implementation of the proposed algorithm using a cloud platform current load determined from measurements is used as a reference for the model predictive controller. For the same reason, the model predictive controller uses a smaller step size of 30 min compared to the simulation with a perfect forecast.
After calculating the required number of active BHEs, a decision is made on which BHEs to activate based on the simulated ground temperatures and operation requests. The optimal controller then calculates the valve positions and the pump speed to set the nominal volume flows while minimizing the pump power. The controller module uses the HTTPS protocol to submit the optimal set points to the cloud platform in the next step. Since the 41 st BHE was subsequently installed, its valve is not integrated into the building automation system and cannot be controlled remotely. For this reason, the 41 st BHE is disabled in the real-life demonstration.

Simulation Results
The proposed methodology is tested in a simulation to verify its functionality in the first step. Therefore, we use monitored data to initialize the simulation model and provide the simulation with data for the building side mass flow, inlet temperature as well as the thermal demand of the building. With the simulation model, we first analyze the short-term performance of the field. Afterward, we look at the long-term effects of the optimized operation strategy.

Improvements in short term operation
To evaluate the effects on short-term operation, the methodology is tested on data from mid-July to mid-August 2020. This cooling period is determined by high cooling demand of the building, which was met entirely by the geothermal field since the heat pump was under maintenance at this time. The results of a simulation with the optimized operation are presented in Fig. 10 (The subscript meas stands for measured). It becomes clear that almost consistently, less than half of the installed BHEs need to be operated to cover the load of the building ref Q  . This is remarkable, considering that the cooling demand in this period is the highest during the year. Only at a short period, characterized by a low building-side input temperature and a high power demand, all BHEs are used. The developed methodology visibly increases the temperature difference between flow and return temperature of the geothermal field. At the same time, at the end of the month, lower return temperatures are achieved, compared to the measured field operation, because the deactivated BHEs have improved regeneration capabilities.
The reduced number of active BHEs results in a significantly lower volume flow than the not-optimized operation, requiring significantly less pumping power. In the optimized operation, the pumps need, on average, 0.82 kW, which corresponds to a total consumption of electrical energy of 593.8 kWh from mid-July to mid-August. In the actual measured operation, the pumps run with a constant speed and an average electrical power of 2.55 kW. This adds up to total energy consumption of 1833.2 kWh in the considered time. Therefore, the developed control methodology saves about 67.6% in pumping energy.
Overall, the proposed control scheme shows good performance considering the short-term efficiency. A significant amount of pumping energy can be saved while satisfying the energy demand. Additionally, at the end of the simulated cooling period, we achieve lower return temperatures than the original operation mode, improving the cooling capability of the field at the end of the simulated period. Since the cooling demand in the considered time was the highest throughout the year, similar or even higher savings are expected for other periods. In the following, a discussion of the expected long-term effects is conducted.

Effects on long-term operation
To estimate the long-term effects of the optimized operation mode, the average temperatures of the RC model's capacities after one month of operation are considered. Fig. 11 presents the temperatures of the field after the cooling period described before for the optimized and the actual field operation. In both cases, the simulation is initialized with the same MHE estimates. Overall, the temperatures increase less in the optimized cooling mode compared to the original field operation. Furthermore, the temperatures in the second model shell T 2 , which serve as the basis for deciding whether to activate a BHE, are almost completely balanced. This also applies to the outermost model shell, whose temperatures are less scattered.
For further analysis of the field's operation, the temperature distribution in the field is considered. Fig. 12 shows the temperatures in the third shell of the borehole model for the single BHEs after the simulated period. The temperature distribution reveals differences between the three shafts. Especially the BHEs in the western and eastern shaft have heated up in the previous operation due to denser BHE placement. After one month of the original cooling operation, the temperature in shaft west is on average 0.6 K higher than in the southern shaft.
With the optimized operation scheme, this difference is reduced to 0.37 K.
The improved load balancing of the field becomes even more apparent when the temperature difference between the beginning and end of the simulated period is evaluated in Fig. 13. The temperature increase is around 0.5 K in regular operation for all BHEs. The optimized Temperature difference in the third shell of each BHE after one month of operation strategy prioritizes the southern shaft for cooling and thus allows BHEs in the other shafts to regenerate and even cool down in the considered period. This results in a lower temperature increase in the western and eastern shaft compared to the original strategy. Overall, the optimized strategy seems to benefit long-term field performance. It reduces the ground temperature change and reinforces a balanced field operation. In the following section, we evaluate the real-life demonstration of the presented control strategy.

Real-Life Demonstration
The presented control strategy was tested on the real ERC geothermal field in August 2021. This period was characterized by high cooling demand during the day. The first four days of operation are displayed in Fig. 14. This graph shows the thermal power of the field, the number of BHEs, the pump power, and the pump speed. The improved control strategy was activated on August 12 th at 7:30 am. The controller demonstrates demand-oriented behavior. Only a few BHEs are used at night and in the morning, while the number of active BHEs increases in the afternoon when the heat demand rises. Thus, the pump speed decreases to 50% to 55% during the nighttime.
Before activation of the model predictive controller, the field operated at a constant pump speed of 85% with all valves open, which led to an electrical power consumption of 3.96 kW. After activation, the average pumping power decreases to 1.15 kW, which corresponds to savings of 71% in electrical power consumption. The order of magnitude of the energy savings corresponds well to the simulation results. The volume flow of each BHE during the real-life experiment is presented in Fig. 15. It becomes clear that the BHEs 13 to 22 in the southern shaft and BHE 1 in the western shaft are preferred for cooling operation. This was also observed in the simulations. During the day, initially, the BHEs at the northern edge of the field (BHEs 4-7 and 39-40) are activated additionally to cope with the increased heat demand. As expected, the algorithm prefers BHEs for cooling located at the edge of the field and thus less densely packed. One explanation for the preference for the southern edge could be that the other shafts are located under an asphalted, dark area. This could have caused the near-surface parts of the BHEs to heat up more than the southern shafts' BHEs. What stands out is the frequent activation of BHEs 26 and 33, which are located in the middle of the field. This can be either explained by inhomogeneous underground properties like groundwater flows or, more likely, by estimation deviations of the MHE. Furthermore, Fig. 15 shows that the optimal controller for the valves and pump speed works well, as a flow rate of 23 to 26 L/min is maintained in the active BHEs with few outliers.
Overall, the real-life demonstration confirms the energy-efficient performance of the developed control strategy. The controller strategy saves up to 70% of pumping power while exploiting individual BHEs with better regeneration capabilities at the edge of the field.

Conclusions and Future Work
To reduce the consumption of fossil fuels, renewable energy sources like geothermal fields are becoming increasingly important. The main building of the E.ON Energy Research Center in Aachen uses a geothermal field consisting of 41 BHEs to meet its heating and cooling demand. The field is unique in terms of extensive monitoring and the capability to control each BHE individually. In this paper, we present a methodology to control borehole BHEs individually and evaluate the advantages.
The proposed methodology is based on a nonlinear model predictive controller, moving horizon estimation, and optimal control to estimate the ground temperatures. Since both techniques require a dynamic system model, an RC borehole model is presented and validated on measurement data. We show that a simple model with three thermal capacities is well suited to predict the borehole dynamics for up to one month.
The moving horizon estimator provides the nonlinear controller and a field simulation model with initial values. The NMPC scheme then optimizes the number of active BHEs to meet the building's current energy demand. In the next step, the BHEs to be activated are selected by the calculated subsurface temperatures.
The methodology is tested in a simulation study based on real-life measurements of a one-month cooling period and a real-life demonstration. Compared to the field's previously implemented constant hydraulic balancing, the proposed strategy performs significantly better in the simulation. With the optimized strategy, more than half of the BHEs can be deactivated most of the time, resulting in 67% savings in electrical pumping energy. Furthermore, the temperature change in the ground due to the cooling load decreases compared to the constant strategy. Therefore, the proposed methodology improves both the short-and long-term performance of the geothermal field. For the real-life demonstration, a cloud-based implementation of the strategy is developed. We achieve savings in pumping power of 71% during a four-day demonstration period. At this point, the controller uses BHEs preferably at the edge of the geothermal field with improved regeneration capabilities and gradually adds BHEs from the inner field when the thermal demand rises.
In conclusion, the presented methodology is a promising approach to improve both the short-and long-term performance of geothermal fields by reducing the energy consumption of the hydraulic pumps and improving the regeneration of the boreholes.
In our future work, we will look at the long-term performance of the algorithm in more detail. It is interesting to see if the algorithm leads to a more balanced operation as expected. In this case, the preferred operation of individual BHEs should become less distinct and the activation times more even. However, individual BHEs may be preferred due to different ground properties. Furthermore, the major drawback of this approach is the need for extensive hardware like actuators and sensors for each BHE. Therefore, we will examine if the investment in sensors and actuators is economically reasonable. We will investigate whether new geothermal fields can be dimensioned smaller when applying such an operation strategy to generate savings in construction costs. Finally, the interaction with the building energy system should be optimized. The latter, for example, could be achieved by a holistic system MPC that efficiently uses storage effects. The efficient usage of storage effects can reduce the peak loads and thus the dimensioning of the field even further.