Geochemical Changes Associated with High-Temperature Heat Storage at Intermediate Depth: Thermodynamic Equilibrium Models for the DeepStor Site in the Upper Rhine Graben, Germany

: The campus of the Karlsruhe Institute of Technology (KIT) contains several waste heat streams. In an effort to reduce greenhouse gas emissions by optimizing thermal power consumption on the campus, researchers at the KIT are proposing a ‘DeepStor’ project, which will sequester waste heat from these streams in an underground reservoir during the summer months, when the heat is not required. The stored heat will then be reproduced in the winter, when the campus’s thermal power demand is much higher. This paper contains a preliminary geochemical risk assessment for the operation of this subsurface, seasonal geothermal energy storage system. We used equilibrium thermodynamics to determine the potential phases and extent of mineral scale formation in the plant’s surface infrastructure, and to identify possible precipitation, dissolution, and ion exchange reactions that may lead to formation damage in the reservoir. The reservoir in question is the Meletta Beds of the Upper Rhein Graben’s Froidefontaine Formation. We modeled scale- and formation damage-causing reactions during six months of injecting 140 ◦ C ﬂuid into the reservoir during the summer thermal storage season and six months of injecting 80 ◦ C ﬂuid during the winter thermal consumption season. Overall, we ran the models for 5 years. Anhydrite and calcite are expected mineral scales during the thermal storage season (summer). Quartz is the predicted scale-forming mineral during the thermal consumption period (winter). Within ~20 m of the wellbores, magnesium and iron are leached from biotite; calcium and magnesium are leached from dolomite; and sodium, aluminum, and silica are leached from albite. These reactions lead to a net increase in both porosity and permeability in the wellbore adjacent region. At a distance of ~20–75 m from the wellbores, the leached ions recombine with the reservoir rocks to form a variety of clays, i.e., saponite, minnesotaite, and daphnite. These alteration products lead to a net loss in porosity and permeability in this zone. After each thermal storage and production cycle, the reservoir shows a net retention of heat, suggesting that the operation of the proposed DeepStor project could successfully store heat, if the geochemical risks described in this paper can managed.


Research Overview
A climate-neutral society requires energy transition concepts with increasing shares of renewable energies. Currently, a large percentage of renewable power is derived from fluctuating sources that need large storage capacities. Chemical storage using batteries is limited, but renewable gas storage or power-to-heat-to-power concepts offer large capacities in the Earth's subsurface. In 2012, 24% of the European gas storage volume, i.e., 23 billion m 3 , was found in Germany. Of this, 47% was accommodated in hydrocarbon reservoirs and an additional 2% in aquifers [1]. Coupling renewable energy sectors and balancing seasonal heat and power demand fluxes, in turn, require new concepts in subsurface energy storage. In the case of heat, this involves increasing subsurface temperatures.
High-temperature heat storage in the subsurface is feasible in abandoned hydrocarbon reservoirs [2]. These reservoirs have an added economic advantage in that their geology, geo-and petrophysical properties, depths, and geometries have already been extensively characterized. Sticker et al. [2] concluded that in the Tertiary hydrocarbon reservoirs of the Upper Rhine Graben (URG), a volume of~155,000 m 3 of water at 140 • C can be injected over six months in a doublet system, depending on the reservoir's permeability. This equates to a storage potential of about 12 thermal gigawatt hours GWh th .
At present, thermal storage technologies range from hot water tanks and gravel pits over borehole heat exchangers to aquifer thermal energy storage (ATES) systems. In Europe, there are numerous gravel pits and ATES systems in Denmark and the Netherlands, respectively. With more than 2800 systems in operation, the number of ATES systems is continuously increasing. These systems provide more than 2.5 thermal terawatt hours. TWh th for heating and cooling purposes, saving several thousand tons of CO 2 emissions per year [3]. They generally operate at temperatures below 50 • C.
To overcome current temperature limitations, the Karlsruhe Institute of Technology (KIT) has designed the DeepStor geo-scientific infrastructure for hydraulic, thermal, and chemical experiments at temperatures up to 170 • C. These experiments will be carried out in the water-saturated rim of the former Leopoldshafen oil field. The design, as a scalable multi-storage system, allows the infrastructure to be extended and integrated into the KIT campus's district heating network with an inlet temperature of 110 • C. Renewable heat sources are located in a deep geothermal system underlying the campus [4]. In the nearby Leopoldshafen-20 well, temperatures of about 170 • C were measured at 3000 m depth. Additionally, the Bioliq plant [5] and the low-temperature power plant MoNiKa [6] are available as heat sources on the KIT campus.
Dissolution and precipitation of minerals in the reservoir and surface plant, as well as reactivation of clay particles in the reservoir, are among the most critical challenges for the installation and operation of DeepStor. Of the 2800 ATES systems installed worldwide, only 21 of them are high-temperature (HT-ATES) systems. Only one of these is running (in 2020), six are planned to be constructed in the near future, and the rest are abandoned [7]. A risk assessment shows two dominant technical risks for HT-ATES systems. Exploration risk and scaling and clogging are rated as 'medium' risks but can be even higher for deeper reservoirs (>1000 m). For most of the abandoned HT-ATES locations, scaling and clogging were crucial risks, with relatively high severities and occurrence probabilities. To construct and run a productive storage system, an analysis and prognosis of the changing geochemistry and its effects on the porosity and permeability of the reservoir are needed. Seasonal fluctuations in production that differentiate HT-ATES systems from deep geothermal production sites and hydrochemical kinetics in the reservoirs are almost unexplored. Clogging was also identified in an HT-ATES system at TU Delft due to the precipitation of carbonates while injecting hot water into the reservoir [8]. Investigations on HT-ATES systems in the Netherlands identified clogging, especially at high temperatures, as a crucial problem [9].
The purpose of this paper, therefore, is to provide a preliminary assessment of the geochemical risks associated with the operation of DeepStor. Specific risks investigated are mineral scaling within the surface installations and formation damage within the reservoir. In this particular system, scaling is predominantly induced by the cooling and reheating of geothermal fluids as they are circulated through the plant's plumbing. Formation damage is caused by complex, non-isothermal fluid-rock interactions in the reservoir. We used equilibrium thermodynamics to identify the quantity and species of potential scale-forming minerals at various stages of the plant's operation. We then used a 2D reactive transport model to predict the long-term geochemical, hydrogeological, and thermodynamic evolution of the reservoir. In addition to presenting a comprehensive method of preliminary geochemical risk assessment in geothermal systems, the results from this study will be used to plan the brine filtration and treatment infrastructure and processes for DeepStor, and for predicting the regional-scale geochemical risks associated with this type of geothermal storage concept.

Study Area's Geologic and Thermal Setting
The study area is located in the central segment of the NNE-trending, essentially nonvolcanic Upper Rhine Graben (URG, Figure 1) that represents a major accumulation zone of Cenozoic sediments in the URG [10]. The URG developed from about 47 Ma onwards between the Schwarzwald, Kraichgau, Odenwald, Vosges, and Pfälzerwald escarpments, with maximum vertical displacement along the Western and Eastern Main Border Faults (WMBF and EMBF). Along the latter, two geothermal projects, Riehen (CH) and Bruchsal (GER), are in operation. The graben-filling Cenozoic sediments reach up to 3500 m in thickness [11] and are subdivided into Pliocene to Holocene fluviatile siliciclastics that unconformably overlie late Eocene (Lutetian) to early Miocene (Burdigalian) syn-rift successions. These successions were deposited in a generally low-energy environment under changing marine-brackish-fluvio-lacustrine conditions that contain depositions of carbonates, sands, silt, clays, marls, and locally intercalated evaporites (e.g, [11,12]). Sediment transport occurred at both low and high angles to the rift axis. Early Eocene and Oligocene sediments prograded generally northward from depocenters in the southern and central segments. Early Miocene (Aquitanian, Burdigalian) sediments show pronounced depocenters in the northern segment, with decreasing thicknesses towards the south. These segments are cut by a regional unconformity that developed between about 5 and 18 Ma [12].
Graben-filling sediments have been explored and exploited for hydrocarbons [13]. The related reservoirs such as the Leopoldshafen field are now being considered for hightemperature thermal energy storage. This former oil field is located in the footwall of the Leopoldshafen fault ( Figure 2). This major NNE-striking normal fault [14] is part of the EMBF and displays syn-depositional normal faulting from the Oligocene to the early Miocene, with footwall syn-rift strata dipping shallowly to the NE [15]. Cumulatively, about 188,000 t of oil was produced from a~5 m-thick fine-grained mica-bearing calcareous sandstone reservoir layer embedded in the up to 100 m-thick, generally marly Meletta Beds [13,14]. The Meletta Beds are part of the Froidefontaine Formation and were deposited under marine conditions during a regressive phase. They display coarsening upward trends in the late Rupelian (about 29 Ma), followed by a phase of transgression in which fining upward trends in the upper part of the Meletta Beds grade into the marly, still marine to brackish Cyrena Beds. These beds record another phase of regression and further development into the fluvio-lacustrine Niederrödern Formation [16][17][18]. Cross-section west of the KIT Campus North area across the former Leopoldshafen oil field (modified, own elaboration based on data from [14]). The first-order approximation of the temperature distribution with depth is based on a conductive heat transport with a temperature gradient of about 45 • C/km. Dotted layers mark the Meletta and Cyrena Beds in the Froidefontaine and Niederrödern Formations, respectively.
The Leopoldshafen field was assessed through about 20 boreholes down to~3000 m depth. The temperature distribution in the Leopoldshafen field is estimated in a firstorder approximation with corrected temperature data from Sauer et al. [19]. Corrected temperature data revealing about 100 • C at 2000 m depth include measurements in the sandstone layers of the Meletta Beds in the Leopoldshafen-1a and -15 wells.

DeepStor Experimental Design
The DeepStor scientific infrastructure comprises a test and a monitoring well that will be drilled to a depth of about 1400 m depth through the Meletta sandstone layers. The test well explores several possible target horizons, including the Meletta and Cyrena Beds in the Froidefontaine and Niederrödern Formations. This well will be completed to sustain geothermal brine production at 170 • C. The infrastructure is designed for experiments on transient hydraulic, thermal, and chemical loading and unloading of the reservoir. For this purpose, after hydrocarbon separation from the fluid, the residual fluid will be stored in a 4000 m 3 basin. A bypass with a heat exchanger is installed to carry out tests at elevated temperatures, initially up to 140 • C. This scientific design can be extended to include a scalable multi-storage system that connects to the heat supply of the KIT Campus North ( Figure 3). In this stage, renewable heat from the deeper Mesozoic layers, or other sources such as waste heat from the Bioliq plant, may be stored. Based on results from the Leopoldshafen-20 well, the expected temperature in the deep Mesozoic layers is~170 • C at 3000 m depth. The Bioliq plant provides a range of waste heat temperatures between 100 and 800 • C at the H 2 O/CO 2 separator and the catalytic reactor, respectively. In connection with the existing area-wide local heating network with an inlet temperature of 110 • C, KIT Campus North offers optimal conditions for the extraction, seasonal storage, and distribution of heat from deep geothermal energy. In this study, we carried out reactive transport models on the future circulation at the DeepStor site through a hot and a cold well. In a first approximation, we used a step function with flowrates of −10 and 10 kg/s over six months. During the thermal energy storage period (summer period), fluid was injected through the hot well at a temperature of 140 • C, while fluid at the reservoir temperature was produced from the 'cold well'. In the following recovery period (winter period), the fluid was produced at an enhanced 'reservoir' temperature through the hot well, while 80 • C hot water was injected through the cold well into the reservoir.

Reservoir Parameters and Boundary Conditions
In order to assess the DeepStor reservoir rock, the mineralogy of the reservoir rock was determined using quantitative X-ray diffraction (D8 Discover, Bruker Scientific Instruments; Billerica, MA, USA) on analogue samples from the lower Meletta Beds. Samples were taken from the clay pit 'Dammstücker' about 40 km NE of the DeepStor site near Nussloch (Germany). The sandstones of the Meletta layers occur here in a rift-edge clod and are unconformably overlain by Quaternary fluviatile sands and gravels. They are weakly consolidated, dark gray, micaceous, carbonatic, fine sandstones with a minimum thickness of 2 m and are located in the extreme southwest corner at the lowest level of the clay pit (UTM Zone 32: 476,622 m E; 5,462,029 m N). The results are presented in (Table 1) and reveal predominately quartz (28%), calcite (27%), dolomite (15%), and albite (12%). In addition, a mixture of orthoclase (6.5%), muscovite (4.8%), biotite (3.9%), clinochore (2.5%), and iron oxides (1.2%) is observed. These values were normalized to 20% porosity (see below) and used as reactants in the model. A representative Meletta Beds' brine chemistry of the Leopoldshafen field (Table 2) was used as the initial fluid concentration in the equilibrium reaction models. This fluid was taken from a 1958 oil field brine measurement conducted by C. Deilmann Bergbau GmbH and supplied to the authors (reproduced with permission) by Neptune Energie Deutschland GmbH. The Meletta sandstones in the Leopoldshafen field are separated typically in layers with thicknesses reaching 8.1-9.5 m. Porosities reach 5.8-29.1% [19] but vary in the different wells. In the wells Leopoldshafen−7 and −8, exploitable 15 to 25 vol.% pore spaces with permeabilities of 29-64 mD are observed at 1223-1224 m and 1240-1245 m depths, respectively. In the model, we assumed an initial porosity of 20%. The permeability was derived from the Geochemist's Workbench's (20) empiric relationship between porosity and permeability, Equation (1), that represents a minimum value and thus a conservative approach: where (k x ) is the permeability in the x-direction (m 2 ), the porosity is given as a volume fraction, and A and B are user-defined variables that, in our case, were initially set to 15 and −17, respectively. The coefficients were chosen to make the initial values align. With a 20% porosity, the permeability was calculated as 1 × 10 −14 m 2 , which is at the low end of the range provided by Sauer et al. [19]. The permeability anisotropy (kx/ky) was assumed to be 1, i.e., the permeability is isotropic. This function allows the permeability of the reservoir to change as the porosity changes with precipitation and dissolution during the simulation. Further rock properties were assigned according to the NIST Standard Reference Database (Number 69) [20] and are provided in Table 3. The boundary conditions in the reservoir were set as open flow boundaries. The overall size of the modeled section was large enough such that the boundaries are far enough away from the wells in order to not be affected by the flow regime. Nonetheless, the effects of the open flow boundaries can be seen at the edges of the modeled flow result figures in Section 3 below.

Reactive Transport Modeling Procedure
In order to design DeepStor, a preliminary assessment of the geochemical risks associated with cooling and reheating geothermal fluids was carried out using equilibrium thermodynamics and 2D reactive transport modeling. A combination of the 'React' and the 'X2T' applications in the Geochemist's Workbench ® (GWB, release 14; Aqueous Solutions LLC; Champaign, IL, USA) program was used to identify the potential for mineral scale precipitation and formation damage in the reservoir. In both applications, we used the Lawrence Livermore National Lab's thermo database (thermo.tdat) which calculates ion activity products using the B-dot equation. The B-dot equation is a modified Debye-Hückel equation, which takes into account charged species' interactions with the surrounding water, as well as with other charged species in the solution. Please see the Discussion section (Section 4) for an explanation as to why this ion activity model was used, as opposed to the Pitzer formulas.
The 'React' application was used to identify potential scale-forming minerals in the proposed system's surface infrastructure during seasonal heating and cooling of the produced brine. The time-resolved potential for formation damage was evaluated using the 'X2T' reactive transport application. The GWB's X2T program also employs a finite difference grid with both spatial and temporal resolutions. Reactive transport was investigated over five years, i.e., five cycles with six months each of injection through the hot and cold wells. Data retrieval was carried out at time intervals of two months, summing up to 30 retrieval steps in total.
The Meletta Beds' reservoir was conceptualized as a 600 × 600 × 10 m (x,y,z) box divided into a 60 × 60 matrix of 10 × 10 × 10 m grid cells ( Figure 4). The hot and cold wells were placed at the coordinates 205 m (x) and 305 m (y), and 405 m (x) and 305 m (y), respectively, in the center of the cells, rather than at a node. The parameters and boundary conditions of the models are detailed in Tables 1-3. First, the equilibrium state and charge balance error of the in situ fluid at the initial reservoir temperature of 80 • C were calculated by speciating the measured fluid (Table 2) in the Geochemist's Workbench's React program [21]. Concurrently, these initial conditions of the reactive transport model were set up in the 'X2T' application. As the GWB only considers chemical components and not the specific mineral stoichiometry, the volume fractions of minerals bearing similar components were combined to form a simplified reservoir mineralogy while maintaining the correct balance of chemical components. We used the 20% normalized mineralogy (Table 1) as the initial reference. Amorphous silica, calcite, albite, muscovite + orthoclase, dolomite + clinochore, and biotite (as annite) + iron oxide were entered into the model as mineral species, accounting for silica (SiO 2(aq) ), bicarbonate (HCO 3− ), aluminum (Al 3+ ), potassium (K + ), magnesium (Mg 2+ ), and ferrous iron (Fe 2+ ), respectively, as system components. The remaining chemical components of the system, i.e., sodium (Na + ), calcium (Ca 2+ ), chloride (Cl − ), and sulfate (SO 4 2− ), were entered in the 'initial' pane as aqueous species.
The modeling procedure is illustrated in Figure 5. The measured fluid ( Table 2) was equilibrated and charge balanced at 80 • C to attain the chemistry of the in situ fluid. In cases where dolomite was found to be oversaturated in the fluid phase, the model was re-run with dolomite precipitation suppressed. This adjustment was conducted because kinetic limitations are known to prevent dolomite from precipitating from an aqueous solution. Due to the complexity of the model, the convergence criterion was set to ε = 4 × 10 −9 . To mimic the effects of production and thermal storage, this equilibrated fluid was heated and re-equilibrated at 140 • C. Oversaturated minerals were removed from this re-equilibrated fluid, and the fluid was then injected through the hot well in the reactive transport model at 140 • C and 10 L/s for a six-month period, corresponding to the summertime thermal power storage season. The fluid composition at the cold well after the six-month thermal storage season was then cooled to 80 • C and chemically re-equilibrated. This cooled fluid, with the oversaturated phases removed, was then reinjected into the reservoir through the cold well at 80 • C and 10 L/s for six months, corresponding to the wintertime thermal power consumption season. Then, the composition of the fluid at the cold well was produced, heated to 140 • C, re-equilibrated, and injected into the hot well for a second thermal storage period. This process was repeated cyclically for 5 years of plant operation.
Potential scale-forming mineral species and the amount of precipitation required to return the fluid to equilibrium were determined between each seasonal cycle. Potential scale-forming minerals are those that are predicted to precipitate during heating and cooling of the circulating fluid once it is removed from the reservoir rock, i.e., in the wells and surface infrastructure. Formation damage was assessed within the X2T program, by analyzing precipitation/dissolution reactions in the reservoir and quantifying their effects on the reservoir's porosity and permeability.

Results
The simulation results are presented in three separate sections. In the first section, we describe mineral scaling risks in the system's infrastructure, i.e., wells, pipes, and heat exchangers. We identify potential scale-forming mineral species and the amount of their predicted precipitation, if the produced waters return to equilibrium before reinjection. This first section also includes a presentation of the hot and cool fluids' compositions at the start of each new season. These fluids serve as the seasonal inputs in the reactive transport models. Next, we map the seasonally changing thermal plume. The third section of the results describes the predicted risks of formation damage in the reservoir. We look at seasonal precipitation/dissolution and cation exchange reactions in the reservoir. Finally, we show the effects of these reactions on the porosity and permeability of the reservoir.

Equilibrium State of the Initial System
At 80 • C, the initial fluid is supersaturated with respect to calcite by~52 mg of calcite per kg of fluid (molal; mg/kg). Precipitating this calcite leads to a pH reduction of~0.44. The initial measured solution has a charge imbalance of about −1.6 meq/kg, or about 0.04%.
After charge balancing the initial (measured) reservoir fluid and removing the excess calcite, the resultant fluid is equilibrated with the reservoir mineralogy ( Table 1). The most significant result of this equilibration of the measured fluid with the reservoir mineralogy is the uptake of nearly 0.01 mg/kg of Al 3+. This species was absent from the original fluid measurement. The uptake of Al 3+ is accompanied by small decreases in the Fe 2+ , HCO 3 − , and SiO 2(aq) concentrations. All of the other dissolved ionic species' concentrations (Na + , K + , Mg 2+ , Ca 2+ , Cl − , and SO 4 2− ) are sub-conservative throughout the equilibration reaction between the measured fluid, charged-balanced initial fluid, and reservoir minerology.

Equilibrium Models over Time
The thermal storage and recovery cycles are described in the Methods section. During the first thermal energy storage cycle, the measured fluid, now equilibrated with the reservoir mineralogy, is produced in the cold well, heated to 140 • C, and re-equilibrated. The re-equilibrated fluid is then injected into the hot well. The compositions of the reservoir fluid and the fluid re-equilibrated at 140 • C are shown in Table 4 below. Table 4. Difference in fluid phase chemistry (mg/kg) between the initial input composition and the composition after equilibration and charge balancing. Calcite becomes increasingly oversaturated as the temperature increases. At 140 • C, 100 mg/kg of calcite must precipitate to bring the solution into equilibrium with respect to calcite. Supersaturation of the solution with respect to anhydrite starts above~130 • C and reaches~119 mg of anhydrite per kg of fluid (mg/kg) at 140 • C. The fluid is also supersaturated with respect to antigorite at 140 • C, but due its high, i.e., ≥3, reaction order, it is unlikely to precipitate [22].
To identify potential scaling risks, we culled data from the models as follows: Potential scaling minerals were identified every six months, after heating or cooling on the surface. Heated fluids were injected in the hot well, and cooled fluids were injected in the cold well. On the hot side, albite, amorphous silica, antigorite, calcite, mordenite (K), saponite (Ca and Na), and talc were all identified as potential scale-forming minerals. On the cold side, albite, amorphous silica, calcite, muscovite, saponite (Ca and Na), and talc were all identified as potential scale-forming minerals. As mentioned early, due to the higher reaction orders (≥3), it is unlikely that clays or micas will precipitate from the solution (21). Thus, the major scaling concerns discovered by this study are amorphous silica, anhydrite, and calcite. The extent of supersaturation of the fluid with respect to these minerals for each time step is summarized in Table 5.  Figure 6 shows the evolution of the thermal plume over the initial thermal storage and production seasons, i.e., months 0-12. The initial reservoir temperature is 80 • C. Upon injection of 140 • C pre-heated fluid during the thermal storage period (months 0-6), a thermal plume appears around the injection well, i.e., the hot well. This plume is subcircular, with a tail towards the production well (cold well) appearing in month 4. The plume reaches a maximum radius of about 100 m after month 6, with the tail extending another~30 m in the direction of the production well. After month 6, the cycle is reversed. The hot well becomes the producer, and the cold well becomes the injector. The plume is reproduced, leaving only a diamond-shaped area of residual heat in the reservoir within ã 50 m radius of the hot well. The temperature in this region increased from the ambient 80 to about 100 • C.  Figure 7 shows the evolution of the plume for the remainder of the simulation, i.e., through month 60. The cycle described in the previous paragraph repeats on a yearly basis, with a thermal plume out to about 100 m from the hot well forming during the thermal storage period, and the plume being retracted again through the hot well during the heating season. The size of the main plume does not change significantly during the storage period, although the tail towards the cold well does become thicker and marginally longer. The diamond of residual heat around the hot well at the end of the cycle increases on a yearly basis. The interior of the warmed area is heated to about 120 • C out to a distance approaching 75 m from the wellbore. This area is surrounded by a 10-20 m boundary area of temperatures above~100 • C.

Overview
In order to identify the potential for formation damage in the reservoir, reactive transport phenomena were investigated over a period of five years. Anhydrite, daphnite, dolomite, minnesotaite, muscovite, quartz, saponite-Ca, saponite-Na, and siderite were identified as reservoir mineral alteration products. In the following section, we show the temporal and spatial occurrence of these minerals over the course of the five-year simulation. For ease of viewing, we have placed these minerals into two groups. Group 1, shown in Figure 8, contains all of the non-phyllosilicate minerals. Group 2, shown in Figure 9, contains the phyllosilicates.

Non-Phyllosilicates
The anhydrite presence spikes at the hot well from 0 to 35 kg per m 3 of reservoir rock (kg/m 3 ) during the initial thermal storage cycle (months 0-6). In the subsequent period, the concentration reduces to about 0.25 kg/m 3 at the hot wellbore. After 12 months, there is no more formation of anhydrite in the system. Throughout the entire five-year period, there is no anhydrite formation in the region around the cold well.
Initially, there is 450 kg/m 3 of dolomite in the reservoir. Throughout the simulation, this value increases to 1000 kg/m 3 at the hot wellbore by month 60. Immediately adjacent to the hot well, the concentration decreases to 400 kg/m 3 . The dolomite mass decreases slightly in the region adjacent to the cold well throughout the duration of the simulation.
Calcite displays an opposite behavior to dolomite. Whereas the dolomite concentration spikes at the hot wellbore and is depleted adjacent to the cold wellbore, calcite is strongly depleted at the hot wellbore and enriched above its starting mass adjacent to the cold wellbore. The calcite presence starts at 580 kg/m 3 in the reservoir. After 6 months of injection through the hot well, the concentration at the wellbore is 0 kg/m 3 . This value increases to a maximum of 650 kg/m 3 20 m from either side before reducing to initial values about 100 m from the wellbore. The cold well is unaffected during this initial injection. At the end of month 12, the concentration of calcite at the cold well increases to about 650 kg/m 3 20 m from either side of the well before reducing to initial values about 60 m from the well. During the same period, the maximum value of 660 kg/m 3 adjacent to the hot well is reduced to about 625 kg/m 3 , with the concentration at the wellbore itself being unaffected; calcite is absent at the wellbore itself.
The initial amount of quartz is 595 kg/m 3 . At first, the amount of quartz adjacent to and within the hot well increases to~700 kg. Initially, at the cold well, quartz remains at the initial value. After month 12, the quartz concentration around the hot well decreases for the remaining duration of the simulation. Concurrently, the quartz concentration at the cold well increases for the remainder of the situation, upwards of 700 kg/m 3 by month 60.
Siderite is absent from the initial reservoir assemblage. It remains absent until a concentration of 10 kg/m 3 appears at the cold well at the end of the first heating period (month 12). After this spike, siderite is removed from the system.
The initial concentration of albite in the reservoir is~240 kg/m 3 . At the beginning of the simulation, albite becomes strongly depleted within 50 m of both wellbores. Over time, the albite concentration at the hot wellbore rebounds. By the end of the simulation, it reaches >600 kg/m 3 at the hot wellbore itself and remains elevated within~25 m of the wellbore. Beyond this, albite remains depleted in relation to its initial concentration to a distance of~50 m from the hot wellbore. After its initial depletion, the albite concentration remains close to zero within 50 m of the cold wellbore for the duration of the simulation.

Phyllosilicates
The initial amount of muscovite is 210 kg/m 3 . After six months, the amount of muscovite reaches >300 kg/m 3 at the hot wellbore and~290 kg/m 3 within a 50 m radius from the hot well. Subsequently, the amount of muscovite in the hot wellbore itself decreases to~0 kg/m 3 over the remainder of the simulation. The muscovite concentration remains elevated over the initial values within 50 m of the hot well. Muscovite at the cold well increases from its initial mass to~350 kg/m 3 within the first injection cycle, and the radius of this additional muscovite mass also increases to~20 m around the cold wellbore throughout the duration of the simulation.
During the initial injection at the hot well, annite (initial concentration of 130 kg/m 3 ) reduces to zero between x = 180 and 250 m and increases to~135 kg/m 3 between x = 250 and 320 m. After month 12, the annite concentration between the wells in the section x = 250-350 m gradually increases to 140 kg/m 3 . During the first injection at the cold well, the concentration of annite spikes a little bit in the space x = 380-420 m before decreasing to 0 kg/m 3 within the wellbore over the course of the simulation. This depletion extends to a radius of~20 m around the wellbore.
Daphnite initially does not occur in the reservoir. At the end of the first injection cycle, the daphnite mass spikes dramatically to >120 kg/m 3 adjacent to the hot wellbore. The value decreases over time back to the initial mass (0 kg/m 3 ) at the hot well. There is no appearance or disappearance of daphnite at or adjacent to the cold well throughout the simulation.
The starting concentration of minnesotaite is 5 kg/m 3 . After month 6, this concentration decreases to 0 kg/m 3 within the hot well but spikes to~140 kg/m 3 within 50 m of the hot well. A similar phenomenon is observed at the cold well, with the exception that minnesotaite is also present within the well. The growth of minnesotaite extends to a radius of~20 m around the cold well. Between the wells and outside of these respective radii, the minnesotaite concentration is reduced to 0 kg/m 3 .
Saponite-Ca and -Na have initial concentrations of 0 and 1 kg/m 3 , respectively. Both species spike to~350 kg/m 3 in the hot wells, with less pronounced increases (up tõ 50 kg/m 3 ) within 50 m radii of the hot well. These increases occur first with the calcium saponite species and later on with the sodium saponite species. There is small spike (up tõ 10 kg/m 3 ) of sodium saponite at the cold well towards the end of the simulation. There is no change in the calcium saponite mass at the cold well throughout the simulation.

Effects of Mineral Alteration on Reservoir Hydraulic Properties
The combined effect the mineral alteration in the reservoir is a marked increase in porosity and permeability immediately adjacent to the wellbores, followed by a decrease in these hydrogeologic parameters distally from the wellbores. The changes in porosity and permeability over time in the x-direction of the simulation are shown in Figure 10.
The initial reservoir porosity is 0.20. At the hot wellbore, a maximum porosity of 0.28 is attained at the end of the simulation, i.e., month 60. The spike in porosity around the hot well extends to a radius of about 10 m from the wellbore. Outside of this radius, the porosity decreases slightly below the starting value to a minimum of~0. 19. The decreased porosity zone continues out to a radius of about 75 m from the wellbore but is most pronounced between about 10 and 30 m from the wellbore. A similar pattern of porosity reduction is seen around the cold well, too. Here, the porosity is reduced in the wellbore to~0. 18, which extends about 50 m to either side of the cold well.
The reservoir's changes in permeability mirror the change in porosity. The initial reservoir permeability is 0.15 × 10 −9 cm 2 . The permeability spikes at the hot well, where it attains a maximum value of~1.6 × 10 −9 cm 2 after 60 months. This spike tapers off withiñ 10 m of the hot well. At the cold well, the permeability decreases only slightly withiñ 10 m of the wellbore. Beyond~10 m from either wellbore, no change in the reservoir permeability is observed throughout the simulation.  We have presented a thermo-hydraulic-chemical model of scaling and formation damage risks associated with the operation of a seasonal geothermal storage system in Germany's Upper Rhine Graben. Due to the complexity of the model, several factors were considered to find a balance between the model's accuracy and the time and computing power required to calculate it. These factors include the selection of an ion activity model; consideration of the system's pressure; and consideration of mineral precipitation, dissolution, and alteration kinetics.

Ion Activity Models
We used the Geochemist's Workbench's Lawrence Livermore National Laboratory thermodynamic database (thermo.tdat). This database employs a B-dot ion activity product model, which is a modified version of the Debye-Hückel database. We used this database as opposed to either of the two Pitzer databases (thermo_hmw.tdat; thermo_phrqpitz.tdat) included with the Geochemist's Workbench. Whereas the Debye-Hückel ion activity models consider electrostatic interactions between dissolved ions and the polarity of water, the Pitzer model considers interactions between ions of like and unlike charges, as well as interactions between charged ions and water. The two models will speciate aqueous complexes differently, thereby potentially leading to discrepancies in the calculation of mineral saturation states. These discrepancies can lead to either false positive identification of potential scale-and formation damage-forming minerals, by incorrectly identifying the fluid as supersaturated with respect to these minerals, or failure to recognize potential scaleand formation damage-forming minerals due to under-calculation of a fluid's saturation state with respect to these minerals. Additionally, varying calculations of saturation indices can lead to improper determinations of the extent of mineral precipitation required to bring the solution to equilibrium.
The brine investigated here has a starting salinity of 119 g/kg total dissolved solids, which is the equivalent of a nearly 2 molar sodium chloride solution. This concentration is well above the threshold for employing the Pitzer formulations at atmospheric conditions. We chose to use the Debye-Hückel model for two reasons. First, all commercially available Pitzer databases have limited chemical species as inputs. Typically, they contain only seawater species and are thus missing some key inputs for our study. Of particular interest in this study is the role of aluminum in mineral alteration in this study. The presence of aluminous silicates in the reservoir, as well as the potential for clay mineral alteration, requires selecting a database that contains trivalent aluminum (Al 3+ ). Furthermore, Pitzer parameters are typically only calculated at standard state conditions. The Geochemist's Workbench explicitly warns against their use at elevated temperatures. Thus, Pitzer databases currently available are not suitable for geothermal energy research. Applying the Pitzer model to our simulation would require experimental parameterization of additional ions at geothermal-specific conditions. This fact alone necessitates the use of the B-dot (Debye-Hückel) ion activity model in this study. Nevertheless, the Debye-Hückel model can be effectively used to identify potential mineral dissolution, precipitation, and alteration products, and to predict at an order-of-magnitude precision the total change in mineral mass. The study presented here effectively identified these minerals on a preliminary basis. More sophisticated models can be used in future studies to improve the accuracy and reliability of the results.

Pressure Effects
We chose to ignore the effects of pressure on mineral solubility in this study for two reasons. First, the equilibrium constant for each identified mineral would need to be manually calculated and entered into the model at each time step. Such action would dramatically increase the time required to execute the simulations. Furthermore, a plethora of studies [23][24][25][26] show that, in comparison with temperature, pressure plays a relatively small role in determining mineral solubilities at lower pressures. Many potential scaling minerals' solubilities are not substantially affected by pressure changes until the system's pressure is in the range of a few tens of megapascals. Given that the proposed DeepStor project is operating at depths less than two kilometers, it is unlikely that considering pressure will significantly alter the results of this preliminary study. Pressure effects may influence the total mass of precipitated minerals, or the specific species of clay mineral alteration, but they are not expected to influence the overall processes involved in mineral precipitation, dissolution, and alteration at the depths involved in the DeepStor project.

Reaction Kinetics
The third key parameter we did not include in this simulation was the reaction kinetics. The kinetics of mineral precipitation, dissolution, and alteration at geothermal conditions are condition-specific. They must be established for each mineral present in the system for the entire range of the system's temperature, pressure, and composition. A true, predictive reactive transport would include the reaction kinetics for all of the reactions identified in the process flow. A partial list of such reactions for this system is found in Section 4.2.3 below. Determining the kinetics of these reactions requires significant laboratory experimentation.
The purpose of this study was not to deliver a completely accurate, predictive model of scaling and formation damage risks during operation of the DeepStor project. Rather, our purpose was to identify minerals and processes involved in creating these risks. To this end, we identified several minerals that will play a role in scaling and formation damage, as well as potential pathways for these reactions to occur. These results set up phase 2 of our research, in which the kinetics of these reactions may be experimentally determined. The fact that the kinetics were not included in the simulation means that (1) the amount of precipitation predicted to return a solution at any given time step to equilibrium is over-estimated, and (2) it is possible, if not likely, that the kinetics of some of the predicted clay mineral alteration reactions will be too slow to occur in nature. Therefore, the results of this preliminary study should be taken as a worst-case scenario.

Potential Scale-Forming Minerals
A major objective of this study was to identify potential scale-forming minerals within the surface infrastructure of the proposed DeepStor project at the KIT campus. Scaling poses a risk to many geothermal plants, as it reduces the efficiency of heat exchangers by blocking passages and reducing heat transfer between fluids. Additionally, in extreme cases, scaling can become so pervasive that it completely blocks the movement of fluid to, from, or along the surface, at which point the system must be shut down for cleaning. This cleaning process is time-consuming, and taking a plant offline should be avoided in order to preserve the productivity of the system. Scaling risks are well known from geothermal energy production sites, but less so from HT-ATES reservoirs. The high temperature of the injection fluid at the DeepStor site is above almost all other HT-ATES systems' temperatures [7] and therefore needs to be investigated due to the retrograde solubility of calcite.
We identified calcite, anhydrite, antigorite, and amorphous silica as potential scaleforming minerals. Calcite and anhydrite are expected to form during fluid heating in the summer thermal recharge season. Amorphous silica may precipitate during the winter heating season, i.e., when the thermal fluid cools. Antigorite's precipitation kinetics are believed to be too slow to occur from a fluid. No effort was made to differentiate dolomite from calcite precipitation due to the low rate of precipitation of dolomite from natural systems (22). Anhydrite only occurs as a scaling phase during the period from 0 to 6 months at the hot well, after which no more anhydrite precipitates at the surface. The fluid is supersaturated with respect to amorphous at the cold well throughout this period, with 30 mg/kg required to be precipitated to bring the fluid back to equilibrium. If the fluid is able to equilibrate before reinjection, all of these minerals will be present in the surface infrastructure and will need to be removed.

Scaling Risk Mitigations
Anhydrite is a well-known geochemical risk in geothermal systems [27,28]. It poses a unique challenge, as it has retrograde solubility with respect to temperature, becoming less soluble with increasing temperature. This leads to precipitation during interaction within high-temperature portions of geothermal processes [28,29]. Anhydrite scaling also poses a problem in that acid treatments are often not effective at removing scale [28]. Instead, alkaline solutions and chelating agents have been used which effectively remove scale from pipes at geothermal conditions [29]. These chemical methods may be used in lieu of mechanical removal of scale.
Calcite is also a well-known scaling risk associated with geothermal systems [28,[30][31][32]. Calcite may be precipitated through decompression, reaction with pipe walls, boiling, and degassing of CO 2 [32]. In high-temperature geothermal systems, boiling is often a major concern [32,33]. Similar to anhydrite, calcite may precipitate in the DeepStor system during the summer thermal recharge season, due to its retrograde solubility. Several methods currently exist to deal with the scaling issues that calcite presents. One of these methods is controlling wellbore pressures with a deep pump. This method reduces calcite scaling and allows for control over the depth of scale deposition [33]. Additionally, chemical techniques have been developed which include the use of phosphate molybdenum inhibitors [32]. These treatments carry environmental concerns. Inhibitors based on polycarboxylates, however, have been shown to be effective at reducing calcite scaling risks under geothermal conditions [32,33]. Chemical inhibitors have been shown to be cost-effective because they mitigate the need for mechanical cleaning of the wellbore. Thus, they are a popular method of controlling calcite scaling [30,33].
No effort was made to differentiate dolomite from calcite precipitation due to the low rate of precipitation of dolomite from natural systems [22]. Calcite precipitation at the surface is variable during the summer thermal recharge, ranging from~1.4 to nearly 100 mg/kg. Anhydrite only occurs as a scaling phase during the period from 0 to 12 months at the hot well, after which no more anhydrite precipitates at the surface.
Amorphous silica is present as a scaling phase in this system, as it is in many geothermal systems. Silica species can be present in geothermal systems as amorphous silica, chalcedony, opal-A, and opal-CT [28,34,35]. Several chemical mitigation techniques have been developed for silica scales, many of which involve the modification of specific properties of the brine such as pH, or inhibition of colloid formation through rapid brine cooling [35][36][37]. There are also chemicals which inhibit the growth of silica colloids by adhering to the surface of a particle, thereby slowing further growth [34]. The dominant method of preventing silica scaling thus far has been to drop silica out of the solution at specified points within the system, e.g., designed settling ponds that often utilize seed crystals to accelerate the precipitation process [35][36][37].

Formation Damage
The potential for formation damage in this system is caused by mineral alteration in the reservoir. Anhydrite, daphnite, dolomite, minnesotaite, muscovite, quartz, saponite-Ca, saponite-Na, and siderite are all alteration products found in the reservoir during the simulation. Of the minerals forming in the reservoir, anhydrite, daphnite, saponite-Ca, saponite-Na, and siderite are not part of the original reservoir mineralogy. Their presence is the result of non-isothermal fluid circulation through the reservoir. Dolomite, minnesotaite, muscovite, and quartz occur naturally in the reservoir. Their concentrations also vary throughout these simulations. These minerals appear to have local phenomena affecting their precipitation and dissolution, leading to both elevated and depleted levels of each mineral at differing distances from wellbores in the model. Some minerals have a localized depletion point. Others have spikes of formation at discrete locations. The wellbores of this model act as local centers of mineral alteration.
Anhydrite forms only during the first 6 months of the simulation and is removed from the reservoir completely by 18 months. The initial formation of anhydrite is due to its retrograde solubility and the heating of the reservoir during the first injection periods. Its subsequent disappearance is due to its dissolution during a thermal extraction phase (cooling) and subsequent precipitation as mineral scale. The formation of anhydrite as mineral scale at the surface presumably removes enough sulfate from the circulating fluid to prevent additional anhydrite formation in the reservoir for the duration of the simulation.
All other alteration product minerals are present during every step of the simulation. Albite, annite, and calcite occur naturally in the reservoir, but their extent and location of their presence vary as the simulation progresses. This model shows albite and annite altering in the reservoir to clay minerals, namely, minnesotaite, saponite, daphnite, and muscovite. This alteration of albite and annite is detailed in the reactions found in Table 6 below. The products of these reaction products may serve as reactants in subsequent reactions. It is unlikely, however, that minerals are completely dissolved and precipitated. Rather, these reactions occur largely due to cation leaching and exchange reactions. Such reactions are also seen in the cation exchange between calcite and dolomite, which, combined with iron leached from annite, subsequently leads to the formation of siderite. Additionally, the porosity and permeability both increase at the hot well in each time step, likely due to the lower levels of alteration and removal of components from the surrounding medium at that location. Albite is altered to clay minerals during this simulation through reactions 1, 2, 3, and 4. Reactions 3 and 4 produce excess SiO 2 , resulting in potential silica deposition around the hot and cold wells whilst altering to daphnite and muscovite at the hot well. This is the only reactions shown above that produces SiO 2 ; therefore, we conclude that albite alteration is the source of the additional silica in the reservoir around the wellbores. This excess silica can also be consumed in reaction 2 to produce saponite. The alteration extent from the wellbores is the same as the alteration extent of daphnite, minnesotaite, muscovite, and saponite; therefore, albite is altered to all four of these minerals in these simulations.
Reactions involving annite are defined in this system by reactions 5, 6, 7, and 8 above. Annite may alter to the same minerals in this system as albite. This is evidenced by the spatial extent of alteration being the same for these four clay minerals and albite. Additionally, annite alteration is also evidenced by the increased annite concentration beyond its depleted zone around the wellbores. This annite increase corresponds to decreases in minnesotaite and muscovite at the same locations. Thus, annite is altered to clay minerals proximal to the wellbore, whilst clay minerals alter to annite directly adjacent to this zone. Partial dissolution of albite and annite provides the components that produce the clay minerals seen in these assemblages.
Calcite and dolomite display a leeching and exchange of magnesium, as well as the addition of bicarbonate to the solution. The presence of iron released by annite weathering, accompanied by the bicarbonate released by calcite and dolomite, leads to the presence of siderite. Siderite occurs around the cold well, reaching a peak after the first 12 months and subsequently reducing. Dolomite follows an inverse trend to siderite, increasing after each injection cycle at the area closest to each wellbore, while calcite is removed from these areas. Distally from either wellbore, calcite is the dominant carbonate phase. From these observations, we conclude that calcite is altered to dolomite at the wellbore, and dolomite is altered to calcite in areas further from the wells. No other carbonate minerals exist to produce the necessary ions for alteration mineral formation; thus, calcite is the sole mineral reactant that produces dolomite and siderite, as shown in reactions 9 and 10 above.

Conclusions
Here, we presented a preliminary scaling and formation damage risk assessment for the operation of a seasonal geothermal energy storage system in Germany's Upper Rhine Graben. We used equilibrium thermodynamics and reactive transport modeling to predict the composition and extent of mineral scaling in the surface infrastructure and formation damage in the reservoir. The identified scale-forming minerals included calcite and anhydrite during the thermal storage periods, when geothermal fluids are heated at the surface before injection into the storage aquifer. In the thermal consumption period, when geothermal fluids are cooled at the surface before injection, the identified scale-forming species was amorphous silica.
Formation damage in this simulation consisted of a series of alteration reactions, predominantly involving annite and albite. Silica and aluminum leach from each of these minerals, subsequently, along with the addition of iron from annite, leading to the formation of a constellation of clay minerals. These clay minerals may include daphnite, minnesotaite, calcium-and sodium-saponite, and muscovite. Cation exchange reactions between calcite, dolomite, and siderite were also found to occur in the reservoir. Overall, formation damage-causing reactions were found to increase the reservoir porosity and permeability immediately surrounding the hot wellbore, out to a distance of about 10 m. The permeability and porosity near the cold wellbore were both reduced.
The goal of this study was to identify minerals involved in scaling and formation damage in this system and to calculate an initial estimate of the total mass of these minerals' precipitation, dissolution, and alteration. In future studies, we will investigate the kinetics of these reactions to create true quantitative predictions of scaling and formation damage risks in this deep geothermal energy storage system. High-temperature flow experiments with synthetic reservoir fluids in a realistic temperature and pressure regime can verify the modeling results.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.