Coupled Hydromechanical Modeling of Induced Seismicity From CO2 Injection in the Illinois Basin

Injection of CO2 for geologic carbon sequestration into deep sedimentary formations involves fluid pressure increases that engage hydromechanical processes that can cause seismicity by activation of existing faults. In this work, we use a coupled multiphase fluid flow and geomechanical simulator to model spatiotemporal fluid pressure and stress changes in order to study the poroelastic effect of CO2 injection on faults in crystalline basement rock below the injection zone. The seismicity rate along features interpreted to be basement faults is modeled using Dieterich's rate‐and‐state earthquake nucleation model. The methodology is applied to microseismicity detected during CO2 injection into the Mount Simon formation during the Illinois Basin—Decatur Project. The modeling accurately captures an observed reduction in seismicity rate when the injection in the second well was into a slightly shallower zone above the base of the Mount Simon formation. Moreover, the modeling shows that it is important to consider poroelastic stress changes, in addition to fluid pressure changes for accurately modeling of the observed seismicity rate.

With the recent surge in seismicity attributed to anthropic activities in United States midcontinent, post-mortem numerical modeling studies serve to increase understanding of the mechanisms underlying induced seismicity (Choy et al., 2016;Ellsworth et al., 2015). Most studies assume that fault reactivation is primarily driven by pore pressure diffusion and thus neglect injection-induced poroelastic stress changes. However, several recent numerical studies indicate that poroelastic effects must be captured by numerical models in order to correctly forecast fluid-induced seismicity (Barbour et al., 2017;Zhai et al., 2019). Most of these studies are related to wastewater injection or enhanced geothermal systems (Barbour et al., 2017;Hakimhashemi et al., 2014;Norbeck & Rubinstein, 2018), while only a handful of sites reported fluid-induced microseismicity (i.e., < 2 , not felt by humans) associated with CO 2 injection, namely at In Salah, Algeria (Rutqvist et al., 2016;Verdon et al., 2015); Otway, Australia (Myer & Daley, 2011;Siggins, 2010), and the Illinois Basin-Decatur Project (IBDP), United States (Bauer et al., 2016;Kaven et al., 2015;Will, El-Kaseeh, et al., 2016;Will, Smith, et al., 2016).
The IBDP is the first carbon capture and sequestration project in the United States that injected commercial volumes of CO 2 into a deep saline aquifer for GCS (Finley, 2014). One million tons of CO 2 was injected over a 3-year injection period from November 2011 to November 2014 at the well CCS1 into a high permeability Mount Simon Sandstone interval at a depth of around 2,140 m. Nearly 20,000 induced microseismic events were detected with most events located within the underlying crystalline Precambrian basement . Identified clusters of microseismic events form semilinear features oriented within 30° of the direction of the maximum horizontal principal stress (azimuth N068°) and indicate that the seismicity at the IBDP is occurring along pre-existing basement faults (Goertz-Allmann et al., 2017). In April 2017, CO 2 started to be injected in the CCS2 well in a zone less than 50 m shallower than the injection zone in the CCS1 well. Injection was into the Lower Mount Simon in both wells, but with a higher injection rate in the CCS2 well compared to that of CCS1. Yet, there is very little microseismicity occurring during injection into CCS2 .
Here, we demonstrate an approach for modeling the induced seismicity observed at the IBDP along basement faults using multiphase fluid flow and geomechanical model simulations coupled with a rate-and-state nucleation model. The organization of the study is as follows: In Section 2, we describe the computational model used in this study. In Section 3.1, we show the result of the earthquake catalog declustering that is used to calibrate the rateand-state parameters in Section 3.2. The modeling results for both injections in wells CCS1 and CCS2 are then detailed in Section 3.3. Finally, we use our model to study the effect of varying injection scenarios on predicted induced seismicity in Section 3.4.

Numerical Model
In this work, we consider a three-dimensional domain that includes discretized faults within the basement inferred from the microseismic clusters observed at the IBDP site. We simulate the CO 2 injection and stress evolution using the coupled multiphase flow and geomechanical model and apply the rate-and-state seismicity model to study the response of the basement faults to the CO 2 injection. LUU ET AL.

Computational Model
We consider a simplified version of the subsurface structure at the Decatur site with a three-dimensional layercake model geometry consisting of 10 homogeneous geological layers with the top layer representing the primary seal Eau Claire formation (1,540 m depth below ground surface [bgs]) and the bottom layer representing the crystalline basement (2,100-3,000 m bgs; Bauer et al., 2016). The Mount Simon sandstone formation is divided into six different layers in the model, representing from bottom to top, the Lower Mount Simon A-lower zone, Lower Mount Simon A-upper zone, and the Mount Simon B, C, D, and E zones. The injection interval is located in the Mount Simon A-lower zone which has been divided into three sublayers to improve flow modeling within the reservoir. A thin continuous mudstone layer is included to honor multilevel pressure data recorded at the IBDP site which shows that vertical migration of the CO 2 plume that formed after injection into CCS1 is limited by discontinuous low-permeability layers that inhibit vertical fluid flow within the reservoir (Senel et al., 2014;Strandli et al., 2014;Williams-Stroud et al., 2020). Hydromechanical properties of the geological layers are summarized in Table 1.
Sixteen microseismic clusters are identified using the DBSCAN algorithm (Ester et al., 1996) and used to map faults in our model. Faults are discretized as finite-thickness elements within the basement and are displayed in Figure 1 (middle). This type of fault model is conceptually similar to conduits with along fault flow dominated by flow in a highly fractured damage zone (Caine et al., 1996). All faults are about 20 m thick and uniform in properties (i.e., fault core is not distinguished from damage zone). A detailed microseismic analysis showed that the basement faults at Decatur are hydraulically connected to the reservoir (Goertz-Allmann et al., 2017). Therefore, the faults discretized in our model vertically extend from the bottom of the reservoir (2,146 m) to the bottom of the model (3,000 m). We consider the faults to be hydraulically conductive with permeability logarithmically decreasing with depth from 1 mD at the top (2,146 m) to 0.1 mD at the bottom (3,000 m). This type of permeability variation has been reported to be associated with critically stressed crystalline basement faults (Barbour et al., 2017;Townend & Zoback, 2000). Mechanically, faults are assumed to be transversely anisotropic where the Young's modulus in the direction normal to the plane of isotropy is equal to 80% of that of the plane of isotropy (Glamheden et al., 2007). As the locations of induced seismicity do not clearly outline the faults and for the sake of simplicity, we consider that all faults are vertical (dip angle = 90 • ) and have the same elastic properties as the host rock units they transect (i.e., only the permeability of faults are different).
Figure 1 (left-hand side) shows the computational mesh with the applied boundary conditions. The mesh consists of 200 × 200 × 50 (2 million) hexahedral elements uniformly discretized horizontally and refined vertically in the vicinity of the injection zone. Lateral and bottom boundaries are open to fluid flow with only the top boundary being closed to flow. We apply fixed stress conditions at lateral and top boundaries, and rollers at the bottom (no vertical displacement). Following Senel et al. (2014), we assume an initial hydrostatic gradient for pore pressure (10.15 MPa/km) and vertical geothermal gradient for temperature (18.2°C/km). The system is initially  100% brine-saturated with salinity of 20% and hydrostatic initial fluid pressure. Initial in situ stress conditions are defined according to Bauer et al. (2016) and correspond to a strike-slip faulting system with > > ℎ (Figure 1, right-hand side). In situ stress measurements show that the maximum horizontal stress direction has a fairly constant azimuth and is oriented N068° (Bauer et al., 2016;Williams-Stroud et al., 2020). The minimum horizontal stress gradient in each formation is estimated based on measurements (obtained using hydraulic fracturing, overcoring and borehole pressure meter methods), whereas the maximum horizontal stress gradient is calculated assuming that the host rock is near critically stressed conditions for instability for a friction coefficient = 0.6 . Figure 2 shows capillary pressure and relative permeability curves used in the multiphase fluid flow simulation. We follow Mehnert et al. (2019) and use van Genuchten capillary pressure model (Genuchten, 1980) with fitting parameter = 0.55 , residual liquid saturation = 0.6 , saturated liquid content = 0.999 and maximum capillary pressure = 6.9 MPa. Relative permeability curves are constructed using the van Genuchten-Mualem model with fitting parameter = 1.36 , residual liquid saturation = 0.65 and residual gas saturation = 0.01 .
For simplicity, all coordinates shown in this work are relative to injection well CCS1 (i.e., well CCS1 is located at = = 0 m and its top is at = 0 m).

Rate-and-State Seismicity Model
We model the seismicity rate using a hybrid approach where seismicity rate is calculated from time-dependent pressure and stress changes simulated by our coupled hydromechanical model (Hakimhashemi et al., 2014). We use Dieterich's rate-and-state earthquake nucleation model to assess the evolution of seismicity rate due to injection-induced stress changes along basement faults (Dieterich, 1994). The rate-and-state seismicity model estimates the number of independent events in response to a change in stress on a set of faults and is described by the following ordinary differential equation (Dieterich, 1994;Segall & Lu, 2015) = ̇ ̇ 0 − (1)  where is the ratio between the seismicity rate relative to the background rate, =̇0 is the characteristic relaxation time, and 0 are the Coulomb and background stressing rates, respectively. We solve the ordinary differential equation using a fifth-order adaptive time step Runge-Kutta-Fehlberg algorithm (Fehlberg, 1969) with a relative error tolerance = 10 −6 .
We note that the rate-and-state seismicity model is only applicable if optimally oriented faults are already critically stressed prior to injection (Chang & Segall, 2016;Zhai et al., 2019). In addition, the theory only relates to earthquake nucleations (mainshocks) and does not account for the physical processes involved in aftershock sequences. More specifically, while the geomechanical model accounts for stress transfer from injection pressure changes and poroelastic stress propagating ahead of the pressure front, stress changes induced by seismic slip of individual fractures or faults that can trigger another event are not included. Thus, earthquake catalogs must be declustered (removal of aftershocks) to be able to compare observed seismicity rates with results of the rate-andstate model. This model limitation also implies that it does not forecast magnitudes of earthquakes. However, physics-based seismicity rate models can be combined with the Gutenberg-Richter law to calculate the probability of occurrence of an earthquake of magnitude (Navas-Portella et al., 2020;Segall & Lu, 2015). For magnitudes , the total number of events at location = ( ) at time step is defined as with 0 and being the background seismicity rate and the b-value, respectively. The total number of earthquakes of magnitude at time step is then obtained by integrating Equation 2 over the volume : The number of earthquakes in time interval [ 1,2] is written where and are the minimum and maximum magnitudes simulated. In the following, is set to the catalog's magnitude of completeness and is chosen sufficiently large.
Assuming that earthquake occurrence is described by a inhomogeneous Poissonian process, Zhai et al. (2020) estimates the magnitude probability of exceedance (i.e., the probability of having at least one event of magnitude larger than ) following where ≥ ( 1, 2) is the expected number of earthquakes with magnitude greater than or equal to . Let us define the following cumulative probability distribution as a function of earthquake magnitude: A magnitude-time distribution can be simulated by randomly sampling ( , +1) earthquakes over this distribution for the whole injection period (Zhai et al., 2020).

Stressing Rate Modeling
Dieterich's rate-and-state seismicity model relates changes in Coulomb stress to changes in seismicity rate. We define the Coulomb stressing rate as the change in Coulomb stress ΔCFS per unit of time which is calculated at each time step of the simulation following LUU ET AL.

10.1029/2021JB023496
6 of 19 where is the friction coefficient (assumed to be 0.6 for all faults), Δ is the change in shear stress, Δ is the change in normal stress (positive for tension), and Δ is the change in fluid pressure. Shear stress and normal stress acting on a fault plane can be calculated from the stress tensor following where is the normal vector of a given fault plane and ||⋅|| denotes the Euclidean norm.
We simulate the spatiotemporal distributions of fluid pressure, shear, and normal stresses using the latest version of the coupled fluid flow and geomechanical software TOUGH-FLAC (Rutqvist, 2011;Rutqvist et al., 2002) that sequentially couples the finite-volume multiphase flow simulator TOUGH3 (Jung et al., 2017) and the commercial finite-difference geomechanical software FLAC3D V7. The latest version of TOUGH-FLAC (Rinaldi et al., 2021) integrates all the new features of TOUGH3, in particular the use of PETSc parallel solvers which allows execution of coupled simulations with a large number of grid blocks (here, 2 million elements). By the use of TOUGH-FLAC, we account for full hydromechanical coupling with porosity changes modeled as a function of bulk modulus and volumetric strain (Kim et al., 2011). Fluid pressure and stresses are calculated at discrete time steps controlled by TOUGH3 using adaptive time stepping based on the number of Newton-Raphson iterations needed for each time step. However, we set the maximum time-step size to 3 days to better capture amplitudes of pressure changes due to the numerous shut-in phases. We further fit cubic splines to the simulated pressures and stresses which are used to calculate the changes in Coulomb stress ΔCFS . Finally, the stressing rate is taken as the numerical time derivative of ΔCFS with a time step size = 1 day, followinġ The coupled hydromechanical model generates spatial and temporal distributions of pressure and stress in the whole model. However, we assume that seismicity occurs only along pre-existing critically stressed faults and therefore only calculate Coulomb stress changes at integration points corresponding to the finite-thickness fault elements.

Results
We apply our coupled hydromechanical rate-and-state nucleation model to generate seismicity forecasts for both injections in wells CCS1 and CCS2. Modeled seismicity, especially for the first injection during which most of the seismicity is observed, is compared to the declustered catalog.

Catalog Declustering
At the IBDP, more than 5,000 microseismic events have been located with magnitudes ranging from −2.1 to 1.2, and the magnitude of completeness is = −0.7 (Goertz-Allmann et al., 2017; Williams-Stroud et al., 2020). Earthquake catalogs usually contain independent earthquakes (mainshocks) and earthquakes resulting from stress release after a mainshock (aftershocks). As explained in Section 2.2.2, the rate-and-state theory mainly focuses on the mainshocks and therefore aftershocks must be removed from the catalog to compare the hybrid seismicity model with observed seismicity. We consider a complete catalog and remove events with magnitudes lower than the magnitude of completeness ( −0.7 ). We decluster the earthquake catalog using the nearest-neighbor method (Zaliapin & Ben-Zion, 2020). By this approach, for each earthquake in the catalog, we calculate the nearest-neighbor interevent distance in the space-time-magnitude domain. Given a pair of events and , the nearest-neighbor is calculated following = where and are the rescaled time and distance, respectively, written with the Euclidean interevent distance, the interevent time, the fractal dimension of earthquake epicenter/hypocenter, a weighting coefficient, and the magnitude of event . In this study, we did not consider earthquakes' depths due to uncertainties in their locations. In Equation 11, only and are user-defined while the other terms depend on earthquakes' parameters. Besides, the nearest-neighbor method has two additional parameters, namely an initial cutoff threshold 0 and a cluster threshold 0 . A sensitivity analysis of these four parameters (see Section S1 in Supporting Information S1) showed that the declustering is mainly sensitive to the cluster threshold 0 . Figure 3 shows the 2D distributions of calculated nearest-neighbor distance for the full catalog (left-hand side) and a declustered catalog (right-hand side). An earthquake that yields a low distance is close in space and time to its nearest-neighbor and is thus discriminated as an aftershock.
Following Zaliapin and Ben-Zion (2020), we set = 1.1 (b-value) and = 1.5 (epicenter), we determine 0 using the bimodality of the distribution of earthquake nearest-neighbor proximities, and 0 is initialized to 0 and expanded until the clustered mode is removed (i.e., lower left corner of plots in Figure 3). Thus, we finally set 0 = 10 −4.5 and 0 = 1.9 . The declustered catalog is plotted against the observed catalog in Figure 4 and indicates that most of the clustered events occurred at the beginning of the injection around March and July 2012.
It should be mentioned that declustering algorithms are usually tailored to remove aftershocks in natural earthquake catalogs. Nevertheless, the bimodality of the nearest-neighbor distribution has been observed in induced earthquake sequences (Schoenball et al., 2015;Schoenball & Ellsworth, 2017;Zaliapin & Ben-Zion, 2016).

Parameter Calibration
The rate-and-state seismicity model is governed by three parameters, namely the background stressing rate, the background seismicity rate and a constitutive parameter that controls the characteristic relaxation time. The background stressing rate is usually obtained through geodetic measurements and the background seismicity rate can be estimated by monitoring the seismicity prior to the injection. We use a background stressing rate of 0 = 5 Pa/year as estimated for the Southern Illinois Basin (Hamburger et al., 2010). Continuous microseismic monitoring has been carried out at the IBDP site prior to the first injection during 18 months and eight earthquakes with magnitude −1.5 were interpreted as local events (Smith & Jaques, 2016). Because of the lack of recorded natural earthquakes with magnitude ≥ −0.7 in the area of study, we calibrate the background seismicity rate 0 along with the parameter by manually fitting the modeled cumulative number of events during the first injection to the observed one (Hakimhashemi et al., 2014). We use a global optimization algorithm, namely the CMA-ES (Hansen & Ostermeier, 2001), to further refine the two parameters with a population size of 20  for the evolutionary algorithm, a maximum of 100 iterations and the manually calibrated parameters as initial mean. Eventually, we found a background seismicity rate 0 = 0.35 events/year ( ≥ −0.7) and = 0.0295 MPa. It should be mentioned that the background seismicity rate is estimated so that model outputs (relative to background seismicity) can be compared with the declustered catalog. The calibrated value of the background seismicity rate is fairly uncertain (see Section S2 in Supporting Information S1). Figure 5 shows the calibration result displayed against the observed cumulative number of events. Figure 6 shows the modeled seismicity rate (integrated over all integration points) considering both injections in wells CCS1 (from November 2011 to November 2014) and CCS2 (from April 2017 to April 2018). To model the first injection, we inject CO 2 below a low permeability mudstone layer simplified in the model to represent the discontinuous baffles in the Mount Simon that restrict vertical flow. The two perforated zones (2,121 m and 2,129 m bgs) are modeled as single injection element. For the second injection, CO 2 is injected above the low permeability layer 50 m shallower (2,178 m bgs) compared to the first injection. We note that only the first year of the second injection is modeled. Overall, the modeled seismicity rate follows the average behavior of the observed seismicity rate. More specifically, the modeled seismicity rate is consistent with the declustered catalog in terms of onset timings and peak rate amplitudes for the first injection, which means that the model is able to reproduce the main temporal features of the earthquake sequence. We observe that many of the longest shut-in phases (e.g., September 2012, March 2013, February 2014, October 2014) yield a sharp decrease in the modeled seismicity rate which indicates that the modeled seismicity rate and the injection rate are correlated. After the end of the injection in well CCS1, the modeled seismicity rate progressively decreases and predicts a lower rate than the background seismicity from July 2015 due to negative stressing rates. Despite larger injection rates in well CCS2 (1.7 times the injection rate in well CCS1), the modeled seismicity rate is negligible compared to the seismicity induced by the first injection (about two orders of magnitude smaller in seismicity rate).

of 19
We simulate an earthquake magnitude-time distribution by combining the Gutenberg-Richter law and the injection-induced seismicity rate obtained using our coupled hydromechanical model (Equation 6), Figure 7 displays the simulated distribution compared to the observed one, and shows that the simulated earthquake magnitude time dependence is fairly consistent with the observations. We further investigate the relative contributions of pressure and injection-induced stress changes to the Coulomb stress by tracking their evolutions at different points throughout the model. According to Equation 7, contribution to the Coulomb stress changes are changes in shear stress Δ , changes in normal stress Δ and changes in fluid pressure Δ . The changes in shear stress and normal stress are induced by poroelastic stresses in the system that are in turn due to injection-induced pressure changes in the system. In the following, we define Δ as the pressure contribution and terms Δ + Δ of Equation 7 as poroelastic contribution to the Coulomb stress change. Positive ΔCFS indicates weakening of the fault planes bringing them closer to failure. To monitor the evolution of the pressure and poroelastic stress to the Coulomb stress, we select a first point midway between wells CCS1 and CCS2 to study the near-field, and a second point on the westernmost fault for the far-field. For both points, we display in Figure 8 the evolution of the Coulomb stress change (black) along with the contribution from pore pressure (purple) and poroelastic stress (green) changes at different depths (2,200 m, 2,400 m, and 2,600 m bgs).
During the first injection at well CCS1 (before January 2015) for both the near-field and the far-field cases, the pressure front expands radially around the wellbore until it reaches the high permeability faults causing downward propagation of the pressure front within the basement. Besides, poroelastic stresses have a strengthening effect on the fault planes where their contributions are negative (Δ + Δ < 0 ). Inversely, during the second injection at well CCS2 (after April 2017), we observe that Δ + Δ > 0 which indicates that fault planes are weakened by poroelastic stresses. Indeed, in the first case, the fluid is injected below the impermeable baffle, yet right above the basement where the permeable faults are hydraulically connected to the reservoir. When the fluid flows into the faults, it causes the pressure to build up within resulting in compressive stress. On the other hand, during the second injection, fluid is injected above the baffle layer disconnecting the faults from the zone of injection. In this case, pressure increase within the reservoir causes lateral expansion of the rock formation, inducing extensive stress within the basement beneath. At the top of the basement in the near-field (upper left plot), poroelastic stress effects are not negligible and impede reactivation by reducing the changes in Coulomb stress as ΔCFS < Δ . However, farther from the reservoir, we have Δ + Δ ≈ 0 and ΔCFS ≈ Δ , which indicates that the poroelastic stress impact decreases with depth where direct pressure effects become dominant.
The modeled relative seismicity rates at the selected points are also shown in Figure 8 on a logarithmic scale (blue). We note the exponential relationship between the Coulomb stress and the relative seismicity rate, consistent with the solution to the ODE described by Equation 1 when the stressing rate is larger than the seismicity rate (i.e., ̇̇0 ≫ ). We also consider a case where we neglect the first injection and only model the second injection (black dashed line). For this case, the injection starts at the original reservoir pressure. This is to investigate the relevance of the stressing history on the seismic response. Looking at the modeled seismicity rate (right vertical axis of Figure 8), we observe that in the near-field, it becomes lower than the estimated background seismicity ( 1 ) after the end of the first injection. This is likely due to the post-injection pressure drop yielding a negative pressure rate, and thus a decrease in the relative seismicity rate (Almakari et al., 2019; see Section S4 in Supporting Information S1). Once the pressure rate becomes zero (around July 2015), the relative seismicity rate steadily increases back. This behavior is not observed in the far-field where the seismicity rate goes back to the 10 of 19 estimated initial background value ( = 1 at the end of the first injection). Interestingly, the relative seismicity rate in the near-field at the top of the fault is about one order of magnitude larger for the second injection if we neglected the first one. Nevertheless, while this local decrease in seismicity rate following the shut-in of the first injection may have contributed to the lack of recorded seismicity during the second one, its impact is negligible compared to the overall lower pressure changes acting on the faults.

Effect of Injection Rate on Seismicity Rate
For equivalent total injected volume, Barbour et al. (2017) showed that a variable injection rate may induce more seismicity compared to a constant injection rate. Here, we investigate the effect of four different injection scenarios by comparing the seismicity rate generated by our hydromechanical earthquake nucleation model for different injection rates. We note that we only simulate and compare with the injection in well CCS1 for which fluid-induced seismicity has been observed. More specifically, given a total volume of approximately 1 million tons injected within the span of 3 years, we consider a first constant injection rate at 11 kg CO 2 /s (Scenario A), and a second piecewise constant injection rate increasing from 10 kg CO 2 /s to 12 kg CO 2 /s (Scenario B). For both rates, in Scenarios C and D, we also consider a variant with 2-week shut-in phases every 6 months (equivalent to the longest shut-in period during injection in CCS1). Figure 9 shows the seismicity rates and annual magnitude probabilities of exceedance for the four injection rates using a b-value of = 1.1 (Bauer et al., 2016), and a minimum and maximum magnitudes of = −0.7 and = 4 . Figure 10 shows the same annual magnitude probabilities of exceedance as Figure 9 with the results displayed for each year. The seismicity rate and magnitude of exceedance probability modeled for the first injection in well CCS1 are also displayed for comparison. For the reference case, the probability for exceeding M2 is 24%, 21%, and 18% in 2012, 2013, and 2014, respectively. For Scenario A, most seismicity occurs at the beginning of the injection and decreases over time resulting in a higher probability for exceeding M2 in 2012 (32%). For Scenario B, the seismicity rate steadily increases up to a maximum of 0.9 events/day and followed by a steady decrease, annual probabilities for exceeding M2 are simi- 12 of 19 lar throughout the injection (about 22%). Scenarios C and D show that the shut-in phases induce an immediate drop in the seismicity rate, followed by a larger seismicity rate increase when the injection restarts, compared to scenarios A and B. This behavior is also observed in the reference case where long shut-in phases (e.g., March 2013, February 2014, October 2014) yield an instantaneous drop in seismicity rate which subsequently increases with a time lag. Yet, annual probabilities for exceeding M2 are only slightly lower for both Scenarios C and D.
Overall, in our model, we observe a correlation between the injection rate and the modeled seismicity rate for which the response appears to depend on the amplitude variations of the injection rate.

Discussion and Conclusion
In this study, we modeled the seismicity induced along the Precambrian basement faults by the two CO 2 injection wells at Decatur Illinois, specifically wells CCS1 and CCS2 from November 2011 to April 2018. Our coupled hydromechanical model reproduces characteristic features of the observed microseismic activity. The modeled seismicity rates are comparable to recorded seismicity in terms of onset timings and peak rate amplitudes for the first injection, while modeled seismicity is negligible for the second injection consistent with field observations. Our modeling results suggest that the seismicity at Decatur is strongly influenced by pressure effects. However, modeling of injection in CCS1 indicates that poroelastic stresses are not negligible and tend to impede reactivation, in particular in the vicinity of the injection wells. Because the seismicity rates forecast by the rate-and-state earthquake nucleation model are exponentially related to the pressure and poroelastic stress rates, ignoring poroelastic effects (i.e., only flow modeling and using ΔCFS = Δ ) would overpredict the seismicity rate by approximately one order of magnitude according to our model (using the same parameters for rate-and-state simulation). This result highlights the necessity of coupled hydromechanical modeling to accurately capture the main physical processes related to fluid-induced seismicity, in agreement with recent studies (Barbour et al., 2017;Chang & Segall, 2016;Fan et al., 2019;Zhai et al., 2019).
In addition, the rate-and-state model used in this study estimates the induced seismicity rate relatively to the background seismicity rate. Due to the lack of earthquake with magnitude ≥ −0.7 recorded prior to the injection, we could not derive a background seismicity rate based on recorded data. Hence, we calibrated the background seismicity rate along with the constitutive parameter to match the observed seismicity rate. The background seismicity rate has been calibrated for comparison purpose and the inverted value is fairly uncertain (see Section S2 in Supporting Information S1). An analysis of the sensitivity of the seismicity rate with respect to the background seismicity rate shows that modeled seismicity rates would fit the observed seismicity rate comparably well for background seismicity rates ranging between 0.2 and 0.6 events/year. Outside this range, the main peak amplitudes of the earthquake sequence are not properly captured. Given these uncertainties on the background seismicity rate, the modeled seismicity rates shown in this work must be interpreted within the context of a probabilistic analysis (Barbour et al., 2017). Nevertheless, regardless of the value of the background seismicity rate, our numerical model shows that the seismicity rate after the end of the first injection becomes lower than the initial background seismicity rate ( 1 ), in particular near the injection well CCS1. Similar results have been observed in other modeling studies (Almakari et al., 2019) and is linked to negative pressure rates as the pressure is diffusing out of the faults. Due to the stressing history, seismicity rates forecast for the second injection are lower than if we had ignored the first injection phase. Nevertheless, despite the higher injection rate, the modeled pressure changes induced by the second injection on the basement faults are significantly lower which indicates that the absence of observed seismicity during the second phase is principally due to the injection zone in the CCS2 well location above the low permeability mudstone layer and the higher porosity and permeability in CCS2 injection zone relative to the CCS1 injection zone .
We note that our model domain is fairly simple and consists of a three-dimensional layer-cake model that only includes vertical basement faults inferred from the observed microseismic clusters. Structural faults interpreted in the 3D seismic volume that could potentially impede pressure diffusion have not been modeled. Because of the locations of the induced seismicity do not tightly outline the faults, and also due to the absence of inverted source mechanisms for all the events in the catalog, we assumed that all faults are vertical and that micro-fault planes are aligned with the embedding faults. This simplification can lead to errors in the calculation of the CFS using Equation 7 and thus in the resulting seismicity rate. Langet et al. (2020) showed for one cluster in the catalog that the source mechanisms are mostly dipping sub-vertically with slip direction around −10°. However, source mechanism inversion should be carried out for the whole catalog to obtain more accurate CFS and seismicity rate. For the sake of simplicity, we also assumed that the basement faults are merely hydraulically connected to the lower part of the reservoir, but do not vertically extend across it. Additionally, we only considered a homogeneous set of faults with invariable permeability. Several numerical simulations have demonstrated that location and timing of fluid-induced seismicity is affected by the variations of fault permeability (Chang & Segall, 2016;Zhang et al., 2013). Besides, we considered the rate-and-state parameters to also be homogeneous across the area of study. All these simplifications yield some discrepancies between the model outputs and the observations. For example, in our model, seismicity initiates at the top of the faults and propagates downward into the basement, while in some clusters, the observed seismicity starts within the basement. Despite these disagreements and the low complexity of our model, it is able to reproduce the principal features of the earthquake sequence recorded at the site, implying that the main physical processes involved are captured by our model.
Finally, we used our coupled hydromechanical earthquake nucleation model to study the effect of different injection scenarios on the seismicity rate, assuming equivalent total injected volumes of CO 2 . We found a correlation between the injection rate and the modeled seismicity rate. More precisely, seismicity rate immediately decreases in response to a shut-in phase and increases with a time lag when the injection restarts, the peak amplitude of the seismicity rate depending on the amplitude of the injection rate increase. However, we did not find significant changes in terms of modeled seismicity (total number of events and magnitude probability of exceedance) between the few scenarios tested and the actual injection rate, which can be explained by the already fairly constant injection rate used for the first injection in well CCS1. Additional studies are being planned to improve the model by considering heterogeneity in several model parameters and to identify factors leading to more accurate characterization of the risk of inducing earthquakes in GCS activities.

Appendix A: Calibration of Hydromechanical Model Parameters
We history-matched hydrological model parameters against multilevel pressure data and saturation profiles measured at verification well VW1 located approximately 300 m away from injection well CCS1 (see Figure 1, middle). We inverted porosities and permeabilities of layers Mount Simon A-upper through Argenta by minimizing the joint objective function defined by Equation A1, written where is the vector of model parameters to invert (porosities and permeabilities of layers), subscripts and respectively denote pressure and saturation, obs and obs are the measured data vectors to history-match, and ( ) and ( ) the data vectors calculated by the forward operator . and are coefficients that weigh the contributions of each data set to the joint objective function, and are arbitrarily set to 1 and 2, respectively (with pressure expressed in MPa). The objective function is optimized using the CMA-ES (Hansen & Ostermeier, 2001) which is known to be a robust stochastic global optimization algorithm, especially when the number of parameters to invert is relatively high (Auger, 2016). For the CMA-ES, we use a population size of 20 and 100 iterations, the initial means and standard deviations are summarized in Table A1  To reduce the computational cost of the forward modeling (TOUGH3 simulation), we considered a radially symmetric layer-cake model with the same layering as our 3D computational mesh. Only the pressure data measured in the vicinity of the injection zone (zones 1 through 4) and the saturation profiles measured in March and July 2012 are inverted. Results of the history matching for the best fit model are represented in Figure A1. Poisson's ratio, bulk modulus, Biot's coefficient and pore compressibility are calculated using mechanical conversion functions or empirical models (see Table A2). Figure A2 shows the pressure change front and CO 2 plume modeled for March 2012 using the 3D geomechanical model with porosity and permeability values inverted using the radial layered mesh.  (Settari et al., 2005) Note.Young's moduli are given in Table 1.

Data Availability Statement
The induced seismicity catalog (2019 version) and injection data used in this work were acquired by the Illinois State Geological Survey under projects funded by the U.S. Department of Energy through the National Energy Technology Laboratory. Updated data sets have been uploaded to EDX (https://edx.netl.doe.gov/dataset). The three-dimensional computational mesh is generated using the open-source meshing software LaGriT (https:// lagrit.lanl.gov/). Hydromechanical properties are taken from published literature and the computational model is fully described in Section 2.1. The numerical simulations are carried out using TOUGH3-FLAC3D. TOUGH3 is a fluid-flow numerical simulator developed at Lawrence Berkeley National Laboratory and FLAC3D is a geomechanical simulator commercialized by Itasca Inc. TOUGH3 input and output simulation files are pre-and post-processed using the Python package toughio (Luu, 2020(Luu, , 2022b. The calibration of the model parameters (history matching and rate-and-state) uses the CMA-ES optimizer implemented in the Python package stochopy (Luu, 2021). Earthquake declustering, modeling of rate-and-state seismicity and magnitude-time distributions are implemented in the Python package bruces (Luu, 2022a). These three packages are available on Zenodo at the links below: (a) toughio: https://zenodo.org/record/3961278; (b) stochopy: https://zenodo.org/record/4058008; (c) bruces: https://zenodo.org/record/6422572.