Impact of Mineral Reactive Surface Area on Forecasting Geological Carbon Sequestration in a CO 2 -EOR Field

: Mineral reactive surface area (RSA) is one of the key factors that control mineral reactions, as it describes how much mineral is accessible and can participate in reactions. This work aims to evaluate the impact of mineral RSA on numerical simulations for CO 2 storage at depleted oil ﬁelds. The Farnsworth Unit (FWU) in northern Texas was chosen as a case study. A simpliﬁed model was used to screen representative cases from 87 RSA combinations to reduce the computational cost. Three selected cases with low, mid, and high RSA values were used for the FWU model. Results suggest that the impact of RSA values on CO 2 mineral trapping is more complex than it is on individual reactions. While the low RSA case predicted negligible porosity change and an insigniﬁcant amount of CO 2 mineral trapping for the FWU model, the mid and high RSA cases forecasted up to 1.19% and 5.04% of porosity reduction due to mineral reactions, and 2.46% and 9.44% of total CO 2 trapped in minerals by the end of the 600-year simulation, respectively. The presence of hydrocarbons affects geochemical reactions and can lead to net CO 2 mineral trapping, whereas mineral dissolution is forecasted when hydrocarbons are removed from the system.


Introduction
Geological carbon sequestration (GCS) is a critical component in accomplishing the goal of net-zero or carbon neutrality set by governments and industries [1][2][3], as it provides an enormous estimated storage capacity, and its efficacy has been successfully demonstrated many times by pilot-scale and field-scale projects worldwide [4]. Among the storage options, GCS at depleted oil fields, especially enhanced oil recovery with CO 2 (CO 2 -EOR), has drawn a lot of attention, because of the higher economic incentives (in addition to government subsidies, e.g., the 45Q tax credit in the U.S.) and the lower costs of characterization, construction, and deployment. While most CO 2 trapping mechanism analyses were performed for GCS at deep saline aquifers (e.g., [5]), the essence of the trapping mechanisms is similar for CO 2 -EOR projects [6][7][8]. However, one particular trapping mechanism, the mineral trapping of CO 2 , is often ignored in storage forecasts for CO 2 -EOR projects, such as Jia et al., 2016 [6], even though mineral trapping is the most secure mechanism to sequester CO 2 in the long term. Due to the lack of geochemical modeling in these forecasts, the processes of CO 2 dissolution and its presence in aqueous ions are also ignored. However, understanding the geochemical interactions between CO 2 , in situ fluid, and formation, is critical to ensure long term CO 2 conformance and to mitigate the risks of groundwater contamination due to CO 2 leakage. Mineral reactions are numerically described via two key parameters: chemical equilibrium constant and mineral dissolution and precipitation reaction rate. For a given reaction, the chemical equilibrium constant is a function of temperature, whereas the mineral reaction rate is controlled by further factors. As defined in the Transition State Theory (TST), the mineral reaction rate is a function of reactive surface area (RSA), activation energy, temperature, and equilibrium constant [9]. Particularly, the mineral RSA is crucial as it determines how much mineral is accessible and can participate in the reactions, yet it has been a challenge for it to be fully characterized. There are several approaches to measuring the mineral RSA (or specific surface area) in the laboratory, such as the widely used Brunauer-Emmett-Teller (BET) method [10], and newer image-based methods [11,12]. However, neither of these methods is ideal for characterizing mineral RSA. For example, the geometry-based method is constrained to the assumption of uniform-sized mineral grains, and thus its performance is compromised for the heterogeneously shaped mineral grains; on the other hand, the image-based methods need a drastically large number of images to capture the evolving mineral RSA over time. Therefore, it is common to observe order of magnitude differences in RSA measurements for the same mineral [13]. Such wide ranges of uncertainty pose great challenges in selecting parameters for numerical simulations. Even worse, most GCS simulation tools that include geochemical modules (e.g., CMG-GEM, TOUGHReact) ignore the spatial and geometry heterogeneity of mineral RSA values and simplify the temporal variation of mineral RSA values.
Existing uncertainty studies for GCS primarily focus on uncertainties stemming from porosity, permeability, and operational factors (e.g., injection scheme) [6,[14][15][16][17][18][19][20]. Only a few studies have examined the impact of uncertain reactive surface areas on GCS performances. Qin and Beckingham (2021) compared simulation results of a core-scale model using mineral RSA obtained from different methods [21]. Luo et al., 2012 investigated the effect of the RSA of calcite and anorthite using a two-dimensional (2D) generic deep saline aquifer model and confirmed that these parameters have a significant impact on CO 2 mineral trapping [22]. Bolourinejad et al., 2014 followed a similar approach to evaluate the impact of the RSA of seven minerals on CO 2 trapping for a depleted gas field in the Netherlands [23]. It was found that the RSA of quartz has the greatest impact on CO 2 mineral trapping; however, neither hydrocarbon components nor realistic field operations were taken into consideration in their numerical model. Moreover, both Luo et al., 2012 and Bolourinejad et al., 2014 employed only one well (the injector) in their simulations. Recently, Jia et al., 2019 studied CO 2 trapping mechanisms, including mineral trapping, for an ongoing CO 2 -EOR project with multiple five-spot well patterns, but only briefly discussed the potential impacts of using different mineral RSA values in their simulations [22][23][24].
This work aims to evaluate the impact of mineral reactive surface area on GCS numerical simulations where hydrocarbons are present. The novelty of this work includes: (1) using a three-dimensional (3D) field-scale reservoir model to eliminate the boundary effects imposed by 2D models; (2) using site-specific formation and fluid properties as well as realistic operational activities, i.e., CO 2 -EOR with multiple wells; and (3) taking hydrocarbon components into consideration, thus simulating reactive transport for the three-phase system (water, gas, and oil). The Farnsworth Unit (FWU) site in northern Texas is selected as a case study. As the study site of the Southwest Regional Partnership on Carbon Sequestration (SWP) Phase III, the FWU site has been investigated with many characterization and monitoring techniques (e.g., [25][26][27][28]) and numerical analyses (e.g., [7,8,14,15,24,29,30]).
The rest of this paper is structured as follows: Section 2 describes the FWU reservoir model and the methods to evaluate the impact of mineral reactive surface areas; Sections 3 and 4 present simulation results and discuss how the presence of hydrocarbon components and uncertain RSA values affect reactive transport and CO 2 mineral trapping; and Section 5 summarizes the findings of this work.

FWU Description
The FWU site is located in the Anadarko basin in northern Texas (Figure 1a). Following the primary and secondary production that started in the 1950s and 1960s, the FWU has been currently undergoing CO 2 -EOR since December 2010. As of 2020, more than one million metric tons of CO 2 has been sequestered at the FWU [31]. The CO 2 storage system at the FWU consists of late Pennsylvanian Morrow B sandstone as the storage formation and overlying Atokan Thirteen Finger limestone as the sealing layer. With the support of the U.S. Department of Energy and the National Energy Technology Laboratory, the SWP team has been conducting a variety of characterization; monitoring, verification, and accounting (MVA); simulation; and risk assessment research activities over the past decade. For the sake of brevity, please refer to previous publications for the details of CO 2 sequestration at the FWU. , 2020 addressed uncertainty analysis due to both geological and operational factors and quantitative risk assessment of CO 2 leakage and its impact on overlying underground sources of drinking water [16,17,[37][38][39]. This study is built upon this previous FWU work and focuses on risk assessments associated with reactive transport, in particular mineral reactive surface areas.

Mineral Reactive Surface Area
The reactive surface area of minerals is a dynamic property that varies from mineral to mineral, changes over time as geochemical reactions develop, and depends on the heterogeneous distributions, various shapes, and complex contact interfaces of all minerals in the subsurface. A variety of methods are available to characterize mineral reactive surface areas, such as the well-known BET (Brunauer-Emmett-Teller) approach, analytical methods based on kinetic experiments, and the image-based approaches that rely on scanning electron microscopy (SEM) and computed tomography (CT) techniques [11,40]. These methods have been used for measuring mineral reactive surface areas for GCS purposes, e.g., [40,41]. As expected, the measured reactive surface area of the same mineral drastically varies from method to method and from site to site. For example, Bolourinejad et al., 2014 measured the specific surface area of kaolinite in the Rotliegend reservoir cores in the range of 1.2 × 10 6 m 2 /m 3~3 .4 × 10 7 m 2 /m 3 ; Beckingham et al., 2016 reported different specific surface areas of calcite, ranging from 8.1 × 10 4 m 2 /m 3 , measured by the BET approach, to 7.6 × 10 5 m 2 /m 3 , measured by image-based methods [11,23]. The term specific surface area (SSA) is generally used to describe the reactivity of the pure mineral, while the term reactive surface area (RSA) is usually referred to as the average reactivity of each mineral in the porous media. Therefore, the SSA values need to be converted to RSA values with site-specific mineral volume fractions. At the FWU, there are seven key minerals identified from core analyses. A survey of the reactive surface area of these minerals was performed and is summarized in Table 1. Please note that mineral volume fractions (listed in Table 2) were used to convert the SSA values reported in other studies to FWU site-specific RSA values. The several orders of magnitude difference between the lower and upper values in Table 1 reiterate the primary research question of this study of whether the choice of a reactive surface area value impacts the results of reactive transport simulation and to what extent. Table 1. Ranges of reactive surface area (RSA) of seven minerals from a literature survey [11,[22][23][24]40,41]. While choosing a value from the wide range for a certain mineral is already difficult, describing the changing characteristics of the mineral reactive surface area is no easier. The common practice in reactive transport simulations is to calculate the reactive surface area with the following equation:

RSA (m 2 /m 3 ) Calcite Kaolinite Dolomite Quartz Ankerite Siderite Illite
where A i is the reactive surface area of mineral i at current time step, A 0 i is the reactive surface area of mineral i at time 0 (i.e., initial value), N i is the mole amount of mineral i per unit grid block volume at the current time step, and N 0 i is the mole amount of mineral i per unit grid block volume at time 0. There are two main restrictions on using this approach. Firstly, it assumes a uniform distribution of reactive surface area for the same mineral across the entire model domain (which could be kilometers for GCS sites); in other words, a homogenous reactive surface area is assigned to each mineral. Secondly, the controlling factor of reactive surface area at the current time step, N i , is determined by the current mineral dissolution and precipitation rate that is affected by the mineral reactive surface area at previous time steps, as shown in Equation (2): where r i is the reaction rate for mineral i, k i is the rate constant of mineral reaction i, Q i is the activity product of mineral reaction i, and K eq,i is the chemical equilibrium constant for mineral reaction i. Therefore, deviations between the estimated and real reactive surface areas accumulate over time and may lead to significantly different mineral reaction predictions after hundreds or thousands of years. Given that an exhaustive dataset for the mineral reactive surface area is not available, and neither is an advanced reactive transport numerical model framework, an uncertainty analysis is probably the best approach to investigate the impact of mineral reactive surface area on GCS predictions. Based on the ranges listed in Table 1, a seven-factor Box-Behnken Design was developed to cover the entire uncertain space with 85 cases. The low (−1), mid (0), and high (+1) values for each of the seven minerals are presented in Figure 2. Please note that the mid values are not the average of the low and high values, but rather the median taken from the RSA values of the literature survey. Two more cases, (−1, −1, −1) and (+1, +1, +1), i.e., all low values and all high values, were added to the simulation design. A full list of the cases is presented in Table A1. It is anticipated that the design of these permutations will discover the impact of mineral reactive surface areas, with the mutual dependencies among the mineral reactions taken into account.

Numerical Models
The 3D FWU reservoir model consists of 82 cells in the x-direction, 78 cells in the ydirection, and four cells in the z-direction, with a total of 25,584 grid blocks. Each grid block is 200 ft (or 60.96 m) by 200 ft (or 60.96 m) in x and y directions and about 8.67 ft (or 2.64 m) in the z-direction. This model only includes the Morrow B sandstone, as the integrity of the overlying Thirteen Finger limestone will not be compromised for at least 5000 years [31]. A total of 55 wells (23 injectors and 32 producers) were distributed across the model domain in 5-spot patterns ( Figure 1b). This reservoir model has been calibrated with the FWU field history by the SWP [34,42]. Well schedules include a 10-year CO 2 -EOR period, a subsequent 10-year post-EOR CO 2 injection period, and another 580-year monitoring period. In the CO 2 -EOR period, CO 2 and water were alternatively injected (water alternating gas (WAG)), reflecting the site history. In the post-EOR CO 2 injection period, all producers were shut-in, and CO 2 was continuously injected at 5 × 10 6 ft 3 /day (or 1.42 × 10 5 m 3 /day) through all injectors. In the monitoring period, all wells were shut-in, and no CO 2 was injected into the model. The heterogeneous porosity and permeability distributions of the reservoir model were derived from the site characterization and geostatistical analysis of previous SWP work. Figure 1b presents the heterogeneous porosity distribution of the 3D FWU model and the locations of the wells.
Seven major minerals were identified from FWU cores, as listed in Table 2. The measured volume fractions for each mineral are used as the initial conditions for the reservoir simulations. The SWP team periodically monitored the water quality of samples from the shallow aquifer and produced water collected at the FWU site. In this work, the averaged measurements of the produced water are used as the initial conditions of the aqueous species, as shown in Table 3. The compositions of hydrocarbon components were analyzed from FWU oil samples and are also listed in Table 3.
Seven mineral reactions and three aqueous reactions were modeled in this work, as shown in Table 4. The chemical equilibrium constants (log K eq ) were calculated for reservoir temperature at 70°C with the EQ3/6 database [43], and the activation energy and the standard TST rate constant (log k 25°C ) values were taken from Palandri and Kharaka (2004) [44]. In particular, fourth-order polynomials were fitted to the EQ3/6 database of chemical equilibrium constants of the reactions versus temperature (please see Figure A1). Table 4. Mineral and aqueous reactions and parameters for reactive transport modeling [43,44].

Reaction
Activation Energy (kJ/mol) logK eq at 70°C logk 25°C (mol/m 2 s) Reactive transport models require more computational resources than multiphase flow models because of the additional governing equations for aqueous and mineral reactions. Therefore, simplified modes are often used for reactive transport simulations. A similar approach was followed in this work. Instead of using the abovementioned reservoir model for all 87 simulations, a screening study was conducted first, which used a simplified 3D model to evaluate the performances of all 87 combinations of mineral reactive surface areas. Then the selected combinations, which are either representative or extreme scenarios, were repeated with the full FWU reservoir model. The simplified 3D model effectively represents a batch simulation. In particular, there are three cells in each direction (x, y, and z), with infinite acting aquifers around all boundaries to avoid the boundary effect. A CO 2 injection well was placed in the center of the model to introduce supercritical CO 2 to the model at a bottom-hole pressure constrained injection rate. Except for the grid structure and well configurations, all other model settings are identical to the full-size reservoir model. Simulation results at the center of this simplified model can be used to screen the combinations of mineral reactive surface areas. In addition to the much lower computational cost, the other advantage of this approach is that it can isolate control factors to simulation results, e.g., the domain size and the interference of multiple wells, and therefore focus on discovering the impact of mineral reactive surface area.
All simulations were performed with the 2019.10 version of the CMG-GEM simulation package, the flow equations of which are listed in Appendix B.

Screening the RSA Combinations
The 87 combinations were simulated to study the ranges of uncertainty in mineral precipitation/dissolution and CO 2 mineral trapping caused by using different mineral RSA values. Figures 3 and 4 present the simulation results of the simplified model and emphasize three RSA combinations, Case #85, Case #86, and Case #87, where all mineral RSA values were set at mid, low, and high, respectively (please refer to Table A1 for the detailed RSA combinations). The predicted pH values of all cases exhibit a similar pattern (Figure 3a). A drastic drop from the initial pH (~6.7) to the range of 4.0~4.4 during the first day is followed by a relatively slower recovery through the rest of the simulation period. While most of the curves overlap each other, the dashed dark green line, which indicates the results of Case #86, presents a deeper pH decrease in the beginning and a slower pH recovery rate. Comparing the much faster pH recovery rate of the other cases, it is obvious that the prolonged recovery time of Case #86 results from the lack of buffers, i.e., fewer accessible minerals to react with the H + introduced by CO 2 .  Because the changes of the silicate minerals are much smaller than the total amount of silicate minerals in the model, the mineral reaction rate is a better parameter to illustrate the silicate mineral reactions. The predicted precipitation or dissolution of silicates at the center of the model are shown in Figure 3b-d. For illite and quartz, Case #87 (dashed dark blue line) and Case #86 (dashed dark green line) is the most and the least reactive scenario for these two minerals, whereas Case #85 (dashed dark red line) overlaps with most of the medium cases. On the other hand, changes in kaolinite are more complicated. There are many crossovers and overlaps among results using different kaolinite RSA values, and the deviations between lines in the same color are much greater than those in Figure 3b or Figure 3d. This suggests that the reaction of kaolinite is controlled by not only the RSA of kaolinite but also the RSA of other minerals. Nevertheless, the high RSA cases (blue) predicted prompt and more intense reactions, while the low RSA cases (green) predicted delayed and more restricted reactions.
The simulation results of four carbonate minerals, ankerite, calcite, siderite, and dolomite, are shown in Figure 4a-d, respectively. Please note that all 87 cases started with the same initial mineral amounts, and Figure 4 presents mineral amounts since Day 1; therefore, the differences between the starting points of the lines in each subplot indicate the difference in the reactions that occurred during the first day of simulation. For example, the total amount of ankerite in the model was initially about 295 kmol, and dropped to 74 kmol in Case #87 (dashed dark blue line), whereas it dropped to about 245 kmol in Case #85 (dashed dark red line) and to 292 kmol in Case #86 (dashed dark green line), after the first day of simulation.
For all 87 cases, a monotonic dissolution of ankerite is observed, while all three other minerals show monotonic precipitation within the 600 years of simulation. The low RSA cases (green lines) once again exhibit the least reactivity, i.e., the slowest to deviate from the initial values. Interestingly, the prediction results fall into three clusters for each of these four minerals in the first year. In other words, the RSA value of the particular mineral, being either low, mid, or high, determines the reaction of this mineral in the first year, regardless of what RSA values are used for the other six minerals. Moreover, there is zero or small difference among predicted mineral amounts using the three levels of RSA values for ankerite ( Figure 4a) and siderite (Figure 4c) in the later stage of the simulation period. The differences among the three clusters of lines are more obvious for calcite ( Figure 4b) and dolomite (Figure 4d) after 300 days. At the end of the simulations, the cases with the high dolomite RSA value forecast much greater dolomite precipitation than the other cases.
Compared to the silicate minerals (Figure 3b-d), a stronger clustering effect is observed in the carbonate minerals, especially in the first two years of simulation, suggesting that the reaction of the carbonate mineral (Figure 4a-d) is more controlled by its own RSA value. The interference (or mutual dependence) between mineral reactions is weaker for the carbonate minerals than for the silicate minerals. As a result, Cases #85, #86, and #87 seem to be able to represent the median, minimum, and maximum reaction scenarios of all tested cases for all minerals.
The mineral trapping of CO 2 in the simplified model was evaluated with all 87 cases, and the results are presented in Figure 5. There is a much greater deviation among total mineral trapping amount predictions than the predictions for each mineral, as shown in Figures 3 and 4. No clear clustering is observed in Figure 5. While the impact of RSA on each individual mineral is apparent and easier to interpret, its impact on the CO 2 trapping mechanism is rather complicated. This is likely due to the different abilities of each mineral to sequester CO 2 . For example, one mole of precipitated calcite effectively secures one mole of CO 2 , while for dolomite the ratio becomes 1:2. At the end of the simulation, the amount of CO 2 trapped in the minerals ranges from 2 kmol to 200 kmol. The results of three representative cases (#86, #86, and #87) are very close to the median, minimum, and maximum values of all cases, as shown in the dashed lines in Figure 5. Therefore, these cases were selected for reactive transport simulations with the FWU reservoir model.

Reactive Transport Prediction with the Reservoir Model
As minerals precipitate or dissolve, it can change the porosity of the storage formation and thus affect fluid flow patterns and subsequent mineral reactions at new fluid-rock contacts. Therefore, we analyzed the porosity change due to mineral reactions in the FWU model. Figures 6-8 present the changes in porosity at three critical time steps: the end of the CO 2 -EOR period, the end of the post-EOR CO 2 injection period, and the end of the simulation period, respectively.   After ten years of CO 2 -EOR operation, porosity reduction is observed in all three cases. Figure 6 shows changes of porosity due to mineral reactions at the top layer of the Morrow B sandstone using the RSA values of Case #86 (all low RSA values of seven minerals), Case #85 (all mid RSA values), and Case #87 (all high RSA values), respectively. In general, a dominant porosity reduction was observed for all cases. Specifically, Case #86 (Figure 6a) exhibits the most restricted variation, while Case #87 (Figure 6c) presents the greatest porosity change across the layer, showing a maximum porosity reduction of 1.86 × 10 −3 (or 0.9% of the initial porosity). However, very similar CO 2 global mole fraction distributions were predicted by three models (Figure 7). This suggests that the mineral reactions have an insignificant impact on the CO 2 migration. Given that the CO 2 plume shapes are mostly identical among the three cases, the difference in porosity change is attributed to the mineral RSA values.
The CO 2 -EOR period was followed by a post-EOR CO 2 injection period, during which all production wells were shut-in, and CO 2 was continuously injected via all injection wells for ten years. Figure 7 presents the simulation results of porosity changes due to mineral reactions. There are visible changes in porosity reduction comparing Figures 8b and 6b, and  Figures 8c and 6c. A clear expansion of areas with porosity loss is observed. The greater porosity reductions occur along the edges of the CO 2 plumes (Figure 8c). Specifically, after ten years of CO 2 injection, the maximum porosity loss has been increased to 0.7% in Case #85 and 2.5% in Case #87, from 0.3% and 0.9%, respectively. However, the CO 2 flow patterns are still almost identical for all three cases (figure not shown). Therefore, the impact of mineral reactions on CO 2 flow remains insignificant during the post-EOR CO 2 injection period. Figure 9 presents the estimated porosity loss at the end of the 600-year simulation period. The simulation results of using all low RSA values present almost no change in porosity due to mineral reactions (Figure 9a). On the other hand, using all mid and high RSA values leads to greater porosity reduction, as shown in Figure 9b,c. While there is only a slight expansion of porosity loss areas during the no-injection period, the maximum porosity loss due to mineral reactions increased to 1.19% and 5.04% for Case #85 and Case #87, respectively. It is worth noting that the predicted porosity changes in Figures 6c and 8c are more profound than those in Figure 9a,b, suggesting that using high RSA values leads to dramatically different porosity change predictions, hence mineral trapping of CO 2 , over even a short time period. Comparing the distributions of CO 2 global mole fraction ( Figure 10) and porosity change (Figure 9), it is clear that the mineral reactions had a nominal impact on the forecast of CO 2 migration in 600 years. It is interesting to note that there are lower porosity changes in the centers of the CO 2 plumes, where they exhibit very high (greater than 0.8) CO 2 global mole fractions (see Figures 7 and 10). The areas with greater porosity reduction are located on the edges (or fronts) of the CO 2 plumes, suggesting that there are more intense geochemical reactions and thus more CO 2 trapped in minerals or aqueous ions. This is because mineral reactions require sufficient contacts between the minerals and the formation fluids, which is less likely to be present in areas with very high CO 2 saturation.  The performance of trapping mechanisms that are directly related to the geochemical reactions was evaluated, as shown in Figure 11. As expected, and similar to the simplified model result (Figure 5), there is a significant difference between the amount of CO 2 trapped in minerals in the FWU field-scale reservoir model (Figure 11a). The mineral trapping first appeared as early as about 200 days (for Case #87), accelerated during the post-EOR CO 2 injection period (between the two vertical blue lines), and kept growing at a slower rate after CO 2 injection stopped. At the end of 600 years, the estimated amounts of CO 2 trapped in minerals are 3.25 × 10 6 kmol (1.43 × 10 5 metric tons), 0.8 × 10 6 kmol (3.52 × 10 4 metric tons), and 0.05 × 10 6 kmol (2.2 × 10 3 metric tons) for Case #87, Case #85, and Case #86, respectively. In other words, mineral trapping with all high RSA values is about four times more effective than with all mid RSA values, and 65 times more effective than with all low RSA values. At the end of 600 years, mineral trapping would contribute to 0.15%, 2.46%, and 9.44% of the total sequestered CO 2 at the FWU when using the low, mid, and high mineral RSA values, respectively. However, the mineral RSA values have much less impact on the CO 2 trapped in aqueous ions (Figure 11b). The maximum difference between predictions is only about 0.02 × 10 6 kmol (880 metric tons), and only 0.12 × 10 6 kmol (5.28 × 10 3 metric tons) of CO 2 was sequestered in aqueous ions by the end of the simulation. Nevertheless, mineral reactions and aqueous ions are able to sequester at least around 3000 metric tons of CO 2 (in Case #86), which would be otherwise presented in other forms (i.e., supercritical phase or dissolved) if reactive transport was not taken into consideration.

Discussion
Geochemical reactions, especially mineral reactions, are critical in predicting the fate and behavior of injected CO 2 into the subsurface and are involved in many aspects of a GCS project, such as the risk assessment of CO 2 leakage and caprock integrity. While most reactive transport studies were focused on CO 2 storage in deep saline aquifers, only a few incorporated the geochemical module into their simulations for CO 2 storage in depleted oil fields. Possible reasons for this could include: (1) mineral reactions in a depleted oil field are assumed to be similar to those in a deep saline aquifer, but in a much-restricted approach because of the presence of hydrocarbons; (2) the mineral trapping at depleted oil fields is considered to be less effective than other trapping mechanisms in a short period of time, thus it can be ignored for the purpose of storage estimation; and (3) most CO 2 storage numerical simulation packages are incapable of simulating reactive transport in three-phase systems (oil/water/CO 2 ). The results of this work demonstrate that neither the porosity reduction due to mineral reactions nor the mineral trapping of CO 2 at a CO 2 -EOR field is negligible. Using the mid RSA values of all seven minerals, up to 0.7% of porosity loss was forecasted after 20 years of CO 2 injection, and up to 1.19% of porosity loss was forecasted after 600 years. This is similar to the porosity change estimated by , which reported up to 0.7% of porosity reduction and up to 2.7% of porosity increase after 1000 years due to mineral reactions [45]. In that work, a modified FWU model was used, where there was only one five-spot well pattern and a deep saline storage scenario had been adopted (i.e., no hydrocarbon in the model). Nevertheless, similar results in estimated porosity change suggest that mineral reactions in three-phase systems probably are not as limited as it is assumed.
We have performed a preliminary study comparing the differences of mineral reactions in two-phase ( Figures A2 and A3) and three-phase systems (Figures 3 and 4) using the simplified 3D model. When hydrocarbon components were replaced with water, different predictions were found in most minerals. While the changes of illite and quartz slowed down in the monitoring period and tended to recover to the initial values in the threephase models, the dissolution of illite and precipitation of quartz in the two-phase model was accelerated. Moderate differences were observed in the carbonate minerals. In the two-phase models, an equilibrium state was reached during the first year of CO 2 injection and maintained throughout the entire simulation period ( Figure A3), leading to almost no CO 2 being trapped in minerals ( Figure 12); however, the three-phase models predicted that more carbonate minerals would precipitate in the monitoring period (Figure 4), which would then result in a significant amount of CO 2 being present in the minerals ( Figure 5). There is a clear impact of the presence of hydrocarbons on mineral reactions with this simplified model. Therefore, the assumption that the mineral reactions in three-phase systems are only a scaled-down version of those in two-phase systems is not necessarily correct. When considering the interactions between hydrocarbon components and minerals and aqueous species (which are not incorporated in this work), it would be inevitable to develop a new paradigm for the reactive transport modeling of CO 2 storage in depleted oil fields. In order to identify and quantify the roles of mineral RSA, many other important factors were fixed in this work. The simulation results will be affected by varying these factors, which include geological and hydrogeological properties (e.g., porosity and permeability values), multiphase flow models (e.g., relative permeability, hysteresis, and capillary pressure models), and operational factors (e.g., well spacing and depths, well operational limits, and CO 2 -EOR and WAG configurations). The impacts of mineral RSA will also be subject to other site-specific conditions when transferring the specific findings of this work to other case studies. However, it is believed that the importance of mineral RSA and the necessity of adding reactive transport modeling to CO 2 -EOR simulations are applicable to other GCS projects.
Additionally, reactive transport modeling is an indispensable component in understanding and mastering the coupled thermal-hydro-mechanical-chemical (THMC) processes in the subsurface porous media [46][47][48]. A thorough investigation of mineral RSA is vital for coupling mineral reactions with other critical subsurface phenomena, for example, mechanical and thermal deformation. On the other hand, the thermal and mechanical in situ conditions will affect mineral reactions on the whole. Given the capability of existing simulation tools, improvements in geochemical modeling and a better coupling of THMC processes are in urgent need. The geochemical properties should be treated the same as the flow or reservoir properties (e.g., porosity, permeability, pressure, saturation, etc.). The heterogeneity in mineral distributions and the temporal evolution of mineral properties should be incorporated into the reactive transport modeling. This work employed 87 combinations of mineral RSA values as a workaround to take both spatial and temporal variations in mineral RSA into account. While this approach has been effective in showing the impact of mineral RSA values on mineral reactions, it does not demonstrate the geochemical-hydrological coupling effect primarily due to the fact that the uniform distribution of mineral compositions across the model inhibits the ability of mineral reactions to affect the fluid flow. With the recent advances in RSA measurement techniques, it is vital to make better use of lab observations in numerical simulations for more accurate predictions. The absence of the geochemical-hydrological coupling effect could also be attributed to the spatial resolution (200 ft, or 60.96 m) of the numerical model employed in this work, which might be too coarse to capture reactions and fluid flow at small-scales (from centimeter to meter). Even with paralleled high-performance computing nodes, each of the 600-year simulations on the model with 200-ft resolution took more than eight hours. Refining the spatial resolution will no doubt increase the computational cost significantly, and very likely cause numerical convergence issues as the reservoir heterogeneity will be more complicated after adding more small-scale local geological features. Moreover, obtaining a reasonable down-scaled heterogeneity renders its own challenges and will be tackled in our future research.

Conclusions
In this work, the impact of mineral reactive surface areas on the CO 2 storage forecast was evaluated for the GCS project in the Farnsworth Unit in northern Texas. In order to reduce the computational cost, a simplified 3D model was used to screen representative RSA combinations from 87 cases. Three cases were chosen for the FWU reactive transport simulation. The main conclusions are drawn as follows: (1) The inter-dependency effects of mineral RSA values are stronger in the silicate mineral reactions and almost not observed in the carbonate mineral reactions; (2) The impact of mineral RSA values on CO 2 mineral trapping, on the whole, is more complex than it is on individual geochemical reactions. However, the three selected cases (with all low, all mid, and all high mineral RSA values) are representative for predicting CO 2 trapped in minerals; (3) While the low RSA case predicted negligible porosity change and an insignificant amount of CO 2 mineral trapping for the FWU model, the mid and high RSA cases forecasted up to 1.19% and 5.04% of porosity reduction due to mineral reactions, and 2.46% and 9.44% of total CO 2 trapped in minerals by the end of the 600-year simulation, respectively; (4) The presence of hydrocarbons affects geochemical reactions and can lead to net CO 2 mineral trapping, whereas negative CO 2 mineral trapping is forecasted when hydrocarbons are removed from the system.     Figure A1. Relationship between chemical equilibrium constants and temperature based on the EQ3/6 database (blue diamonds) for determining the chemical equilibrium constants of seven mineral reactions [43].