Geochemical implications of production and storage control by coupling a direct-use geothermal system with heat networks

This paper outlines a method in which the heat production of a geothermal system is controlled in relation to the demand from a district-heating network. A model predictive control strategy is designed, which uses volume measurements in the storage tank, and predictions of the demand, to regulate the production of the geothermal system in real time. The implications of such time-varying production for the reservoir are investigated using a 2D reactive transport reservoir model. As a case study, the Groningen geothermal project is considered. The numerical data generated by the controller, in closed loop with a modelled district-heating network, are used as inputs for the reservoir simulations. The latter make use of discrete parameter analyses to evaluate the effect of pressure depletion, reservoir permeability, flow rate, re-injection temperature and injection pH on the geothermal reservoir, and also mitigate possible risks during development. Using a model predictive control does not create adverse geochemical effects in the reservoir; instead, the controller is able to improve the efficiency of the geothermal heat extraction. The findings pave the way for stronger integration between elements of heat networks and a more sustainable development of geothermal resources. 2017 The Author(s). Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).


Introduction
Currently heat demand constitutes the largest part ($78%) of household energy consumption in the EU [1]. Renewable heat from geothermal systems can aid the reduction of CO 2 emissions associated with conventional heat production [2][3][4]. Geothermal district heating systems have been used since the 14th century [5], but were challenged by the wide use of cheap fossil fuels. However, due to global warming and a revived focus on renewable energy sources, district heating and heat energy networks are gaining importance in the provision of renewable energy [6][7][8]. The spatial topology and integration of a heat network in an urban setting has been recently analyzed in a comprehensive manner [9]. The complexity and challenges related to geothermal heat distribution have been previously outlined [10], while the efficient production and use of geothermal resources has been identified as an important aspect of their sustainable development [11]. Such insights are relevant for all geothermal fields. Nonetheless, the nature and impact of the challenges cannot be fully generalized and should be addressed by means of project-level studies [11].
In a conventional setup, a geothermal system provides the baseload heat supply, while backup systems cover any excess demand [8,12]. Recently, it has been shown that periods of high and low geothermal heat production can lead to more sustainable utilization of the geothermal resource [13]. Global historical data on direct-use geothermal systems suggest that there is a capacity factor drop over time [14]. This drop could be partly attributed to a better utilization of the geothermal heat produced by means of coupling supply and demand.
It is important to match the heat supply and demand [3] in district heating systems, and daily discrepancies can be bridged with the use of a storage component [8,15], while the geothermal production level can be adjusted to match the seasonal changes. However, applying a seasonally variable production rate to the geothermal system can have several consequences at reservoir level, among which are changes in chemical composition, pressure and possibly also temperature. Moreover, a seasonally variable production rate could also affect the cold front breakthrough time of the reservoir, when compared with a constant production level. Salt dissolution or precipitation can affect reservoir permeability [16] and therefore needs to be addressed with location-specific modelling. Recent experimental studies presented geochemical interactions [17,18] in conduction-dominated geothermal settings [19]. Moreover, chemical implications, in the form of salt precipitation during geothermal production using CO 2 as the energy carrier, have also been recently highlighted [20]. It was demonstrated that a dynamic production rate can lead to clogging of the reservoir due to salt precipitation. To avoid this, the production rate of change should be constrained. Obtaining realistic values for these constraints is often difficult, as they are dependent on the characteristics of the geothermal system and therefore case specific. Furthermore, a geothermal system has upper and lower production constraints determined by reservoir properties and engineering specifications.
In order to satisfy all previously mentioned constraints while supplying a time varying heat demand, a storage device can be used to shift loads in time [15]. To provide the geothermal system with time-varying production rates, a controller should be designed that takes the production and storage capacity constraints into account. In case the demand has a periodic structure, an internal model controller can be used, as is shown in [21]. In [22] several other controller designs are presented that do not require a periodic demand, among which well-tuned propor tional-integral-differential (PID) controllers and model predictive controllers (MPC) are the most promising. The PID controllers are very easy to implement and guarantee stability, but cannot guarantee that the constraints are always met. Conversely, an MPC does have the capacity to guarantee that the constraints are always satisfied but the stability of these controllers is hard to prove. Moreover, these MPC mostly rely on ad-hoc tuning and experimental analysis [23]. Despite these drawbacks, MPC received a lot of attention [23][24][25][26] and are also applied to the control of pressure control of geothermal systems [27] and thermal energy storage for buildings [28].
An MPC solves an optimal control problem over a finite discrete time horizon, returning a sequence of control inputs of which only the first one is implemented. After this implementation, the process is reiterated using a new finite horizon that is shifted one step forward. Since the future demand is often unknown, a prediction can be made to solve the optimization problem. These predictions can be based on, for example, historical data and weather predictions. Also, a dynamic model of the system that is to be controlled is required to implement an MPC. Such a model relates flow rates and storage level [29], and is well suited to modelling a districtheating network.
The importance of direct-use, deep geothermal systems for renewable heat supply is recently highlighted [2,4]. However, for the simulation of such geothermal reservoirs, the implications and complexity of the geothermal system are usually simplified [2] or not discussed [4]. Moreover, demand pattern changes are either not taken into consideration [30], or only described by a maximum [31,32] or annual demand level [3,9]. Additionally, direct use geothermal systems often exhibit risks that are difficult to estimate, particularly at the early phases of development [12]. The effect of a variable geothermal production resulting from the coupling between demand and supply to the reservoir geochemistry has not been studied before, possibly impeding the application of such systems.
In this work, the heat demand is delivered using a district heating system that includes a storage device. An MPC is designed that regulates the production of the geothermal system. Although the design of MPC is not new, such a design has not been applied before to a geothermal system. The controller uses a storage level measurement and demand prediction as inputs and takes constraints into account for the production level, change in production level and storage level. A realistic, yearlong demand pattern for an equivalent of 10,000 households is used as input.
The resulting MPC production levels for the geothermal system have a realistic time-varying behavior, which is used as input for geochemical reservoir simulations. In this paper both the implications and complexity of the geothermal system and the changes in demand pattern are taken into account. Moreover, to take uncertainties in reservoir pressure depletion, permeability, flow rate, injection temperature and pH into account, multiple simulations are performed which helps to mitigate possible risks during the development phase.
The reservoir simulations are performed using a 2D model to obtain several insights. Firstly, it is investigated whether the geothermal doublet is able to provide the demanded energy (i.e. feasibility of delivery). Secondly, the long-term effects of a variable, demand driven, seasonal production pattern on the reservoir behavior (i.e. pressure, power, permeability and chemical changes) is compared to constant production rate data. Lastly, the interaction between the chemical and physical parameters of the reservoir is outlined. The analysis makes use of the Groningen geothermal project (NE Netherlands) data and features.
The paper is structured as follows: in Section 2 a controller is designed that regulates the production of the geothermal system, followed by a description of the characteristics of a 2D reservoir located in Groningen. In Section 3 an analysis of the performance of the controller in closed loop with the district heating network is carried out. This is followed by an analysis of the geochemical implications for the reservoir. Finally, a discussion of the findings, and the conclusions that can be drawn are presented in Sections 4 and 5, respectively.

Methods and background
The possibility and implications of a time-varying production of a geothermal system that is controlled in real-time are evaluated. To this end, an MPC is designed that regulates the production of the geothermal system. This controller is connected to a modelled storage device in order to analyze its performance. The demand pattern is predicted based on historical demand data. The MPC uses such predictions, in combination with measurements from the storage device, as its inputs. Additionally, the MPC takes into account the predetermined limitations of the reservoir in the form of constraints on the change in production rate. This control structure is depicted in the upper part of Fig. 1.
To investigate the consequences of this time-varying production, a 2D reservoir model of the Groningen geothermal project is used. The production levels that are generated using the heat network were converted from an hourly to a monthly demand and used as input for a 50-year reservoir simulation. The lower part of Fig. 1 shows the model of the reservoir for which 243 scenarios were considered, each one corresponding to a different uncertainty parameter combination. The results of the reservoir behavior in terms of pressure, permeability, chemical and power changes are investigated.

Controller
In this section, a controller is designed to regulate the production of a geothermal system. It is designed such that it can be applied to any setup that includes a district-heating network, storage device and geothermal system. The controller can also be tuned in order to meet case specific requirements. Firstly, the demand, the model of the system, the constraints and demand predictions are introduced and discussed, after which they are incorporated into the design of the MPC.

Demand
A geothermal system that is connected via a heat exchanger to a storage tank was considered. This storage tank is connected to a network with several consumers. These consumers have an aggregated demand d a ðkÞ given in watt (W) where k denotes the discrete time step. However, the geothermal system has a maximum production rate u max and, in order to guarantee that the demand can be met with an admissible production rate, the demand is split in the demand covered by the geothermal system and the excess demand that is covered by a backup system such as a gas-fired boiler or biomass incinerator. For simplicity, we consider the case in which the excess demand d b ðkÞ can be ignored, which occurs when the backup systems automatically satisfies the excess demand. For optimal coordination of production in the multiple producer case, the reader is referred to [29] and to Section 2.1.6. In order to obtain a discrete time model for the storage device, the demand is converted from W to m 3 /h. This is done using: where d is the demand covered by the geothermal well in m 3 /h, C p is the specific heat in J/kgÁK, q is the density of the water in kg/m 3 and DT is the temperature difference between the production and injector well in the reservoir in Kelvin (K). It is assumed that temperature difference DT is independent of the time step k, which implies that d g ðkÞ depends linearly on dðkÞ. Moreover, the temperature difference in the reservoir and district-heating network are taken to be equal by assuming a heat exchanger with no heat loss and equal flow rates on both sides. Having defined our demand, the model of the storage is introduced.

Storage model
A storage model is considered, which is connected to a district heating grid with demand dðkÞ and a geothermal well with productionuðkÞ (in m 3 /h) at time k. The heated water contained in the pipes of the district-heating network can act as a buffer, but since it is assumed that DT is constant, the temperature dynamics of the fluids in the pipes can be neglected. The resulting model for the storage device takes the form of a difference equation and is given as where V is the volume of the storage tank. The controllable input will be the production u while d is considered to be a disturbance. Same as in (3), the input uðkÞ can be converted back to Watts using: where u g ðkÞ is the production of the geothermal system given in W.

Constraints
In order to keep the volume of the storage tank between the minimal and maximal volume capacity, it is required that: Furthermore, the pumps in the geothermal system have a minimal and maximal capacity. Therefore, the following is also required: where u max is defined as in Section 2.1.1 Demand. Lastly, the porosity, and therefore also the pressure drop, in the geothermal reservoir (see Section 3.2 Geochemical reservoir simulations for more details) can be affected by rapid changes in the flow rate. For this reason, it is desirable for the change in production rate of the geothermal system to have upper and lower limits, to avoid rapid changes, therefore: Du min 6 uðk À 1Þ À uðkÞ 6 Du max ð8Þ for all k P 1. Let k 0 denote the current time step with corresponding input uðk 0 Þ, k 1 ¼ k 0 þ 1 be the next time step and let N c 2 N be a control horizon. It is required that model dynamics (4) and constraints (6)- (8) are satisfied for all time k 2 fk 1 ; . . . ; k 0 þ N c g. In order to write the model and constraints in a more compact form the following state, input and auxiliary variables are defined: where 1 denotes the vector of all ones. Using this notation, (6) and (7) can now be written for k 2 fk 1 ; . . . ; k 0 þ N c g as where the inequality relations are to be intended component-wise.
Moreover, again to facilitate a more compact form, the matrix S 2 R Nc ÂðNc þ1Þ is defined as: from which it follows that: which implies that (8) for k 2 fk 1 ; . . . ; k 0 þ N c g reads as: Again using (11) it is possible to write (4) for k 2 fk 1 ; . . . ; k 0 þ N c g as:

Predictions
Note that at time k 0 , d is not known since it represents future demand. To overcome this, a prediction b d is used for the implementation. These predictions can be based on weather forecasts, historical demand and/or known changes in the network such as planned maintenance or newly connected consumers. The controller will in general perform better when more accurate predictions are utilized.

Objective function
In order to keep the storage from draining or flooding, a performance measure is introduced that penalizes a storage deviation from a predetermined reference value, which is taken as half the storage size. Furthermore, besides the hard constraint (13), the rate of change of the control input is also penalized in order to obtain a sufficiently smooth input. The overall performance measure is defined as a weighted sum of these two performance measures and is given by: where a is a positive dimensionless parameter. Minimizing (16) implies that storage deviations and rate of change of the input are penalized. Which one of the two is penalized more aggressively can be decided by attributing a value to a. Once the predictions, constraints and objectives are finalized, the controller can be designed.

Controller design
As before, let k 0 denote the current time step. First the controller makes a prediction b d of the disturbance over a finite control horizon N c . This prediction is used to solve an optimization problem that results in a vector of inputs u of which only uðk 1 Þ is implemented. After the initial implementation, the system evolves under the conditions caused by actual disturbance, whereupon the time window is shifted to the next time step (k 0 is set to k 1 ). These procedures are then repeated for all further time steps. During each iteration, the following optimization problem is solved: In case one wants to control not only the geothermal system but also the backup system, the backup production costs can be included in (17) and the possible constraints associated with such a backup system can be included in (18). In this way the MPC controls both the geothermal system, as well as the backup facility by minimizing the production costs. However, as mentioned before, control of any backup system is beyond the scope of this paper.

Reactive flow model
In order to evaluate the impact on the reservoir of fluctuating production, a reservoir model was used. Since project-level studies are needed to understand geothermal systems and gain insights, the model is based on data from the Groningen geothermal project [33]. The reactive flow model is discussed in further detail below.

Reservoir model
The model was built in the PetraSim pre-post processor [34], using the TOUGHREACT [35] code for chemically reactive, nonisothermal flow and utilizing the EOS1 (Equation of State) [36]. The TOUGHREACT code uses space discretization through Internal Finite Difference (IFD) [37]. The solver uses a sequential iteration approach for coupling between fluid and reactive flow [38].
The chemically reactive flow model used is largely based on the model geometry and characteristics of the Groningen 3D reservoir model [33], adapted to 2D. Overall characteristics and architecture are presented in Table 1 and Fig. 2, respectively. It should be noted that the EOS1 does not consider the viscosity aspects of the brine, contrary to the 3D reservoir model presented earlier. Therefore, differences in the pressure behavior of the reservoir between the Table 1 Characteristics of the reservoir model. A visual representation of the mesh can be found in Fig. 2  two models can be attributed to the different EOS. Lastly, the model does not consider any effects inside the wellbores or the surface facilities.

RO-4-SLL RO-3-SLL RO-2-SLL
2.2.1.1. Stratigraphy, geometry, wells and mesh. The model characteristics can be found in Table 1. The stratigraphy represented in the model is an adapted version of the Groningen 3D geological and reservoir model [33]. For the purposes of examining reactive flow through the reservoir, a 2D model is considered sufficient since there is no supporting data for a spatial, 3D variation of chemical characterization. Therefore the 2D model assumes radial symmetry around the wells. A 2D model has the added benefit of reducing the simulation time.
Simplifications had to be made to adapt the model from 3D to 2D. The contacts of the stratigraphic layers in the model are now horizontal, using the average depth of the respective 3D model layers; furthermore, the thickness of the layers is also averaged for each respective interval. The wells are vertical, maintaining a distance of 1275 m from each other at reservoir depth. The horizontal discretization of the mesh is symmetrical; it uses a finer mesh (5 m), close to the wells, that gradually coarsens (40 m) towards the middle part of the reservoir (Fig. 2). The Rotliegend Slochteren layers (upper and lower members) were used as the reservoir body for the injection of cold, and production of hot, water. Vertical discretization ensures a constant thickness of the layers inside the producing interval.

Petrophysical inputs and boundary conditions.
Petrophysical data for porosity and permeability were obtained from Panterra [39] and were scaled with respect to the most proximal well to the project area, as presented in previous research [33]. Probability levels 90%, 50% and 10% are derived per reservoir layer and are denoted by P90, P50 and P10 (summarized in Table 2).
The top and bottom cells were used as boundary conditions for the numerical simulation. The values for the temperature and pressure in the boundary cells were computed based on the respective pressure and temperature gradients (Table 1). Their distance to the production layers ensures that they do not interfere with the reactive flow in the reservoir cells. Furthermore, the boundaries remain valid in this geothermal setting, which is dominated by conductive heat transfer [19].

Chemical characterization.
The mineralogy of the reservoir layers was adapted from the X-ray Diffraction (XRD) and Scanning Electron Microscopy (SEM) analysis of core samples [41] ( Table 3). The TOUGHREACT code requires the specification of secondary minerals that might form as intermediate products of the chemical reactions. The primary and secondary minerals considered and their respective reactions as implemented in TOUGHREACT can be found in Appendix A. The kinetic rate parameters need to be specified for both the primary and secondary mineral phases. They are reported in Appendix B, together with the general rate expression and kinetic rate constant formulations used in the TOUGH-REACT code.
The brine composition is derived from gas field data, as presented in previous geochemical studies of the Permian Rotliegend sandstone [41] in the area of Groningen. Concentration levels presented in Table 4 represent the equilibrated values with the initial rock mineral composition. The Carman-Kozeny porositypermeability relation is used in TOUGHREACT to account for changes in the flow properties due to chemical processes.

Simulation parameters
To account for the uncertainty inherent in such geological data, but also to evaluate operational possibilities, a discrete parameter Table 2 Model layer characteristics. Grouping the permeability values in P90-P50-P10 scenarios, the permeability range is taken into account as a worst-middle-best estimation, based on the petrophysical data presented. Vertical permeability is an order of magnitude lower than horizontal permeability for all layers [40].   Table 4 Brine composition. Data obtained from [41]. The initial reservoir pH is set at 7.
Additionally to the minerals shown here, a thin layer of iron was also introduced to the model in the location of the well cells to account for the well casing. As in other geochemical studies, petroleum exploration and production data are usually the only data available to geothermal studies [18] in conduction dominated geothermal fields. analysis is used, similar to recent studies for geothermal power systems [44]. Uncertainty is taken into account regarding the initial reservoir pressure, the reservoir permeability, flow rate control strategy, re-injection temperature and injection pH ( Table 5). The uncertainty classes cover the initial reservoir state (pressure depletion) and geological uncertainty (mainly permeability related), which are determined by the field conditions and are beyond project control. Additionally, the operational parameters of the flow rate control strategy, the re-injection temperature and the injection pH provide some degree of freedom in designing the geothermal system [33].

Results
For the case study of Groningen, a simulation of the MPC, using appropriately chosen parameters, is first performed and presented. This is followed by the results of the reservoir simulations. The effects of the 243 scenarios on reservoir pressure, permeability change, power output chemical changes at the injector and 2D property changes are analyzed.

Controller and MPC
In order to simulate the district heating network in closed loop with the MPC, the demand pattern for the Groningen case study is firstly introduced. This heat demand pattern is depicted in Fig. 3a and is obtained from a statistical model representing 10,000 households.

Setup
Since the re-injection temperature of the geothermal system has a big influence on the energetic performance of the overall sys- Table 5 Uncertainty classes and respective discrete parameters considered. In total 243 unique reservoir simulation realizations are computed (3 5 ). The flow rate levels are unique for each re-injection temperature level as the heat content of the injected water changes. Therefore MPC, constant min and constant max flow rate levels are adjusted accordingly to reflect this (see also . Model heating energy demand data per hour from the geothermal project of Groningen city (a). An upper limit of 15 MW is used in the MPC for the geothermal system that is sufficient for 77% of the demand points; the demand in excess of this limit is expected to be delivered by other means, namely bio-gas burners. The markers represent the monthly averages of full and geothermal supplied demand. Production levels of the geothermal system resulting from the MPC (b). The markers represent resampled monthly averages of the data series and are used as flow rate levels in the reactive flow simulations. The rate constraints in (19) are set close to zero, resulting in slowly varying production levels of the geothermal system. For this reason, a monthly averaging is a good approximation of the production levels. The increased production in July with respect to June and August can be attributed to an overshoot of the controller in combination with the limited prediction horizon and the various constraints. This behavior could be avoided by considering better predictions, relaxing the constraints and adopting a longer control horizon. Storage levels for the respective scenarios (c). The storage helps to ensure that the demand is met by accounting for the rapid fluctuations of the demand. Moreover, the capacity of the storage is not exceeded in the simulation, resulting in a feasible control signal (i.e. a control signal that satisfies all the constraints at all time).
tem, three different simulations are performed, each with a different return temperature (see Section 2.2.2 Simulation parameters). Along with the different temperatures, several other parameters are introduced in these simulations (see Table 6). The parameters in Table 6 are physical constants or control parameters set by the designer. In order to ensure that the input change is minimized Du max and Du min are chosen close to zero while a is large. It should be noted that the storage size chosen is larger than what would be desired in practice. However, depending on the demand pattern and requirements of the controller, the storage size can be smaller (as discussed in Section 3.1.3). Also the accuracy of the expected demand b d plays a big role in the storage sizing. The demand expected during the simulation is solely based on historical data and therefore does not rely on any predictions. As a consequence, the performance of the controller can be improved if, for example, weather predictions are used.
Let k 0 be the current time step at which the demand is predicted. A heuristic demand prediction is made in which the average demand of the past three days is used, i.e.: for all k 2 ft 1 ; . . . ; minðk 0 þ N c ; k 0 þ 24Þg with k given in hours. For k > k 0 þ 24 (since N c > 24) let the prediction be equal to b dðk c Þ as in (17), where k c is the corresponding hour on the first predicted day. Since no historical data is available on the first two days, it is assumed for simplicity that the average demand of these days are known and used.

Simulation results
The optimization problem in each iteration step is solved using CVX [45]. The results of the simulations are depicted in Fig. 3b and c. The production exhibits a low frequent change in spite of the high frequent change in demand. Furthermore, the volume capacity and production rate constraints are satisfied throughout the simulation.

Influence parameters
In this section, the influence of the parameters on the performance of the controller is discussed. As can be seen from (3) the demand d depends on the reciprocal of DT. For this reason, the demand d increases for a decreasing DT, given that d g remains the same. A decrease in a leads to more fluctuating production and a lower deviation from the reference storage value. Such a change implies a smaller storage device can be used without compromising the feasibility of the solution in one of the iterations.
Adjusting the rate constraints Du min and Du max towards zero has a similar effect as decreasing the value of a. However, Du min and Du max only constrain the extreme rate changes and therefore do not affect small variations in the input. An increase on the control horizon N c can, but does not necessarily, result in a better performance. To illustrate the influence of a, Du min , Du max and N c the same simulation as in Fig. 3 is performed, but using a ¼ 10, Du min ¼ À0:3, Du max ¼ 0:3 and N c ¼ 20 Á 24 as input. The result can be found in Fig. 4 in which a significantly smaller storage size suffices for the total demand. However, taking Du max larger and Du min smaller results in an input that has a higher fluctuation which not desired in the geothermal reservoir. Furthermore, increasing N c leads to an increased simulation time caused by the larger number of states in the optimization. If such a controller is implemented in practice, a necessary condition for the controller to work is that the time to solve the optimization problem (17), (18) should be strictly smaller than the time step, so that it can be implemented in real time. Since the time step is taken as one hour, only large values of N c (e.g. > 10 3 hours and depending on the speed of the device) can result in non-feasible computation times. The minimal and maximal production levels u min and u max are determined by the physical constraints of the pipes, pumps and the characteristics of the well. If the demand pattern exceeds this upper or lower limit for large periods of time, the demand needs to rely on the storage device for the excess demand. Depending on the storage size and state at that time a non-feasible solution might appear and should therefore be treated with caution. Lastly, more accurate predictions result in smaller storage reference value deviations and therefore smaller storages can be used.

Geochemical reservoir simulations
Changes in pressure, temperature, permeability and chemical composition could affect the output of the geothermal system. It is therefore important to investigate the consequences of timevarying production for the reservoir. For these reasons, 2D reservoir simulations are carried out comparing constant and timevarying production. In order to take the uncertainty of the reservoir parameters into account the simulations are performed for several scenarios, as discussed in the methods section. The results shown hereafter include data from all the performed simulations with the exception of the two dimensional plots.
The production rates that result from the MPC simulation are used as input for the reservoir simulation. However, due to the large timespan of the simulation (50 years) and increased number  of simulations (243), a lower resolution is used by extracting monthly averages (Fig. 5). One needs to keep in mind that the resampling might affect the simulation results. That is, due to low frequency resampling, the high frequency behavior of the production flow rate is neglected in the geochemical analysis of the reservoir. However, due to the conservative choices of a, Du max and Du min , the production flow rate has already a low frequency, which implies that the resampling is a good approximation of the original one. The yearly pattern obtained is then used as input for each of the 50 years of the reservoir simulations. In these simulations, the changes in the pressure difference between the producer and injector, the permeability around the injector, the power production and chemical composition in the reservoir are investigated. In the following sections, a comparison between the different flow rate control strategies (MPC, constant max and constant min) with regard to the temperatures in the reservoir is presented.

Pressure difference producer-injector
First the influence of the re-injection temperature, permeability and flow rate control strategy on the pressure is investigated. The pressure difference between the wells is mainly determined by reservoir permeability as seen by the distinct grouping of the results (Fig. 6). A high permeability value (P10) results in a small Dp while for a low permeability (P90) we observe higher levels of pressure difference between the wells. The MPC controlled flow rate levels appear to fall between their respective max and min flow rate intervals. Moreover, the range of Dp is broader for the P90 permeability scenarios; this range is attributed to the other parameters (injection pH and depletion).
A low re-injection temperature (40°C) results in a lower pressure difference compared to a high re-injection temperature. This effect is not very pronounced as the differences are of the order of a few bars. Causally the re-injection temperature effect on the pressure could be explained by the permeability changes in the reservoir, which is presented in the next section (Fig. 8).
All flow rate levels exhibit a slight increase of pressure over time, but this increase is more prominent for the P90 scenarios. Where present, the pressure increase is gradual and never exceeds 50% of the initial pressure. The highest pressure increase is observed in the case of a constant max flow rate and less so for the MPC flow rate levels. The pressure increase for the min flow rate levels is only marginal. These results are consistent with the permeability results presented hereafter.  5. Flow rate levels for the different re-injection temperatures. Since the demand curve is expressed in MW (and it is unique for each re-injection temperature, see Fig. 4b), the flow rates have to be adjusted accordingly when using different re-injection temperatures. Consequently, min and max flow rate control strategy levels for different re-injection temperatures are different for each scenario. Fig. 6. Pressure difference between the wells for all simulations. The data series are colored based on their permeability specifications. The subplot horizontal axis differentiates the flow rate levels (MPC, max, min), while the subplot vertical axis differentiates the re-injection temperature (40, 55 and 70°C). All subplots use the same scale and are therefore cross-comparable.

Doublet power production and energy reserve
As discussed in Section 2.1.1, the power production P (in Watt) of the geothermal system is derived from Eq. (5). However, the temperature and pressure are not necessarily constant which can result in a deviation from the desired power output; this deviation can be attributed to enthalpy changes, according to: where DH is the enthalpy change, DT the temperature change, V the volume, b the coefficient of thermal expansion, T the temperature and Dp the pressure change in the system. The simulations show that the deviations from the desired power output can be attributed to changing pressure levels inside the reservoir (Fig. 7b). This is because these changes ultimately affect the producer T p causing small temperature variations resulting in these minor changes in produced power. The simulation results illustrate the dependency between the flow rate, temperature and power output (Fig. 7). A first observation is that the power output in all simulations that use the MPC flow rate control strategy, is contained within the power production levels corresponding to minimal and maximal flow rate (see also Fig. 5). Secondly, a small decrease in the power output is observed, when moving from no depletion to 200 bar of depletion in the initial reservoir conditions (Fig. 7b). Although small (between 1 MW and 1.7 MW), this decrease needs to be considered in designing and sizing the installation to accommodate for the uncertain levels of pressure depletion. Lastly, the re-injection temperature levels cause minor changes ($0.5 MW) in the power output (Fig. 7c).

Permeability and chemical changes around the injector
The permeability in the reservoir is not necessarily constant; changes in permeability are related to changes in mineral volume fraction, which are in turn affected by the simulation input parameters like the initial mineral compositions and their reaction mechanisms (see Appendix A). Therefore, permeability changes and chemical changes are presented together in this section. Due to the large number of minerals, only the ones with significant changes are discussed. Permeability and chemical changes around the injector well are presented for all the simulations to derive comprehensive insights.
The changes in permeability or mineral volume fraction are, however, not limited to the area around the injector well, but also exhibit spatial differences within the reservoir; the two dimensional changes in permeability and for the minerals that proved significant around the injector are discussed for a sub-set of input parameters in Section 3.2.4. Nonetheless, it is evident that the injector well is the first point at which major physical (e.g. temperature) and chemical (e.g. composition, pH) changes are taking place.
All simulations share an initial permeability increase of about 15% within the first 3-4 years of injection, irrespective of the injection pH and the re-injection temperature (Fig. 8). The temporal volume fraction of Anhydrite closely matches this time interval, where an increase in permeability is observed for all simulations (see Fig. 8). The volume fraction of anhydrite simultaneously reduces to almost zero for all simulations (Fig. 9); it therefore presents a direct link to the permeability increase, since no other mineral volume fraction changes in any of the simulations within the first 3 years. There exists a clear positive correlation between the rate of Anhydrite dissolution and the flow rate control strategy:  higher flow rates lead to faster dissolution. However, these differences only affect the time at which the final volume fraction is reached; the value of the final volume fraction is the same for all flow rate control strategies.
The changes in Iron concentration levels are very low (6th decimal) and appear to be affected by flow rate levels.
Beyond the initial permeability increase, simulation results show that the permeability is mainly affected by the pH value and temperature of the injected fluid (Fig. 8). Notably, following the initial increase of permeability, which occurs in all simulations, the changes caused by an injection pH of 4 and 7 are of the same order of magnitude, but in opposite directions. In both cases, a permeability change of an additional 17-27% takes place. For a pH of 4, the permeability continues to increase, reaching values between 32% and 45% higher than the initial values. For a pH of 5.5, the changes of permeability over time are relatively small; For an injection pH of 7 and after the initial increase, shared by all simulations, a sharp decrease follows. Between production years 20-40, the permeability change decreases. The effect is large enough to reduce the overall permeability to À5% to À15%, with respect to the starting values.
The characteristics of the Dolomite plots, together with its high initial volume fraction, indicate that Dolomite could be the cause of the permeability changes observed in all simulations following the first increase caused by Anhydrite dissolution. Dolomite concentration exhibits a linear change over time, which is caused by the pH level of the injected fluid (Fig. 10). Also, Dolomite is the fourth most abundant mineral in the reservoir (Table 3). Additionally, the volume fractions changes of Dolomite are the highest of all minerals compared to the initial concentration, indicating significant relative changes in its volume.
With an injection pH of 4 Dolomite dissolves, leading to further increases in permeability (in addition to the trends observed during the first 2-3 years, Fig. 8). An injection pH of 5.5 causes small increases or decreases of the initial Dolomite concentration and corresponding effects (i.e. both decreases and increases respectively) to the permeability. A pH of 7 triggers Dolomite precipitation, which is in turn reflected in the permeability decrease; this effect follows the initial increase (attributed to Anhydrite dissolution) observed in all simulations.
Beyond the initial changes common to all simulations and the subsequent differentiation of permeability caused by the injection  pH, a third factor of influence can be distinguished based on the reinjection temperature (Fig. 8). Consistently, lower re-injection temperatures lead to greater increases in permeability than higher re-injection temperatures, for each injection pH.
For an injection pH of 4, the range of increase is influenced by the re-injection temperature; lower injection temperatures lead to higher permeability increases. For a pH of 5.5, the changes in permeability over time are relatively small; following the initial increase of about 15%, minor increments lead up to an overall level of circa 16-17% at a re-injection temperature of 40°C. For a reinjection temperature of 55°C and 70°C, a slight permeability decrease compared to the initial increase of 15% takes place. This reduces the permeability change to an overall increase of 12-13% and 10-11% for the 55°C and 70°C degrees re-injection, respectively. Lastly, for an injection pH of 7 an overall decrease in permeability is observed at the end of the simulation time. These levels are once again causally linked to the re-injection temperature, with higher re-injection temperature (70°C) resulting in a stronger decrease ($15%).
The volume fraction of Albite and Illite is dependent on the re-injection temperature (Fig. 11). Albite dissolves under all scenarios, but changes are more drastic for 70°C and minimal for 40°C of injection temperature. Contrary to Albite, Illite precipitates under all scenarios. However, the changes in Illite volume fraction are very low (5th decimal). An injection temperature of 40°C causes only minor changes, while 70°C more significant ones. The volume fraction changes caused by the re-injection temperature do not seem to be of sufficient magnitude to explain the changes in permeability associated with the reinjection temperature. Moreover, the marginal increase or decrease that is affected by the re-injection temperature for an injection pH of 5.5 would imply dissolution or precipitation affected by the temperature level. A more plausible explanation is the marginal changes in the volume fraction of Dolomite, for an injection pH of 5.5, which are further differentiated by the re-injection temperature (Appendix C).
Minerals with high initial abundance (e.g. Quartz, k-Feldspar and Kaolinite) exhibit little volume fraction change over the 50 years of the simulations. Quartz and K-feldspar volume fractions are mainly affected by the injection temperature (Fig. 11) while the change in Kaolinite volume is caused by the injection pH level (Fig. 10). For all three minerals, however, the changes can be considered minor, in comparison to their initial levels.
Barite and Galena are the only minerals that are not present in the initial state of the reservoir. Barite precipitates under all simulations and the mechanism of precipitation appears to be affected by the re-injection temperature (Fig. 11). Although Barite precipitates in all the scenarios considered, its volume fraction is lower for higher injection temperature. This is in line with observation in other Rotliegend geothermal sites, where lower temperatures favor Barite precipitation [46]. Similarly, Galena also precipitates under all scenarios, although the volume fractions for Galena are the smallest of all minerals (7th decimal). It is also notable that Galena precipitation seems to be independent of temperature, pH level and only small differences can be explained through flow rate levels (Fig. 9).

Two-dimensional property changes
In order to compare the differences between the flow rate levels of the MPC and the max scenarios in terms of lifetime expectancy of the geothermal system, the 2D heat distribution in the reservoir is analyzed. Under the MPC control strategy, and after 50 years of production, the cold front propagates to about a quarter of the distance between the injector and producer wells (Fig. 12a). Using the max flow rate levels for the same re-injection temperature, the cold front propagates further to about a third of the well distance (Fig. 12b). The difference between the two scenarios ( Fig. 12c) highlights the benefits in terms of a more sustainable use of the geothermal resource when utilizing the MPC flow rate levels. Nonetheless, the reservoir can be used for 50 years without a cold front breakthrough under any of the considered scenarios.
Similar to permeability changes, alterations in mineral composition are not limited to the injector but propagate further inside the 2D space of the reservoir model. Due to the large number of simulations, a selection is made to illustrate the changes in the two-dimensional space (see Fig. 13).
The dissolution of Anhydrite in the injector well as the affecting mechanism for the permeability increase is further corroborated by the 2D reservoir plots. A direct relationship can be observed between the spatial distribution of Anhydrite volume fraction change and the changes in permeability (Fig. 13a). A decrease in Anhydrite volume fraction (Fig. 13b) results in a permeability increase (Fig. 13a). This connection also exists between the Dolomite volume fraction (Fig. 13c) and the changes in permeability, but it is less strongly correlated spatially. Changes in Dolomite vol-  ume fraction also appear at a later stage in time (see Video 1) and therefore have a second order effect to the permeability.
Moreover, precipitation (positive volume fraction change) of Anhydrite occurs further away from the injector. This spatially coincides with the changes in permeability, which decreases in this part of the reservoir. It appears that precipitation is stronger at the contact between the upper, less permeable, layers than the middle ones. Similar behavior is observed for Dolomite. Thus, it can be concluded that Anhydrite and Dolomite are dissolved by the water injection and are partly re-precipitated further inside the reservoir. In these simulations, we have not observed a cold-water breakthrough, but in the event that production continues for a longer period we could see the deposition of these two minerals close to the injector. In such circumstances, the permeability around the producer will be negatively impacted.

Discussion
The coupling of the controller to the reservoir model provides new insights on the energy network to which they are both connected. These insights can facilitate a stronger integration between the different parts of the system, if they are analyzed and developed in tandem. Furthermore, possible incompatibilities can also be mapped out. The production of heat becomes constrained by the geothermal system specifications, while the geothermal system itself is also more efficiently utilized as it can respond to changes in demand.
The designed controller is able to generate real-time production levels, using demand predictions and storage level measurements as inputs. It is guaranteed that these productions levels will satisfy several hard constraints. These include storage capacity constraints, minimal and maximal production levels and production rate constraints. The latter constraint is not known at the design phase of a project. Moreover, if such bounds would be available, it is most likely that they would depend on the characteristics of the geothermal system and therefore be case specific. Fortunately, the controller is designed such that updated values can be included once they are available. This implies that the methods and design presented in this paper should be considered as a proof of concept.
It is well known that the stability of a MPC is hard to guarantee [25]. The stability depends on the constraints, demand pattern, objective function, control horizon and demand prediction. A sufficiently large storage device results in a set of constraints that is easier to satisfy and therefore in a stable system. However, due to the high costs involved, this is undesirable in practice. For this reason, additional tuning (see for example [47]), better predictions (see for example [48]), and the integrated control of additional heat sources should be addressed in future work.
With regard to the geothermal system, the relative changes in pressure levels can be explained by the geochemical changes in the reservoir. The changes in pressure over time appear to be minor; this is in line with the mostly increasing permeability around the injector well. Only when using an injection pH of 7 we observe a decrease of the initial permeability levels, and this change only occurs after 25 years of production; thus, such an issue can be identified in time before it affects production. The increase under all other scenarios is attributed to the initial dissolution of Anhydrite.
The geochemical behavior of the reservoir is impacted primarily by flow rates, followed by injection pH and injection temperature. Notably the variable production rates of the geothermal system (resulting from the MPC control strategy) alter only the rate and not the nature of the changes in the chemical reservoir properties. A prominent example is the dissolution of Anhydrite, which appears to be affected by the flow rate. The dissolution or precipitation of other minerals is mainly affected by either re-injection temperature or pH. Nonetheless, the effect of the variable production is beneficial with regards to the utilization of the geothermal resource. This is exemplified by the volume of rock that remains less affected by production using an MPC production rate as opposed to a constant production level (see Fig. 12). The absence of adverse geochemical effects combined with the improved efficiency of heat extraction opens up possibilities for a more sustainable development of geothermal resources [11].
Initial pressure depletion in the Groningen case history does not appear to have a significant effect in any of the results. Especially the geochemical behavior and the precipitation or dissolution of minerals is unaffected by possible pressure depletion. This is unlike the importance of the pressure depletion to the physical aspects of the 3D reservoir in previous research [33]. These differences can be ascribed to the different model setup, namely the inclusion of methane and the use of a different Equation of State in the 3D model compared to the 2D model presented in this work. The 2D model does not include the physical implications of the brine viscosity, as this was a limitation of the EOS1 that was used. This could change the behavior of some of the minerals discussed here.
Due to the nature of the coupling between the two models (i.e. the input for the reservoir simulations is pre-defined for the whole simulation period of 50 years), the effect of varying either the injection pH or the injection temperature during the production lifetime could not be evaluated. This could, however, be of interest for a further study. Several of the minerals discussed are affected by either the injection pH or the injection temperature and altering the injection temperature could be an interesting scenario from the operator perspective. Analyzing a geothermal system with a time-varying production can lead to a better utilization of the reservoir and helps balancing demand and supply in a heat network, resulting in lower CO 2 emissions. The analysis in this paper allows for identifying potential geochemical risks in the development of heat networks which utilize a deep geothermal system. That is, at an early stage of project development and in the absence of production data, this study provides an indication of the temporal behavior with regards to reactive transport at the reservoir level. Moreover, the discussed geochemical implications expand on previously presented risk analysis in literature [33]. The model setup and results could be refined once data from the exploration well become available and the system development progresses. The presented method can assist in the implementation of a demand driven heat network utilizing a geothermal system.

Conclusions
The combined use of an MPC controller and scenario analysis enable a more efficient integration of a geothermal system in a heating or energy network. An MPC is well suited to controlling the production of a geothermal system due to its ability to take several hard constraints into account. The performance of the controller depends greatly on the demand pattern, the tuning of the controller and the quality of the prediction. Furthermore, the performance is also affected by the lower and upper bound for the rate of change of a geothermal reservoir. Since these bounds are usually unknown, obtaining a method to find these bounds is an interesting open problem. The controller is designed such that new values can easily be implemented once they become available. Since a better performance leads to a smaller storage size it is desirable to investigate how this performance can be optimized.
With respect to the energy production, no cold front breakthrough is encountered after 50 years of production under any considered scenario. The findings further suggest that for the case study of the Groningen geothermal project in Rotliegend sandstone, the use of a variable production rate has no adverse geochemical effects on the reservoir. Moreover, it enables a more efficient use of the geothermal resource by limiting the heat extraction to levels dictated by existing demand.
Reservoir geochemical behavior is affected primarily by flow rate levels, followed by injection pH and injection temperature.
The key minerals that affect the injector area are Anhydrite during the first years of production and Dolomite in the following years. Anhydrite dissolution is strongly correlated to an increase in permeability around the injector and is the only considered mineral for which the rate of change is influenced by the flow rate control strategy; higher flow rates lead to faster dissolution but all flow rates eventually cause the same volume fraction to dissolve. Dolomite affects the evolution of the permeability and its change rate is primarily affected by the pH and secondarily by the temperature. An acidic pH favors dissolution and neutral pH precipitation of Dolomite, while lower temperature reduces precipitation and increases dissolution. Consequently Dolomite becomes a crucial mineral for the temporal system behavior.
The geochemical results are representative of a geothermal system in Rotliegened sandstone in the area of Groningen, NE Netherlands. Nonetheless, the controller design, the simulation of the storage device, the simulation of the reservoir model and the analysis that integrates the control engineering and geochemistry domains is readily applicable to other geological contexts and can aid in a more widespread integration of such systems.
with k 25 the rate constant at 25°C, R the gas constant, T the absolute temperature, a the activity of the species and n a constant exponent.
Superscripts or subscripts nu, H and OH indicate neutral, acidic and base mechanisms respectively.
The kinetic rate parameters need to be specified for the both the primary and secondary mineral phases that are constrained to kinetic conditions, in order to evaluate Eqs. (21) and (22). The kinetic data presented are obtained from [50] and are presented in Table 8. Reactive surface area and grain radius data is obtained from [35,41]. Since the reaction rates of anhydrite and calcite are fast compared to the modelling time, they are assumed to react at equilibrium [51]. Additionally, the thermodynamic properties of the chemical species are obtained from the EQ3/6 database [52]. The database data cover a temperature range of 0-300°C, a pressure range up to several hundreds of bars and is applicable for a brine concentration equivalent of up to 6 molal. Thermodynamic data for the minerals anglesite, barite and lead were further added, obtained from the Thermoddem database [53].

Appendix C
Dolomite dissolution/precipitation is primarily affected by the pH (Fig. 10) by secondarily also by the injection temperature (Fig. 14). Higher injection temperature leads to more precipitation and less dissolution. Table 8 Relevant parameters for the evaluation of the kinetic rate laws according to Eqs. (21) and (22). The rate constant k25 is in mol=m 2 , the activation energy Ea in kJ=mol, the reactive surface A in cm 2 =g. Where the acidic and/or base mechanisms are not present the corresponding reaction mechanisms are not considered. Data are obtained from [50].