Development of a Long-Term Operational Optimization Model for a Building Energy System Supplied by a Geothermal Field

In order to reduce energy consumption and CO2 emissions in the building sector, more and more renewable energy sources are integrated into energy systems. Especially geothermal fields combined with heat pumps are able to supply buildings with heat and cold at low carbon emissions. However, using geothermal fields as heat and cold source influences the ground temperature. Consequently, the ground temperature can change dramatically over a building’s lifetime, leading to less efficient operation of the energy system. Therefore, a sustainable operation is required to ensure the long-term efficiency of geothermal fields. In this paper, we develop an optimization model to derive operating strategies for an efficient long-term operation of a building energy system coupled to a geothermal field. The investigated energy system is the main building of the E.ON Energy Research Center in Aachen, Germany, which includes a heat pump, two boilers, a combined heat, and power unit, a glycol cooler, and a geothermal field with 41 probes. For each component, we develop energy-based sub-models, which are connected to form the overall system. The geothermal field is modeled by using a g-functions approach as well as a simplified resistance-capacitance approach. To achieve short computing times and realize an optimization horizon of several years, the optimization problem is formulated as mixed-integer linear programming (MILP). The developed model is optimized regarding two different objectives: the minimization of energy costs and the minimization of long-term temperature changes in the ground. Conclusions for an efficient and sustainable operation of the field, especially for the cooling supply, can be derived from the optimization results. It is shown that a state of equilibrium should be aimed to achieve an energy-efficient operation, in which the temperature of the field is close to the initial ground temperature.


Introduction
The building sector is responsible for 28% of the energy-related CO 2 emissions [1]. One possible way to reduce carbon emissions is the integration of renewable energy sources in building energy systems. Especially, the use of geothermal fields as an energy source has increased over the last decades [2]. The fields are usually connected via heat pumps to the energy system and can serve as a source for both heating and cooling by injecting heat into the ground or extracting heat from the ground. However, the extracted or injected heat flows affect the ground temperatures around the boreholes [3]. In particular, if the amount of heat extracted and injected is unbalanced, the ground temperature can decrease or increase dramatically over a building's lifetime, which results in reduced efficiency of the heat pump and the overall energy system [4,5]. To avoid this, a sustainable operation is required to ensure the long-term efficiency of geothermal fields. To improve the performance of geothermal heat sources, Paly et al. investigate an operational optimization for a field with 25 boreholes [6]. Compared to a non-optimized operation, the optimized operation reduces the average temperature decrease in the ground by 18%. Bayer et al. study a geothermal system that can also be used for cooling [7]. The results show that an optimized operation can reduce the temperature change over 15 years compared to an equal load operation. Additionally, Capozza et al. figured out that a balanced heating and cooling demand is beneficial for a sustainable long-term operation [8]. Especially using the borehole field as long-term energy storage is promising [9,10]. Further, Li et al. consider a load optimization as well as the geometric arrangement of the boreholes to improve the efficiency [11]. However, these studies focus on the load distribution for the boreholes, while the whole building energy system, which includes components such as heat pumps, boilers, air-handling units, or chillers is not part of the optimization. However, an improved operation of the building energy system can lead to savings of 20%-30% for non-residential buildings [12]. Thus, the overall energy system including the geothermal field has to be investigated for reducing carbon emissions over long time horizons.
In this paper, we present an optimization model to derive operating strategies for an efficient long-term operation of a complex building energy system coupled to a geothermal field. The investigated energy system is the E.ON Energy Research Center's main building in Aachen, Germany, which includes a heat pump, two boilers, a combined heat, and power unit, a glycol cooler, and a geothermal field with 41 probes. Generally, the proposed approach can be transferred to other buildings and the developed models can be used for that purpose. However, if other buildings include further components that are not considered in this paper, these components need to be modeled according to the proposed energy-based approach. The structure of this paper is as follows: in Section 2, the main building of the E.ON Energy Research Center is described in detail and mixed-integer linear models for the components are presented. The models are based on energy flows between the components. In Section 3, three different modeling approaches for the geothermal field are compared and the results of the long-term operational optimization are presented and discussed. Section 4 concludes this work.

Methodology
This chapter presents the investigated building and its energy system and the modeling approach based on mixed-integer linear programming (MILP).

The E.ON energy research center main building
The main building of the E.ON Energy Research Center is an office building with additional laboratories and has an area of 7222 m 2 [13,14]. The building includes a multi-functional energy system, depicted in Fig. 1.
The system provides heat for consumers in a high temperature circuit (HTC) and low temperature circuit (LTC) as well as cold on two different temperature levels. The high temperature consumers include two Fig. 1 Energy system of the E.ON Energy Research Center main building: The core of the system is a heat pump for cold and heat supply. Two boilers and a chp unit provide additional heat and a geothermal field can be used as a heat and cold sink.
air-handling units that supply the building with fresh air. Each air-handling unit has a nominal air volume flow rate of 12.000 m 3 /h. Further, the building can be heated up and cooled down by a concrete core activation. The low temperature circuit supplies the concrete core activation for heating. The first cold circuit (cold 1) distributes the cold to the air-handling units and the concrete core activation, while the second cold circuit supplies the static cooling demand of server rooms. The heart of the system is a turbo compressor heat pump with a nominal electrical power of 51 kW. The maximal heating power is 254 kW and the maximal cooling power is 203 kW. The coefficient of performance ranges from 4 to 5.7 depending on the part load rate and temperature level. The heat pump is connected to a hot and a cold water storage with volumes of 4 m 3 and 5 m 3 respectively. A combined heat and power unit (chp) based on an internal combustion engine with 30 kW and two natural gas condensing boilers with 120 kW each are used to provide heat on the high temperature level. The heat from the high temperature circuit can be transferred to the low temperature circuit via a heat exchanger (HX) if needed.
In the case of high cooling demand, an air cooler can be used to dissipate heat from the hot side of the heat pump.
If the outside temperature is below 10°C, the air cooler can be directly used for cooling the cold consumers. Further, a geothermal field consisting of 41 probes serves as a heat source for the heat pump as well as a heat sink for the cold consumers. Each probe includes a double-u pipe and has a depth of 100 m. The probes are arranged in three shafts around the building: a west, a south, and an east shaft. The three shafts are connected via one heat exchanger to the building. The fluid in the hydronic circuit of the probes is a water-glycol mixture and the volume flows, the inflow and outflow temperatures at each probe as well as the mixed temperatures and volume flow of the three shafts at the heat exchanger are measured. The ground temperature has to be kept in a certain temperature range to achieve a sustainable and efficient operation of the geothermal field and the overall system. In the next section, an optimization model of the system is presented, which is used to optimize the long-term operation.

Optimization model using mixed-integer linear programming
In this section, we develop the model equations for the optimization problem. Each component is modeled separately based on energy balances and connected to the corresponding circuit via heat flows. To realize an optimization with a time horizon over several years, the influence of mass flow rates and the energy consumptions of hydronic pumps are neglected.

Geothermal field
A common approach to model boreholes of geothermal fields is based on three-dimensional response factors, also called g-functions, as presented in Ref. [15]. This approach assumes an infinite line source heat flow [16]. The change of the borehole temperature from the undisturbed ground temperature ΔT due to an average heat flow q  at a discrete-time period depends on the average heat flows per borehole length and a factor g: The g-function depends on the time t and geometric parameters p geom . The parameters include the geometric features of the borehole, such as radius and length. The parameter k s denotes the ground thermal conductivity. As described in [15] multiple boreholes can be modeled using the superposition principle. Thus, the g-function approach results in a system of linear equations for multiple boreholes.
In addition to the g-function-based approach, we investigate a thermal resistance-capacitance-based (RC) model as presented in Ref. [17]. However, we use a simplified approach with one node, resulting in one capacity and one resistance. Thus, the geothermal field is assumed to be a mass with a homogenous temperature: The parameter C represents the thermal capacity of the geothermal field, and the parameter R is the resistance to the undisturbed ground. The variable T denotes the temperature of the field; T 0 is the temperature of the undisturbed ground and Q  the heat flow added or extracted from the field. The parameters C and R are estimated in a calibration with measured data of the field, resulting in C=478.8×10 9 J/K and R=0.001 31 K/W.
In Section 3.1, we compare the g-function approach with the capacity-resistance model. Further, we investigate the g-functions with one representative borehole considering all 41 probes. Therefore, the following three models are formulated and compared: (1) probe g-function model; (2) 41-probe g-function model; (3) node RC-based model.

Heat pump
The heat pump model considers the enthalpy flows of the low temperature circuit and the cold circuit. An energy balance around the heat pump results in Eq. (3). Here, ev Q  denotes the evaporator heat flow on the cold side; con Q  denotes the condenser heat flow on the hot side of the heat pump, and P el is the consumed electrical power.
The efficiency of the heat pump is characterized by the coefficient (COP) of performance according to Eq. (4).
The COP is variable and depends on the temperature of the condenser T con and evaporator T ev . Thus, Eq. (4) is a non-linear expression. To allow a fast computation, we approximate the temperature dependency with a linearized formulation: The parameters c 1 to c 4 in Eq. (5) are obtained by a regression of measured data. The binary variable δ on is one, if the heat pump is operating, or zero otherwise. The product of a binary variable and a continuous variable is linearized according to Ref. [18]. We assume that the evaporator and condenser temperatures are equal to the storage temperatures of the neighboring hot and cold storages.

Storages
The cold and hot water storages that are connected to the heat pump are modeled as one volume with a homogenous storage temperature T s each. The energy balance around a storage leads to: Here, c w denotes the thermal capacity of the water and ρ w denotes the density of the water. The volume V s is 4 m 3 for the hot storage and 5 m 3 for the cold storage. In the case of the hot storage, the incoming heat flow in Q  corresponds to the heat flow provided by the condenser of the heat pump. The outgoing heat flow out Q  supplies the low temperature circuit since the hot storage is also connected to the low temperature circuit (cf. Fig. 1). In the case of the cold storage, the incoming heat flow in Q  comes from the cold circuit and out Q  is the evaporator heat flow of the heat pump.

Boilers and combined heat and power unit
The condensing boilers and the chp are modeled assuming a constant efficiency η th in order to calculate the gas consumption gas m  . The parameter HHV denotes the higher heating value. The efficiency for the boiler is 0.984 and the thermal efficiency of the chp is 0.62.
The chp further produces electricity P el depending on the electrical efficiency η el and the consumed gas according to Eq. (8): The electrical efficiency of the chp is assumed to be a constant value of 0.3.

Air cooler
The air cooler can be used to dissipate heat from the hot side of the heat pump (re-cooling) if the heat pump is needed for cooling and no heat is required. Additionally, heat from the cold circuit can be dissipated to the air if the outside temperature is below 10°C. The outside temperature is a parameter and not a decision variable of the optimization. Thus, the following linear equation without any integer variables can be formulated depending on the outside temperature T amb and the maximal heat flow cooler,max Q  : amb cooler cooler,max cooler if 10 C : 0

Consumer
To calculate the energy demand, the building is split into two zones supplied by an air-handling (ahu) unit and a concrete core activation (cca) each. The air-handling units are connected to the high temperature circuit for heating and connected to the cold 1 circuit for cooling. The concrete core activation is connected to the low temperature circuit for the heat supply and is connected to the cold 1 circuit for the cold supply. The zones are modeled based on an RC approach, where the air and the concrete are a capacitance each. The heat transfer between the room and the concrete and the heat transfer from the room to the environment is modeled as resistance and depends on the temperature of the concrete T cca and the heat transfer coefficient k cca between the concrete and the zone. Additionally, the zone is supplied by fresh air with the constant volume flow rate supply V  at the temperature T supply . The heat loss to the ambient is based on the heat transfer coefficient multiplied with the outside wall area k wall and the ambient temperature T amb . Further, the internal gains emitted by persons and electrical appliances int Q  and the solar radiation sol Q  is included in the energy balance of the zones: Here, the parameter V zone is the air volume of the zone; ρ air denotes the density of the air and c air is the specific heat capacity of the air. The zone temperature T zone has to be in the range of 21°C to 23°C. We assume constant internal gains int Q  of 5 W/m 2 and a heat transfer coefficient to the ambient k wall of 0.3 W/(K·m 2 ). The required energy to heat up or cool down the supply air temperature of the air-handling units can be calculated by an energy balance around the air-handling unit: Here, h,ahu Q  is the heating, and c,ahu Q  the cooling energy flow that has to be provided by the energy system. The parameter ε denotes the efficiency of the heat recovery system in the air-handling unit, which is 0.75. The heat flows of the concrete core activation h,cca Q  and c,cca Q  are calculated in the same way: The parameter c cca denotes the specific thermal capacity of the concrete and the parameter m cca is the mass of the concrete.
The cold 2 consumer includes the server rooms and has a constant cooling demand of 25 kW.

Objective function
The model equations above form the basis for the operational optimization of the energy system concerning different goals. Optimization aspects can be, for example, reducing costs, reducing CO 2 emissions or, a sustainable operation. In this work, two objectives are considered: Minimizing energy costs and minimizing the change of the ground temperature over the runtime. We scalarize the formulation to a single-objective optimization to solve the multi-objective problem. For this purpose, a combined objective function J is formulated, which allows a trade-off between efficiency and sustainability: The first term sums up the energy costs for all components C t over the optimization horizon, whereas the second term penalizes the temperature difference ΔT abs between the ground temperature at the first time step and the mean ground temperature of the last year. We assume specific costs of 0.07 EUR/kWh for the consumed gas and 0.22 EUR/kWh for the consumed electricity. The produced electricity of the chp is sold for 0.12 EUR/kWh. The weighting parameters α and β are chosen in three different ways according to Table 1. The cost-based objective minimizes the total operating costs and the temperature-based optimization minimizes the long-term temperature deviation ΔT abs in the ground. The multi-objective optimization combines both aims. It is necessary to scale the costs to be in the same order of magnitude as the temperature deviation. The costs are estimated to be in the range of 1×10 5 to 5×10 5 EUR whereas the temperature deviation is estimated to be in the range of 0 to 10 K. Thus, we choose α=1×10 -5 and β=1 for the multi-objective function.

Validation of the geothermal field models
To validate the three modeling approaches for the geothermal field, we simulate a period of three years (2015 to 2018) and compare the results with measurements of the main building of the E.ON Energy Research Center. The measurements include the inflow temperature, the return temperature, and the volume flow of the water-glycol mixture that flows through the probes. We assume that the return temperature corresponds to the mean ground temperature of all probes. The heat flow, which is added to or extracted from the ground is used as an input for the simulation and is calculated based on the temperature and volume flow measurements of the field. The undisturbed ground temperature is assumed to be 10°C. The graphical comparison of the three models is shown in Fig. 2. Here, the measured data corresponds to the mixed outflow temperature of the water-glycol mixture of all probes. The temperature of the RC-based model corresponds to the temperature of the capacitance. In the case of the g-function models, the temperature of the one representative probe is for the 1-probe model and the mean temperatures of all probes are for the 41-probe model. As can be seen in Fig. 2, the measurements have some anomalies in the winter of 2017. At these times, the geothermal field has not been operated. Thus, the volume flow of the water-glycol mixture is zero and the temperature measurements run against the ambient temperature of the building. The resistancecapacitance-based 1-node-model deviates from the measurement. Especially in the first year, the temperatures were underestimated. The measurement is incomplete in the first half of 2016, leading to a nearly constant interpolated temperature profile. In 2017 the temperatures of the 1-node-model are closer to the measurements. The general behavior of the RC-based model corresponds to the measured data in the long term. However, short-term effects like the temperature increase and decrease in the first year are not covered. This behavior is due to considering only one node. This approach assumes a homogenous temperature and does not distinguish between the ground temperature in the direct surrounding of a probe and the ground temperature in a few meters distance. We assume that the accuracy could be improved by using multiple nodes in the RC model. The g-function-based approaches are able to cover the dynamic behavior of the field over the whole simulation period. However, the 41-probe model is less accurate than the 1-probe model, although the 41-probe model is more detailed. Especially in the last year, the 41-probe-model overestimates the temperature. One reason for this could be that the heat flow is distributed equally on all probes in the simulation. Thus, the same thermal power is imposed on all 41 probes. However, the ground temperature and the inflow temperature determine the heat flow in the real system leading to different heat flows for each probe in reality. The inflow temperature is equal for all probes; thus the thermal power is different between the probes if the surrounding soil has different temperatures. Furthermore, it can be assumed that the ground has not a uniform temperature in reality. However, a homogeneous soil temperature is assumed at the beginning of the simulation. Overall, the validation shows that the 1-probe model of the g-functions, assuming equally distributed heat loads on the probes, leads to the most accurate results. Additionally, the calculation effort is less compared to the 41-probe-model. Therefore, the 1-probe-based g-function approach is used for further optimization.

Optimization results
The operation optimization is performed for a horizon of 15 years with a step size of five days. Fig. 3 shows the resulting costs and the mean temperature deviation of the ground in the last year for the cost-based, temperature-based and multi-objective optimization. The cost-based optimization leads to the lowest costs with 329 146 EUR. However, the mean ground temperature increased by 2.06 K after 15 years. The temperaturebased optimization leads to the highest costs (379 967 EUR) but has no temperature change after 15 years. The multi-objective approach leads to 1% higher costs (332 865 EUR) than the cost-based optimization. However, there is no change in the mean ground temperature in the last year of the optimization as well.

Fig. 3
Comparison of the energy costs and ground temperature deviation for the cost-based, multi-objective and temperature-based optimization Fig. 4 shows the yearly mean ground temperature curves in the upper plot and the ground temperatures at every time step in the lower plot for the cost-based, multi-objective and temperature-based optimization. Since the temperature-based optimization only considers the deviation of the start temperature and the mean temperature in the last year, the ground temperature does not follow any logical pattern in this case. The energy system has only to fulfill the demands of the consumers. Thus, an infinite amount of different operation decisions over the optimization horizon would lead to the same objective value of 0 K temperature deviations. We assume that considering the temperature deviation for every year in the objective function would lead to more deterministic results.
In contrast, the cost-based and multi-objective approach results in a similar yearly periodic temperature curve except for the last year. The cost-based optimization runs asymptotically against a yearly mean ground temperature of 285.5 K in the long term. The ground temperature oscillates every year (lower plot), and the maximal ground temperature is nearly constant and does not exceed 288 K, while the minimal ground temperature slightly increases every year. Thus, the cost-based approach leads to a slightly increasing yearly mean ground temperature. This behavior can be explained by the fact that a high decrease in the ground temperature would also result in a less efficient operation of the heat pump, and a high increase would reduce the cooling capacity in the summer. Additionally, an increasing ground temperature causes a higher heat flow from the geothermal field to the undisturbed ground, which counteracts the temperature rise in the long term. Thus, a state of equilibrium in the yearly mean ground temperature will be reached in the long term. The long-term temperature deviation of 2 K of the cost-based optimization is acceptable and would be an improvement compared to the current operation of the building shown in Fig. 2. In contrast, the multi-objective optimization decreases the mean ground temperature to 283.15 K in the last year to minimize the temperature penalization in the objective function. We assume that this rapid temperature change in the ground would have no negative effects as long as local policies, e.g. for groundwater temperatures, are not violated.

Fig. 4
Ground temperature curves for the cost-based, temperature-based and multi-objective optimization The heat flows of the components for the cost-based and the multi-objective optimization for the last four years are illustrated in Fig. 5. Both optimizations result in a similar yearly repeating operation of the energy system. In the winter month at the beginning of each year, the heat demand is covered by the chp, the heat pump, and the boilers. The chp and the boilers supply the air-handling units in the high temperature circuit. Fig. 5 shows that the chp runs at maximal power of 30 kW if the heat demand is high enough. This is primarily the case during the winter. If the maximal power of the chp is reached but further heat is required, the boilers are used to cover the remaining high temperature energy demand. Here, the total required energy demand of the high temperature circuit can be above 140 kW to heat the supply air. Additionally, the heat pump is used to supply the low temperature circuit and thus the concrete core activation, which is additionally used to heat the building.
Further, the heat pump covers the cooling demand of the server rooms of the cold circuit 2 during that time. However, the waste heat of the server rooms is not enough to supply the evaporator of the heat pump with a sufficient amount of energy. Thus, the geothermal field is used as an additional heat source for the heat pump and heat is extracted from the geothermal field. Therefore, the ground temperature increases in winter (compare Fig. 4). When the heat demand decreases after the first month of each year, the heating and cooling power of the heat pump is reduced and the geothermal field is used as a heat sink to cover the cooling demand of the cold 2 circuit. During the summer in the middle of each year, only the geothermal field is used to provide the cold for both cooling circuits, which results in a rising ground temperature. However, an increased ground temperature leads to a more efficient operation of the heat pump during the heating period. Thus, we can conclude that the field should be used as long-term storage serving as a heat source in the winter and heat sink in the summer for the cost-based and multi-objective operation.
The main difference in the operation between the cost-based and multi-objective approach can be seen in the last year of the optimization. While there is no changing behavior in the cost-based operation, the multi-objective optimization extracts heat from the geothermal field and dissipates the heat to the ambient with the air cooler (free cooling) to decrease the ground temperature. This recovery process starts at the end of the heating period in the year 14 because an earlier temperature decrease of the field would lead to a decreased efficiency of the heat pump. Since the air cooler can only dissipate heat to the ambient if the air temperature is below 10°C, the field cannot be cooled down by the air cooler in the summer of the last year. Thus, the recovery process of the field starts again at the end of the last year. Taken together, the optimization results show that a cost-based operation leads to constant mean ground temperature deviation if the time horizon includes several years (see Fig. 4). In this case, the ground temperature is kept in a certain range to realize an efficient operation of the heat pump. Thus, the cost-based operation leads to a sustainable operation as well. In contrast, the optimization of the long-term deviation of the ground temperature leads to no temperature change at the end of the optimization horizon. However, the operation costs are significantly higher. The best results can be obtained by a multi-objective optimization including the costs and the long-term temperature deviation in the ground.

Conclusions
Geothermal fields combined with energy systems can supply buildings with heat and cold at low carbon emissions. However, extracting and adding energy to the field influences the ground temperature. For a sustainable operation over the building's lifetime, an efficient long-term operation is required for the overall energy system. This paper presents an optimization model for a complex energy system coupled to a geothermal field, which allows deriving long-term operation strategies. We use the main building of the E.ON Energy Research Center as a case study, which includes a reversible heat pump, boilers, a chp, and a geothermal field consisting of 41 probes. We presented mathematical models for all system components based on mixed-integer linear programming (MILP) in order to realize an optimization horizon of several years. Further, two g-function based models and one resistance-capacitance model for the geothermal field are investigated and compared with measured data.
The results show that the g-function-based models of the geothermal field lead to a smaller deviation from the measurement than the considered RC-based model. Furthermore, the difference between modeling all 41 probes compared to a model with one representative probe is marginal for the g-function approach. Thus, the g-function model with one representative probe is used for the long-term operational optimization. The developed model is optimized regarding a cost-based objective, a long-term ground temperature-based objective, and a combination of both. The results show that the cost-based objective leads to the lowest costs of 329 146 EUR, but also leads to a temperature increase in the ground. However, the change of the ground temperature is below 2.1 K in the long-term and the temperature curve of the ground is periodic for each year. The temperature-based objective keeps the ground temperature at the same level, while the costs increase by 15% compared to the cost-based optimization. The combined optimization leads to low costs, which are 1% above the costs of the cost-based optimization, while no temperature deviation occurs in the ground at the end of the optimization horizon. The analysis of the cost-based operation and multi-objective optimization shows that the geothermal field should be used as a heat source for the heat pump in winter and as a heat sink to cover the cooling in summer. In future work, we will develop a control strategy based on the optimization results and investigate the control in a dynamic simulation.