Porosity Dynamics through Carbonate-Reaction Kinetics in High-Temperature Aquifer Storage Applications

While near-surface geothermal energy applications for the heating and cooling of buildings have been in use for decades, their practical adoption is limited by the energy transport rates through soils. Aquifers provide a means to use convective heat transport to improve heat transfer between the building and the aquifer. However, the solid matrix in the aquifer is carbonaceous in nature, and calcification prevention techniques in the heat exchangers for the building also lead to dissolution of the aquifer matrix. Due to the Arrhenius nature of the reaction, dissolution rates may decrease with increasing temperature. An effective medium model is derived for the energy, calcium species, and fluid transport through a dynamic calcite porous medium which undergoes a reaction between the matrix and fluid. To better discern how these competing phenomena affect thermal transport in the aquifer, a two-dimensional Cartesian system is considered, where the vertical axis is parallel to the borehole axis, and flow is in the horizontal direction. An effective medium model is derived for the energy, calcium species, and fluid transport through a dynamic calcite porous medium which undergoes a reaction between the matrix and fluid. Since the fluid velocity decays algebraically with radial distance from the borehole axis, two flow regimes are considered. In one regime, far from the borehole where flow rates are small, conductive thermal transport acts faster than the species transport, leading to a case where precipitation dominates and regions of the smallest porosity contract to limit energy recovery. In regions with larger porosity, moderate advection of the species is sufficient to prevent significant pore closures over the time scale of exploration. The second regime, closer to the borehole, larger flow rates reduce species concentrations sufficiently to dissolve the solid phase between pores. In this second regime, Taylor dispersion effects in both energy and species transport compete, but thermal conduction acts more slowly than advection, promoting dissolution. The critical limitation in modeling the long-term evolution of the aquifer structure is the in situ dissolution rate.

transport acts faster than the species transport, leading to a case where precipitation

Introduction
Thermal regulation of buildings is a sizable portion of the energy budget worldwide, with residential and commercial sectors accounting for nearly 40% of the total U.S. energy consumption in 2018 according to the U.S. Energy Information Administration. (https://www.eia.gov/tools/faqs/faq.php?id=86&t=1). Finding an economically competitive technology that has a minimum carbon impact is crucial to reducing CO 2 concentrations in the atmosphere. One option is geothermal energy, and while largescale harnessing of energy near geologically active sites is a more familiar approach to capturing thermal energy from the earth, a lower-power alternative captures excess ambient heat in the summer and stores this in the earth, to be retrieved in the winter for heating. Small residential systems have been installed nationwide for nearly 40 years, and recently the start-up Dandelion Energy, started from Alphabet's X division, has begun to install these systems in the US (see https://x.company/projects/dandelion/). However, these systems are limited by the rate that energy can be delivered from the earth to the residence, so that if geologically the only means of energy transfer is conduction, a larger and more expensive piping network is required to retrieve or deposit the energy.
Aquifers provide a natural and promising solution to this thermal bottleneck. The thermal mass of the matrix allows for comparable energy storage as would be found for these conduction-dominated systems, with the benefit that the fluid connectivity throughout the aquifer provides a means to use the fluid as a coolant, greatly improving the heat transfer rates. One such system is in the Upper Malm formation within the Bavarian Molasse Basin in Germany. The formation is a carbonate-karstified groundwater aquifer which attains a vertical thickness up to 600 m (Prestel et al. 1991). While the Pullach site takes advantage of the higher-fluid temperatures at the deepest parts of the aquifer for geothermal energy production (Baumann et al. 2017), our work here centers on a carbonaceous aquifer in the northeastern part of the Bavarian Molasse Basin. Here, ambient temperatures are 22 • C, and have tight caprock formations above and below the aquifer (Dorsch and Pletl 2012;Fritzer et al. 2014). Karstification in this region allows better flow transport compared to the porosity of the calcite material found in the solid fraction. The goal for this site would be for the energy to be stored in the water fraction (Lee 2013).
At this site a research borehole was drilled, and a push-pull heat storage test was conducted as described in Ueckert and Baumann (2019). The heat storage test was performed at a flow rate of 15 L/s and with injection/production cycles ranging over 2.5 days (Ueckert and Baumann 2019). Nearly 48% of the injected energy was recovered. During the injection phase of the cycle, CO 2 was added to the injected fluid during heating of the water, before the heat exchangers. This injected solution was slightly supersaturated with respect to calcite at this elevated temperature, but was undersaturated at the ambient temperature of the aquifer. This led to significant dissolution of carbonates from the aquifer.
During the initial stage of this test, the time scales over which the temperature and species concentrations equilibrated were much slower than the advective transport near the borehole. The full reaction problem near the fluid/matrix interface is likely far from equilibrium in this case. Most experimental models of this reaction depend on the local species concentrations of the reactants and products, the coefficients of which depend nonlinearly on temperature which is spatially uniform Plummer et al. (1978), and under the assumption of weak gradients in species concentrations and temperature. Further, these dissolution rate models change in form radically over pH, while a consistent behavior of these models is that the dissolution rate decreases with increasing temperature and increasing calcium concentrations. In our model that follows, we assume that the net effects of the nonequilibrium effects over the timescales of temperature equilibration are negligible compared to the much longer time scale of the heat storage and retrieval test and application.
In particular, data related to the composition and flow properties of the Malm aquifer site centered on core samples from the drilling of the research borehole, and hydraulic testing once the borehole was established. Core samples consisted primarily of calcite and dolomite, and the pore scales of these samples were nearly impermeable to water. However, there was a 100 m section of the Malm aquifer for which no data were collected, and it is surmised that this section may consist of karst topography through which the matrix may be quite permeable. The hydraulic tests later confirmed this hypothesis, leading to a net transmissivity of 0.1 m 2 /s (Ueckert and Baumann 2019). The maximum local hydraulic conductivity of 6 × 10 −3 m/s, with a porosity of nearly 20% suggests a pore radius of nearly 1 mm.
Given that a high-temperature aquifer storage system is expected to last decades due to its high installation cost, a mathematical or computational model is needed to better understand the long-term impact of dissolution on the system performance. Any modeling description needs to represent the impact over meters of the dominant physics on the pore scale of 1 mm. The time scales of these processes range from minutes to days, and the goal would be to understand the net behavior over years. These disparate length and time scales show that direct computational approaches are simply not feasible to achieve all of these requirements. Further, given that the flow rate in the aquifer is inversely proportional to the radial distance from the borehole, this modeling approach needs to be applicable for a wide range of flow rates. The central scientific question of this work is How does calcite kinetics with fluid flow affect the thermal performance and structural integrity of the aquifer?
In this work, a mathematical description is derived on the meter-length scale (macroscale), which captures the average effects of hydrochemistry, hydrology and thermal transport over the fine-scale (microscale). To facilitate the modeling, a simple geometric structure to the solid phase of the aquifer is assumed, which allows for convenient analytical solutions of the unknowns in terms of variations over the macroscale. Closure of the macroscale dependence of the unknowns is determined through a modification of the asymptotic technique of homogenization (see Mei and Vernescu 2010 for an overview). By choosing an initially simple geometry for the solid phase, the permeability tensor components can be determined analytically based on the local pore size, which dynamically evolve over long time scales. While limited to the simple geometries, our approach allows for a quantitative understanding of the coupling of different physical processes and their impact on the energy storage performance of the aquifer. These reduced models are then solved computationally. Since the computational models are resolved on the macroscale, run times are significantly reduced, allowing for the exploration over a range of parameters. This research approach provides insight into the hierarchy of dominant mechanisms contributing to heat and mass transport within the system. Several reactive transport modeling approaches have been tried in the past decade. Classical homogenization techniques have centered on the reactant transport through a reactive media whose structure does not evolve during the process (for example, Allaire et al. 2010). These models require that the flow profile within the pore be laminar, or more precisely, that inertial terms in the momentum equations be much smaller than viscous effects in the flow. As the Péclet number of the system increases, however, secondary advection effects become pertinent in species or thermal transport. Two approaches in this limit are the extension of Taylor dispersion to porous media flows (see Salles et al. 1993), and hybrid homogenization approaches, such as the review from Scheibe et al. (2015), or more specifically for reacting flows, the hybrid model by Battiato et al. (2011). In this latter model, nonuniform reaction rates are determined over a mesoscale between the pore scale and the continuum scale. The relations between pore scale, CO 2 concentration, pore evolution in calcite and flow have been explored by Agrawal et al. (2020) for a single isothermal pore, experimentally and computationally. They find that not capturing the pore evolution on the 100 μm scale properly leads to underpredicting flow conductivities up to 80%, if their results are incorporated into a pore network model. While recent computational advances to multiscale reactive transport phenomena in porous media have been significant (e.g. see Molins and Knaber (2019) and Molins et al. (2020) for reviews of several approaches), these techniques are computationally intensive and not amenable to large parametric explorations. In all of these applications, energy transport is not included, and how effects on fluid flow during the reaction (e.g. dissolution and precipitation) mitigate energy transport is critical to the evaluation of the aquifer's potential for use in energy storage.
The remainder of this paper is organized as follows. In Sect. 2 the derivation of the effective medium model in two regions is presented: far from the borehole (species diffusion-limited) and near the borehole (species transport is advection limited). In Sect. 3 the results for both of these cases are presented, focusing on the dynamics during an initial injection phase, along with behavior of the system to periodic cycles of injection and productions. Conclusions are given in Sect. 4.

Problem Formulation
Our research methodology here is to formulate a canonical two-dimensional model of a meter-scale section of the aquifer that captures the dominant behavior of species transport and energy transport. Once this model problem has been formulated, the asymptotic techniques of homogenization and a regular asymptotic analysis are used to derive a reduced system of equations that describe species mass conservation, energy conservation, momentum conservation, and the effects of dissolution and precipitation on the porous structure. This reduced model is then solved computationally over the macroscale, and by assuming no significant variations in the solutions in the vertical direction, quantitative results can be achieved over a range of relevant porosities and flow rates.
Consider a series of channels within solid calcite, as shown in Fig. 1. Within each channel is a solution of water and ion species, whose composition depends on the local solution temperature and species concentrations that evolve depending on the dominant aqueous carbonaceous reactions. Further, along each solid-liquid interface, surface reactions occur which correspond to the carbon-based kinetics of dissolution of calcite. These reaction rates depend on the aqueous chemistry within each channel, the amount of CO 2 dissolved in solution, the local pH of the solution, and the local temperature. The goal of this modeling formulation is to incorporate the dominant aqueous and surface chemistry and energy conservation in order to understand the laminate dynamics within high-temperature aquifer storage applications. Reaction rates are necessarily temperature dependent, and knowing how each channel evolves over many storage cycles can aid in determining the energy storage capabilities of an aquifer.
A cell problem is defined as one fluid channel and one calcite laminate. The lower boundary of a characteristic channel is denoted y * = −H h 1 (x * , t * ), the upper boundary of this characteristic channel as y * = H h 2 (x * , t * ), and the upper boundary of the calcite laminate as y * = H h 3 (x * , t * ), where H is the characteristic thickness of the cell problem. Let h = h 2 + h 1 be the nondimensional channel thickness, or the liquid fraction of the cell, and h s = h 3 − h 2 be the solid fraction of the cell. For all of these functions, variations in x * are assumed to occur over a longer characteristic scale than the cell length (L H ); this assumption is intrinsic to the modeling approach for the effective problem derived below. In terms of the structure of the media, a volumetric constraint between liquid and solid fractions of the medium is assumed where λ is the modulation of the period between solid laminates and fluid-filled channels.
i , the molar concentration of species i, and they are governed by advection, diffusion, reaction equations where u is the fluid velocity, D i (T (1) ) = (ρ 1 k B T )/(6πμR), is the Stokes-Einstein relation for mass diffusion, with k B as Boltzmann's constant, and R is a characteristic radius of the ionic pair, T (1) is temperature of the solution, R * i is the reaction rate for the i th reaction per unit volume, which can depend on the other species concentrations and temperature. For simplicity, species diffusivity is assumed to be independent of temperature in the following.
From the reduced chemistry model Plummer et al. (1978), the following aqueous reactions take place in each channel using the the notation H 2 CO 3 * = CO 2(aq) + [H 2 CO 3 ].
Along each channel wall, dissolution of the solid matrix takes place based on the local temperature and local solid composition of the matrix. The aquifers of interest consist of a mixture of calcite and dolomite, and in this model, our focus is on the evolution of the calcite reactions. The forward reaction is modeled from Plummer et al. (1978), which is pH-dependent Note that the reaction times of (3)-(4) are on the order of seconds, compared to the reaction time of the surface reaction (6), which is on the diffusion time scale on the order of hours. Hence, the aqueous reactions are assumed to be quasi-steady, which results in the following equilibrium relations for the molar concentrations (Plummer et al. 1978) where K i is the equilibrium 'constant' of reaction i, where i = 1, 2, 3, 4. Each of these K i depends on solvent temperature T (1) in the form In solution within each channel, [Ca 2+ ] does not undergo reaction, so the concentration is governed by advection and diffusion The solution in the channels is electrically neutral, and the number of positive and negative ion species balance uniformly. Since the charge density for each species is the species concentration times the charge per ion, balancing positive and negative charge gives the relation 3 = [CO 2(g) ], and c (1) 4 = [H + ], and where the other species can be found from these quantities through (7).
The fluid mixture is convected in the positive x * -direction by the application of a time-dependent pressure drop of amplitude ΔP. The fluid temperature at the medium inlet is elevated compared to the ambient temperature T o by a temperature magnitude of T i −T o . Of interest is to capture the evolution of the channel walls due to the reaction kinetics between the solid matrix and the solution, while incorporating the temperature and flow dynamics over time on macroscopic length scales. These processes change the pore geometries, whose characteristic length scale across each channel, H , is much smaller than the macroscopic characteristic length scale along the channel, L.
For simplicity, the fluid mixture within the channel can be described mechanically by an incompressible viscous fluid (water), where u = ui + vj is the fluid velocity, ρ 1 is the fluid density, p is the fluid pressure, and μ is the fluid dynamic viscosity. Energy transport within the channel is dominated by convection and diffusion where c p1 is the heat capacity of the fluid, and k 1 is the thermal conductivity of the fluid. The fluid properties are assumed to be unaffected by the chemistry and the energy transport within the pore.
In the solid phase, energy transport is governed by conduction where ρ 2 , c p2 , and k 2 are the density, specific heat, and thermal conductivity of the calcite, respectively. Let c (2) 1 = [CaCO 3 ] as the molar volume fraction of calcite in the matrix. For mass conservation, per unit axial length in each channel where M (2) 1 is the molecular weight of calcite, and R * = k * K eq,s c (1) 1 /(K eq,s c (1) 1 + 1) corresponds to the surface reaction of dissolution of calcite, with k * being the reaction rate, Plummer et al. (1978). 1 follows that for Langmuir-isotherm models in surface-based reactions, with K eq,s as a surface equilibrium constant, in units of cm 3 /mole.
Finally, each channel wall is dynamic due to the reaction (15), which gives the kinematic boundary condition for the dynamics of the channel thickness requiring that the reaction rate at the lower boundary of the fluid channel be the same for each x * as the reaction rate at the upper boundary of the laminate. At each channel wall, the no-slip boundary condition, continuity of temperature, the balance of normal heat flux, and the mass balance of calcium for the reaction are applied (6), where q (1) j is the flux (diffusion and convection) of species j in the fluid. The vectors t and n denote the unit tangent and normal vectors along the pore wall Finally, the boundary conditions in the streamwise direction of cell problem, which depend on the flow conditions of the application, are given by Injection: Heated fluid from the well enters the channel and displaces fluid Production: Fluid is pulled from the system x * = L : k 1 T (1) In these conditions, P i is the applied pressure at the inlet during injection, P o < P i is the ambient pressure of the system (assumed to be applied at x * = L), and T i is the Table 1 List of characteristic dimensional values of parameters.

Quantity
Value (cgs) injection temperature of the fluid, while C i is the injection concentration of the calcium ion species. During production, P p < P o is the production pressure to remove the fluid from the system, and T o , C o are the values of the equilibrium temperature and calcium activity in the channel. Note that these conditions may not be representative of the actual conditions of the system during these processes, but do represent a reasonable first approach in order to understand the dynamics of the pore shape during injection and production. An asymptotic approach is used to reduce this problem to one that is computationally feasible on the time and length scales of interest. The derivation can be found in the Supplemental Materials document on the Journal's website. The resulting equations depend on the distance the domain of interest is from the borehole, since the characteristic velocity of the region is inversely proportional to the distance from the borehole based on mass conservation. Thus, far from the borehole, Pe c is small, and here is defined as Pe c = Pe c , where Pe c = O(1). The momentum balance in the fluid, species mass conservation, energy conservation, and channel wall dynamics are described by, to leading order h ∂c where Da is the Damköhler number,ᾱ is the effective thermal conductivity of the medium in the x 1 −direction,α is the effective thermal conductivity of the medium in the x 2 −direction, and v 12 is the calcium density ratio of the solution to the solid matrix. Our interest is in how the dynamics of the medium depends on different stages of energy storage and retrieval. These stages can be described in terms of energy injection: fluid at an elevated temperature and prescribed flux molar fraction of calcium ions entering the pore at x 1 = 0 at a prescribed velocity; energy production: fluid removed at a particular velocity with the temperature of the entering fluid x 1 1 determined by the local matrix temperature. Note that these boundary conditions need to be averaged in the x 2 -direction. Mathematically this gives us the following set of boundary conditions in x 1 for p o , θ o and c (1) 1 Injection: Heated fluid from the well enters the aquifer and displaces fluid Production: Fluid is pulled from the aquifer and enters the borehole For the dynamics closer to the borehole, the model (26) -(29) is not relevant, since the values of the Péclet number are larger than O( ). In this situation, interactions between the solutal or thermal gradients and the flow through the pores result in an additional streamwise diffusive effect in the mean which is commonly called Taylor dispersion (Taylor 1953). Extensions to porous media were initially done by Brenner (1980), but we follow Salles et al. (1993) in the following formulation, along with using the notation in Ortan et al. (2009).
When convection effects are promoted compared to the previous model, time needs to be rescaled in order to capture the effects of convection, For the boundary conditions, the new time scale yields the following conditions during injection and production cases: Injection: Heated fluid from the well enters the aquifer and displaces fluid Production: Fluid is pulled from the aquifer and enters the borehole

Results
Porosity changes in the aquifer are determined by the dissolution rate, which depends on the temperature and species concentration. The dissolution rate R depends on the concentration c (1) 1 and the temperature θ . The number of nondimensional parameters in the dissolution rate (see Supplemental Materials, (5)) is large, so these parameters are kept fixed throughout this study. In Fig. 2, a contour plot about the ambient temperature θ = 0 is shown. For R = 0, the equilibrium concentration of the calcium ions decreases with increasing temperature. Thus, if excess calcium is found in the entering solution, then precipitation within the porous medium is expected. Inlet conditions considered here are in the dissolution region of Fig. 2, which is consistent with scaling prevention within the borehole.
Our choice of ΔC is based on the applied temperature difference ΔT . In the application, the inlet species concentration of the fluid is held below equilibrium conditions in order to prevent scaling. ΔC is chosen so that if the temperature of the borehole fluid is elevated, the species concentration is below the equilibrium value at that elevated temperature. The relevant input concentration depends, critically, on the difference between the equilibrium concentration and ΔC. However, this equilibrium concentration depends on temperature, which changes based on the temperature difference ΔT of the incoming fluid and the ambient aquifer temperature. Because of this, the units for ΔC are either the fraction of the borehole equilibrium concentration at the borehole fluid temperature ΔT , or the fraction of the species concentration in the aquifer at ambient temperature.

Computational Approach and Assumptions
The model (26)-(29) is investigated subject to boundary conditions (30)-(31) or (32)-(33) using finite difference methods. To better understand the different physical mechanisms at work during injection and production, the simpler case when all of the unknown quantities h , p, θ, c (1) 1 are independent of x 2 is considered. Physically, this corresponds to a uniform distribution as a function of depth. The parameter values of the simulation are shown in Table 2, noting that the Péclet number ranges over three decades. These fixed parameters in the dissolution rate R, as defined in (5) in the Supplemental Materials document, result in a dissolution rate shown in Fig. 2. Note that near the equilibrium contour, dissolution increases with decreasing species concentration and temperature. This figure serves as a guide to see if a system at a particular location and time is in dissolution or precipitation.
Spatially, two grids are introduced on x 1 , each with spacing Δx 1 . One grid {Δx 1 , 2Δx 1 , . . . , N Δx 1 } where Δx 1 = 1/(N + 1) is called the Dirichlet grid, and our finite-difference representation of the pressure p is defined on this grid. Our second grid, called the Neumann grid, includes the endpoints of the domain x 1 ∈ [0, 1], with grid spacing Δx 1 , and on this grid our finite-difference representations of the remaining dependent variables h, θ, c The evolution of the system is discretized in time using a fully implicit Crank-Nicolson scheme on (27)-(29) with adaptive time-stepping. The algebraic system is linearized about the current time step, keeping all terms that are linear in terms of the delta formulation for a given time step Δt, which gives a truncation error of O(Δt 2 ). Then, the same problem is solved for two time steps each of duration(Δt)/2 each. This latter calculation is used for the starting value of the next time step, and this value along with coarser estimate of the unknowns at the larger time step is used to determine the next time step in order for the solution to be within a given tolerance.
For simplicity in this study, a common set of initial conditions is used. The pressure p and temperature θ are initially set to zero, h = 1, while the initial data for the c (1) 1 is picked to be the following form where c 1,eq is the equilibrium species concentration at the ambient temperature θ = 0 (far-field equilibrium within the aquifer), and ΔC is the inlet value of the species concentration. This choice of initial conditions for c 1 and of ΔC in (43) allows for a smooth, short, transient time over which these initial conditions relax to a relevant local quasi-equilibrium without inducing numerical instabilities.
Another aspect of this problem for which potential numerical instabilities need to be addressed is the transition between injection and production. To prevent instabilities, the following function of time is used as a "switch": This function has a transition region − log (10)/σ < t − t o < log (10)/σ over which f (t) undergoes a transition between 0.01 and 0.99. The choice of σ is flexible, but this value is chosen so that the transition between stages is short, but long enough in order for the simulation to be temporally well-resolved. This transition from injection to production, for example, is regulated by changing the inlet pressure from 1 to −1 over the transition period, and then the boundary conditions for θ and c 1 are then modulated based on the volumetric flow rate at that given time. One can extend this technique in the determination of ΔT , for example, when looking at an example similar to the heat storage test described in Ueckert and Baumann 2019 (see Sect. 3.3.2).
The code is designed to run over a fixed time interval, but there are limits to which the current model is physically pertinent. The system (26)-(29) is valid only when 0 < h < λ, or if the alternating vertical structure of fluid galley/calcite laminate is maintained. Hence, conditions are introduced at which the script terminates if h < h min = 0.001 or h s < h min . In the results shown below, this ending time of the simulation is a good metric for determining the dominant physical processes under certain operating conditions.
Lastly, Our results are presented in terms of the porosity φ and permeability K of the medium over time. To find these properties in terms of the variables of the current model, the extension to the capillary pore model are used (e.g. see Bear 1988;Costa 2016

P e c = O( )
Small values of the Péclet number correspond to low flow rates through the porous medium. These low flow rates, based on the anticipated flow conditions of the energy storage system, correspond to the response of the reacting matrix far from the borehole (hundreds of meters from the borehole). Although this location may not, at first, seem relevant to the short-term behavior of the system, over decades this section of the aquifer will change and these changes may affect the overall ability to transfer energy to and from the aquifer. Further, this limit provides a convenient framework for describing how the interacting processes of energy transfer, species transfer, and reaction interact, providing a conceptual framework for the description of the Pe c = O(1) results.
Recall that the reduced Péclet number Pe c is defined as Pe c = Pe c , Pe c = O(1), and our results are presented in terms of Pe c .

Injection
For the results in this subsection, an injection case is simulated, with the time of the simulation set to 0 < t < 8, which corresponds to a dimensional running time  Fig. 3 presents the porosity dynamics for this case when 0 < t < 6. In Fig. 3a, the spatial mean of the porosity shows a transient which increases slightly for small times but evolves to a steady decrease reflecting that precipitation dominates in the medium. Figure 3b shows the evolution of the spatial variation in porosity relative to its mean value. For short times, dissolution appears to be faster at the inlet of the flow compare to the outlet, but over longer times, the porosity at the inlet appears to finer compared to the outlet.
The mechanisms related to this change are due to the different rates of energy and solutal transport during the simulation. In Fig. 4a, the pressure relative to a prescribed linear pressure profile based on the boundary conditions at x 1 = 0 and x 1 = 1 is shown. Notice that for t = 6 the pressure deviation is sufficiently large so that there is no net flow beyond x 1 = 0.5. This suggests that a circulatory flow has developed in the region x 1 < 0.5, resulting in a reduction of solutal and thermal transport by advection. While the convective portion of thermal transport is reduced, conductive thermal transport persists. Further, the circulatory flow limits solutal transport through the medium, resulting in locally larger concentrations near the inlet. This local increase in concentration promotes precipitation near the inlet. As a confirmation of these different rates of thermal and solutal transport, in Fig. 4b, the evolution of the spatial average of the temperature scaled on ΔT (solid curve) and the spatial average of the species concentration relative to the far-field equilibrium concentration are shown. Note that the elevated temperature compared to equilibrium is significantly larger compared to the drop in species concentration, suggesting that precipitation would on average dominate over the medium. Hence, the lack of solutal transport compared to thermal transport results in precipitation. Secondly, consider the case where advection plays a more significant role. Consider the case Pe c = 40, φ o = 1/3 and operating conditions ΔT = 0.3 and ΔC = 99% of the equilibrium value. Figure 5a shows the evolution of the spatially averaged porosity over 0 < t < 6. In this case, porosity increases over this simulation, and note in Fig.  5b that the porosity of the medium is larger at the inlet compared to the outlet. Figure  6a shows the pressure deviation over this range, and note that the magnitude of this disturbance pressure is only 6% of the peak pressure, so flow through the medium is unimpeded. This higher permeability allows for faster solutal and thermal transport, which is seen in Fig. 6b. Note that the spatially averaged temperature equilibrates faster than the solutal average. This solutal average shows the net species balance between dissolution and advective transport. The combination of the elevated temperature and the reduced species concentration results in dissolution throughout the entire medium, which results in the eventual removal of the solid matrix at a finite time.
From these examples, the rates of solutal and thermal transport dictate the behavior of the porosity over time. In order to understand how these effects vary over different and φ o = 1/100, 1/50, 1/25, 1/20, 1/10, 1/8, 3/20, 1/5, 1/4, 1/3. This gives us 100 parameter value pairs, and all of these simulations are run with the same initial conditions described above for 0 < t < 8. After collecting these data, level curves of the mean porosityφ and the ratio of the minimum value of permeability compared to its initial value Figure 7 shows these parameter plots for the average porosity (a) and the scaled minimum permeability (b) for ΔT = 0.3 and ΔC at 99% of the inlet equilibrium value. Note that after t = 8, the value of the average porosity is slightly reduced for the porosity with Pe c = 100, but that minimum permeability decreases for φ o < 0.1. For values of Pe c < 70, the minimum scaled permeability is reduced for all values of φ o considered. This observation is consistent with the final porosity profile in Fig. 5a, since the minimum porosity at the outlet is about 0.25, which is significantly less than φ o = 1/3. In terms of flow and dissolution, this constriction leads to a lower pressure gradient for values of 0 < x 1 < 0.8, and then a steeper gradient for 0.8 < x 1 < 1. This lower gradient suggests a recirculation flow in the background of the applied pressure drop in the pores.
In Figs. 8 and 9, consider the simulations at the same Pe c and φ o , but now with ΔT = 0.2 and ΔT = 0.1, respectively. From Fig. 2, the dynamics are expected to promote dissolution compared to precipitation, and these plots reinforce this expectation. Note especially the middle figures that denote the scaled minimum permeability, where the regions for which K min /K o are greater than unity now drop below Pe c = 30 for φ o = 1/3. The figures on the right, denoting the average dissolution rate, also show that precipitation is confined for the smallest Péclet numbers and the smallest initial porosities. These curves suggest that poor solutal transport in the porous system relative to thermal transport is the driving mechanism for precipitation.

Long-Term Behavior
The benefit of our modeling choice is that long time scales pertinent to the field application are computationally feasible. Of particular interest is how the inlet temperature varies over time. During injection, a value near that of ΔT is expected, and during production, this value should decrease as fluid at an ambient equilibrium temperature enters the computational domain from the outlet. Poor advective transport could also result in less energy being produced, and this would be seen in our model as an inlet temperature which does reach the values of ΔT . Consider several cases of ΔT and ΔC for times over 0 < t < 80, which corresponds to a time range of 5 years. Injection and production periods are assumed to be the same 6-month duration. The data of interest centers on the behavior of the porosity, temperature, pressure and species concentration varies over time at the inlet. The same operating parameters for Pe c and φ o are considered as described in the previous section, and we use the same metrics to determine how permeability and average porosity varies over this time frame.
(a) (b) Fig. 9 Parameter plot of final average porosity (a) and minimum permeability relative to the initial permeability (b) for 10 < Pe c < 100 and 1/100 < φ o < 1/3, for ΔT = 0.1 and ΔC kept at 99 % of the equilibrium value at the inlet Figure 10 shows (a) the inlet pressure value; (b) the spatial average relative to ΔT or c 1,eq for the temperature and species, respectively; (c) the transient of the spatial mean of the porosity; and (d) the spatial profile of the porosity relative to its spatial mean over time for Pe c = 10, φ o = 0.01, ΔT = 0.1 and ΔC = 99 % of its equilibrium value. The pressure plot shows the seasonal variation for injection ( p(0, t) > 0) and production ( p(0, t) < 0). What is interesting to note here is that the porosity, temperature, and species concentration have nearly the same value at the inlet and the outlet, suggesting that the spatial variations of these quantities are minimal after transients evolve. Note that equilibrium values for the spatial mean are attained after approximately 1 year, while the porosity value continues to be reduced over nearly three years. Physically, the poor advective transport results in the energy being balanced by conduction, and the evolution of the species concentration and the porosity is governed by the reaction rate. Hence, poor permeability in the matrix results in a reduced retrieval of energy from the medium over time.
In contrast, Fig. 11 shows the case when Pe c = 40, φ o = 1/3, ΔT = 0.1, and ΔC = 0.99 of the equilibrium value at the inlet for (a) the evolution of the spatial mean of the porosity; (b) the spatial variation of the porosity with respect to the mean; (c) the evolution of the relative spatial mean of the temperature and species concentration over time. Notice that over the simulation the porosity is steadily increasing to above 80%, and the temperatures reach nearly the value of ΔT during the injection cycles. Further, the species concentration remains in the dissolution regime for values of θ shown. Hence for larger initial porosities and for sufficiently large reduced Péclet numbers, solutal transport in the void regions are sufficient to maintain a dissolution reaction until locally a solid layer is breached.
A summary of these results can be seen in Fig. 12, where (a) shows the mean porosity at t = 80 over a parameter range of 10 < Pe c < 100 and 1/100 < φ o < 1/3. Notice that for the Pe c = 10 that the mean porosity is reduced for values of φ o < 0.15, approximately. For Pe c > 10, dissolution appears to be the dominant behavior in this system, and the scaled minimum permeability in these cases results in much larger values of permeability. One anticipates from these figures that over much longer time scales, the dissolution would continue for the larger reduced Péclet numbers; this suggests that a porosity coarsening of the aquifer would occur as a result at these values of ΔT, ΔC. Regions with smaller initial porosity result in smaller values, since the solutal transport in these regions are limited. In regions where the porosity is larger and where convective transport of the species is sufficient, dissolution takes place unchecked by any contrary physical process.
In the case of ΔT = 0.1, as shown in Fig. 12, the cases are simulated for the full duration of the time interval 0 < t < 80. For larger values of ΔT , however, this is not the case, and the simulation for some Pe c and φ o terminate early. In Figs. 13 and 14 shows the mean porosity (a), and the scaled minimum permeability (b) when the simulation terminates. The time at which the simulation terminates is shown in the right-hand figure in both Figs. 13 and 14. Notice that for regions in the lower left-hand side of the domain (small Pe c and small φ o ), the termination of the program is due to precipitation within the pores, while in the upper right-hand side of the domain, dissolution is dominant. These regions, however, do depend significantly on ΔT . Comparing the regions where the solution ends before t = 80 for ΔT = 0.2 (Fig. 13) and for ΔT = 0.3 (Fig. 14), these regions of dissolution and precipitation both increase.  Fig. 12 Parameter plot of final average porosity (a) and minimum permeability relative to the initial permeability (b) for 10 < Pe c < 100 and 1/100 < φ o < 1/3, for ΔT = 0.1 and ΔC kept at 99 % of the equilibrium value at the inlet over 0 < t < 80

P e c = O(1)
While the same computational techniques described in Sect. 3.1 can be implemented to solve (35)-(38), the disparate balance between the advection and conduction terms in the energy equation (37)  taking small time steps in order to ensure stability of the solutions. The sampling of data in parameter space Pe c = O(1) require significantly more time and computational resources. Hence, the results presented here are more limited, but they reinforce the behavior seen in the previous section. Please note, for this section, the time variable t has been rescaled by 1/ , so a single time unit corresponds to around 33 minutes, compared to the time unit of 23 days in Sect. 3.2.

Injection
In this section, consider the integration of the system (35)-(38) subject to the boundary conditions (39)-(40) over the domain 0 < t < 5000. The initial conditions are the same as used in Sect. 3.2. The goal is to determine the evolution of the mean porosity and the minimum permeability, scaled on its initial value, after this evolution. Figure 15 shows an example of the evolution of species and energy transport for Pe c = 10, φ o = 1/100, ΔT = 0.3 and ΔC = 99 % of its inlet equilibrium value. In Fig. 15a, the evolution of the scaled temperature at t = 10 (solid blue curve) and t = 20 (dashed blue curve) are shown, superimposed with the scaled species concentration at (a) (b) Fig. 15 a Short-time evolution of scaled temperature profile at times t = 10 (solid blue curve) and t = 20 (dashed blue curve) and of the scaled species concentration at t = 10 (solid red curve) and at t = 20 (dashed red curve). Note that the transport of species occurs much more quickly than that of energy. a Longer-time response for the scaled spatial mean of temperature (solid blue curve) and species (dashed red curve) for 0 < t < 1000. Here Pe c = 10, ΔT = 0.3 and ΔC = 0.99 the same times. Note that profile of the scaled species curve translates nearly uniformly over this time interval, during which the temperature profile shows the slower migration of energy into the medium. These disparate transport rates are better displayed by looking at the spatial average values shown in Fig. 15b over 0 < t < 1000. The averaged species concentration decays nearly linearly from t = 0 to about t = 40, at which point a steady concentration profile takes shape. However, the scaled average temperature evolves on a much longer time scale, since for φ o = 1/100, energy transport is conduction limited. As noted in Sect. 3.2, larger initial porosities allow for more effective removal of the species due to better advection, which would lead to potentially better conditions for dissolution. However, this argument rests on rates of energy transport being larger than rates of species transport. With Pe c = O(1), however, larger initial porosities allow for better energy transport, so that larger Péclet numbers result in an initial transient for temperature and species to a nearly uniform profile, after which the porosity dynamics is governed primarily by the dissolution rate. With this argument, dissolution effects are expected to be more pronounced for smaller Péclet numbers and where porosity is small: regions in parameter space where energy transport is limited, as shown in Fig.  16. The trends for the final mean porosity (left figure) all point toward dissolution, but the changes with respect to Pe c are not significant. Changes in the scaled minimum permeability (Fig. 16b) are largest for the smallest values of initial porosity and the smaller values of Pe c .

Heat Storage Test
Lastly, consider a situation similar to that found in Ueckert and Baumann (2019) for the heat storage test. In this case, the system was run at a flow rate of 15 L/s over an exposed aquifer depth of nearly 300 m. Table 3 shows the nondimensional versions of the different stages of their heat storage test: the start and stop times, and the scaled temperature difference applied to the entering fluid. Note that the ambient temperature (a) (b) Fig. 16 Parameter plot of final average porosity (a) and minimum permeability relative to the initial permeability (b) for 10 < Pe c < 100 and 1/100 < φ o < 1/3, for ΔT = 0.3 and ΔC kept at 99 % of the equilibrium value at the inlet over 0 < t < 5000 of the fluid in the aquifer is 10 • C, and so ΔT chosen here may vary compared to the values at the site. In order to interpret the results, consider two characteristic values. The first represents ΔT which, as shown in Table 3, varies depending on the injection stage. In order to compare the injection and production temperatures for cases with different initial porosity, a characteristic value of ΔT = 0.3 is used for displaying the data only, so on a scaled temperature plot, stages 0 & 1 cannot achieve a relative temperature rise of 1, and for stages 2 & 3 it is possible to achieve a relative temperature rise above 1. This characteristic temperature scale then dictates the characteristic calcium concentration scale, since this is the equilibrium value of the species for temperature ΔT = 0.3. This characteristic concentration is denoted c 0,eq .
Only two different values of ΔC are considered in this section, the basis of which depends on the dissolution rate R. The first value, ΔC = 99 % of the equilibrium value, is chosen where ∂R/∂θ < 0, that is, higher temperature result in lower dissolution rates. This is the behavior expected based on the Arrhenius form of the reaction rate. The second value, ΔC = 50 % of the equilibrium value, corresponds to increasing dissolution rates for increasing temperatures, or ∂R/∂θ > 0. While a full parameter exploration would be ideal, the computational cost of such a study is still prohibitive, even with this reduced model.
1 as a function of x 1 at the start of each stage, as denoted in Table 3. ΔC = 99 % of the equilibrium value at the inlet, and φ o = 0.02. b Corresponding temperature profiles θ at the same times and conditions as (a) Our computational domain and approach here includes a domain on the field scale. The domain 0 < x 1 < 100 is chosen, which corresponds to a domain on the order of 10 m. Our resolution for the domain is coarser than our computations above, with a mesh size of Δx 1 = 0.2, and the time step determined adaptively to maintain second-order truncation errors. Further, we modify the initial conditions for the species concentration so that a transition region between the inlet concentration and the farfield concentration is well delineated This profile, and the resulting species profile at the beginning of each stage, are shown in Fig. 17a for the case when ΔC = 99 % of the equilibrium value, φ o = 0.02. Note that the species values change minimally outside of the transition region, due to the better advection transport. The corresponding temperature profile θ at each of these times is shown in Fig. 17b. Note that the thermal boundary layer near the inlet is gradually increasing at each stage, however, the domain over which the temperature varies significantly is much smaller than the domain over which the species concentration varies. This discrepancy impacts the reaction rate over this test, as we describe below. The case when ΔC is 99% of its equilibrium value for ΔT = 0.3 is considered in Fig. 18, where other conditions are prescribed in Table 3, for initial porosity values of φ o = 0.02, φ = 0.05, φ = 0.12, φ = 0.2, and φ = 0.33. Figure 18a shows the relative change of the average porosity over time for these cases, and overall these changes are minimal. An explanation of this result can be seen in Fig. 18b, where we show the final distribution of the porosity profile. These profiles are nearly identical except near the inlet, where the majority of the variation between the final profiles can be seen. Effectively, the species transport controls the dissolution rate away from the inlet, since the temperature has not deviated significantly from equilibrium. Near the inlet, however, the temperature variation is is significant. In Fig. 18c, we show the inlet temperature value over time for the duration of the heat storage test for each of the initial porosity cases. We note that the more porous cases result better energy transport, which leads to larger deviations of the pore structure near the inlet for  Table 3. Here, ΔC = 99 % of the equilibrium value of the solutal concentration c 0,eq at ΔT = 0.3. a Evolution of the relative inlet porosity growth; b Scaled final porosity profile; c Evolution of the relative inlet temperature; d Evolution of the inlet dissolution rate those cases. Finally, for this case we note that at equilibrium ∂R/∂θ < 0, and larger temperatures result in lower rates of dissolution, and this is reflected in Fig. 18d. While the dissolution rate is nearly uniform for the case φ o = 0.2, the cases corresponding to the largest values of φ o evolve to a lower relative inlet porosity compared to the smaller φ o values. The reason for this is that at later stages of the heat storage test, the dissolution rate becomes negative near the inlet region, corresponding to precipitation, and then the porosity is reduced near that region. These changes in the porosity provide a window into the energy production capabilities which are shown in Fig. 18c. Note that for the smallest porosities that the increase in productivity is slowly increasing over each injection/production cycle. Injection stages are where the inlet temperature is increasing, and production cycles correspond to reductions in the inlet temperature. Only when initial porosities are above φ = 0.2 do we see that the inlet temperature approaches the peak temperature for that stage. Hence, the capability of an aquifer to be used effectively for energy storage depends significantly on the local porosity.
In Fig. 19, the results of the heat storage test are presented in the case when ΔC is 50 % of the equilibrium value of c o,eq when ΔT = 0.3. In this situation, the dissolution rate increases with increasing temperature (∂R/∂θ ). There are two ways in which this difference in the dissolution rate behavior with temperature affects the pore dynamics. Firstly, there are no local conditions at which R < 0, so the matrix in the entire  Table 3. Here, ΔC = 50 % of the equilibrium value of the solutal concentration c 0,eq at ΔT = 0.3. a Evolution of the relative inlet porosity growth; b, Scaled final porosity profile; c Evolution of the relative inlet temperature; d Evolution of the inlet dissolution rate domain is dissolving over the heat storage test. This global variation with respect to the initial porosity is shown in Fig. 19a, since the thermal penetration depth at the inlet increases with increasing porosity, dissolution is then enhanced due to the locally higher temperatures, and hence larger porosity variations on average are observed. Secondly, as the thermal penetration depth increases from the inlet, higher dissolution rates in this local region are expected, and this effect is shown in Fig. 19b. Finally, the inlet evolution of the temperature, shown in Fig. 19c and the inlet dissolution rate, shown in Fig. 19d, are in phase, suggesting that the dissolution rate at the inlet is governed primarily with the changes in the local temperature.
dispersion effects are pertinent, a robust description of the competition between the dominant mechanisms is found. A direct numerical simulation of the full problem is simply not feasible. From an order analysis, a direct simulation of the full problem would require 1/ more mesh points per spatial dimension compared to the effective model, even with neglecting Reynolds stress terms in the Navier-Stokes equation. The current computational scheme of the effective model does not lead to a banded matrix for the Crank-Nicolson solver, suggesting that the time step would have to be reduced by 2 in order to maintain the same accuracy as the effective model. Given that = 10 −3 , 10 6 more steps per time unit would be required, and each time step would take between 10 3 and 10 6 times the computational operations of the effective model to solve. Since this problem is a transient problem, parallelization, either through multiple processors or the use of GPUs, would be of limited use to speed up the simulation.
The ability to use aquifers to store and recover thermal energy depends on the capability of harnessing the energy and the integrity of the overall carbonaceous structure. Both items depend on the interaction of energy, species and momentum balances taking place during operation of the geothermal energy site. The species balance depends on the rate of dissolution, and comparatively the rates of diffusion and advection of the species. The energy balance relies primarily from this modeling on convection and conduction through the medium, while the momentum balance in this study pits viscous stresses with normal velocities induced by the pore dynamics instigated by the dissolution kinetics. Hence, our focus in this work has been on the relative rates of species, energy and momentum transport compared to the rates of the pore dynamics and dissolution.
From this model, the relative rates of energy and species transport dictate the dynamics of the pore structure to first order in the injection stage. When energy transport is dominated by conduction and the species transport is limited, precipitation can take place, preventing further reductions in the local species concentration. Physically, this corresponds to small values of the Péclet number and small porosity. For larger values of porosity, number, species concentrations are balanced by dissolution, advection and diffusion, where the dissolution source can be limited as the temperature within the porous medium increases. The reduced calcium species concentration of the entering fluid plays a significant role in these dynamics according to this model. So, for small Pe c , the smallest pores become smaller while the larger pores become larger, and potentially, change the representative length scale over which the porous medium derivation takes place. For significantly reduced inlet concentrations compared to its equilibrium value, the model's predictions are limited due either to porosity reductions by three orders of magnitude, or to nearly complete removal of the solid phase.
For larger flow rates (Pe c = O(1)), Taylor dispersion plays a role in increasing diffusive transport in both energy and species. This increased transport affects the species transport more quickly than energy transport, resulting in dissolution over shorter times until the temperature in the porous medium equilibrates. The effects on the pore dynamics on the fluid flow are reduced in this case as well. From our model, higher rates of advective species and conductive energy transport of energy and species result in a reduced impact of the evolution of the pore structure.
For this system to be pragmatic for a duration on the order of decades, the long-time evolution of the matrix structure is critical to the determination of energy recovery. In the case of smaller flow rates, viability depends on the degree to which the input species concentration is below its equilibrium value. For inlet species concentrations near its inlet equilibrium value, significantly large permeability changes are predicted over a 5-year cycle. Note that during the production phase of this simulation, the entering fluid is assumed to be at the equilibrium concentration of the initial natural ambient value. It is unclear how the species concentration of the entering fluid is closer to the equilibrium value at the elevated temperature would affect the pore dynamics over these time scales. For larger flow rates, the computational requirements of the current model are significantly higher, and the limit of a few months to one year is the practical limit of our study. Note that our simulation of conditions analogous to the heat storage test of Ueckert & Baumannn Ueckert and Baumann (2019) shows that the effects of the pore structure may be manageable, based on the dependence of the dissolution rate on species and temperature. Future work in this area centers on simulation of the effective model in two dimensions, to see the effect of thermal transport and borehole flow rate variation on the pore structure of the aquifer. Extensions of this model to cylindrical coordinates to capture the effects of the flow rate decaying radially from the borehole would be of interest.
In total, this study suggests that the most critical limitation is the understanding of the dissolution rate under field and operating conditions. The choice here was quite sensitive to the parameter values, and as shown in the heat storage test simulations, slight variations in the modeling choices can result in effects which may lead to significant changes in the porosity evolution. Efforts to a better dissolution rate estimate in this regard would help greatly in connecting how flow conditions can lead to behaviors in high temperature aquifer storage.

Funding
The authors gratefully acknowledge the support from the Bavarian Ministry of Economic Affairs, Energy and Technology and the BMW Group.

Declarations
Conflict of interest The authors declare that they have no Conflicts of interest.
Code availability Computational results performed using custom code.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/ by/4.0/.