Was the 2018 Hokkaido Eastern Iburi Earthquake in Japan Related to Carbon Capture and Storage (CCS)-CO2 Injection? Insights From Geomechanical Analysis

In a recent paper [“Groundwater anomaly related to CCS-CO2 injection and the 2018 Hokkaido Eastern Iburi earthquake in Japan” by Sano et al. (Front. Earth Sci., 2020, 8)], the authors claimed that CO2-enriched fluid may have initially migrated through permeable channels, blocking the fluid flow from the source region, increasing pore pressure in the focal region and triggering a natural earthquake where the brittle crust was already critically stressed. The proposed model is very interesting, but the authors have not shown any quantitative evaluation supporting their conclusion. Here, through geomechanics model analysis, even under extreme conditions, which overestimate the impact of the injection, the impact of the CO2 injection on the Iburi earthquake fault, whether the deep section or shallow part of the fault, is much lower than that caused by Earth tides. In addition, no convincing mechanism exists that would allow fluid channels to heal within a short period of time and block the natural fluid flow along the fault. Therefore, the occurrence of earthquakes was not related to CO2 injection. Geological storage of CO2 is expected to become an effective option for global warming countermeasures, and the assessment of its environmental impact must be carefully conducted.


INTRODUCTION
As one effective option for mitigating global warming, the geological storage of carbon dioxide (CO 2 ) appears promising. Injecting a large amount of CO 2 fluid into an underground aquifer will inevitably cause certain formation deformation and increase fluid pressure. There is therefore a risk of accelerating the activation of pre-existing faults nearby and inducing/triggering earthquakes. In several industrial applications that require the injection of fluids into the Earth's crust, including disposal of wastewater, development of unconventional oil and gas resources, enhanced geothermal development, and well salt production, injection-induced seismicity, including destructive earthquakes with a sizable magnitude (M > 5), has been observed and has attracted widespread attention (e.g., Ellsworth et al., 2019;Atkinson et al., 2020;Lei et al., 2020). At the same time, it also poses challenges for effective operation. At present, some demonstration projects of CO 2 geological storage have been carried out around the world, but there are no reports of important induced earthquakes related to these projects.
The carbon capture and storage (CCS) demonstration project in Tomakomai, Hokkaido, Japan, is a medium-scale CO 2 offshore storage and has injected approximately 300,000 tons of CO 2 into a porous sandstone aquifer at a depth of approx. 1,000 m from April 2016 to November 2019. During this period, on September 5, 2018, at 18:08 (UTC), a JMA (the Japan Meteorological Agency) magnitude 6.7 earthquake occurred in the eastern Iburi region in Hokkaido. The focal depth of the earthquake is 37 km, and the epicenter is approx. 20 km east from the injection location ( Figure 1). Joint inversion of strong motion and geodetic data shows that a large slip occurred at a depth of approx. 22 km, which is shallower than the hypocenter (initial rupture) depth. This indicates that the main seismogenic fault is an east-dipping high-angle reverse fault, which is consistent with the Ishikari-Teichi-Toen fault zone (ITTFZ) (Kobayashi et al., 2019). Based on detailed 3D injection simulation, the predicted pressure buildup at the end of the injection is limited within a few kilometers from the injection well (Kano et al., 2014). By considering that the stress disturbance caused by CO 2 injection cannot reach the ITTFZ, it is generally believed that the earthquake is not directly or indirectly related to CO 2 injection. However, recently, a previous study (Sano et al., 2020) determined isotopes of commercially available mineral water produced from a depth of approx. 100 m in a well located 13 km from the injection site, showing dissolved components of injected CO 2 . Furthermore, the steady flow of deep celestial fluid along the ITTFZ is speculated to be blocked, which increased the fluid pressure in the deep part of the fault and triggered the earthquake. The model proposed in the paper is very interesting, but making such an assumption without any geomechanical analysis may not be valid.
The purpose of the present brief study is to investigate the validity of the above assumption. First, the theoretical solution of a simple two-dimensional radial flow model is used to analyze the decay of fluid pressure with distance. Then, a simplified threedimensional (3D) and coupled thermal-hydraulic-mechanical (THM) model is constructed to quantitatively analyze the influence range and degree of actual CO 2 injection under some extreme conditions.

ANALYSIS BASED ON GEOMECHANICS
Geological Environment of the Tomakomai Carbon Capture and Storage-Injection Site Figure 2 shows a simplified 3D geological model of the Tomakomai CCS-CO 2 injection site and distribution of earthquakes from the JMA catalog. The geological model is modified from the summary of time report of the Tomakomai CCS large-scale demonstration test 300,000 tons injection (METI, NEDO, and JCCS, 2020) METI, NEDO, and JCCS. There are two potential reservoirs, Moebetsu formation and Takinoue formation, for CO 2 storage. The Moebetsu formation is located at a depth of approx. 870 to 1,200 m. The upper part of the formation, which is dominated by silt/mudstone, acts as a cap layer. The lower part of the formation (below 1,070 m), FIGURE 1 | Map view of the Tomakomai carbon capture and storage (CCS)-CO 2 injection site and earthquake epicenters observed from 1990 through 2019 (from the Japan Meteorological Agency (JMA) earthquake catalog) overlaid on geological map (from (Japan, 1995)). The pink square indicates the injection site. Line A-A′ masks the profile in Figure 2.
Frontiers in Earth Science | www.frontiersin.org May 2022 | Volume 10 | Article 873645 dominated by sandstone, is a potential CO 2 reservoir. Another potential reservoir is the Takinoue formation at a depth of approx. 2,400 to 3,000 m, which is a volcanic formation composed of Neogene lava and tuff. Although the reservoir pressure is hydrostatic in the Moebetsu formation, the Takinoue formation is characterized by heterogeneously distributed overpressure zones. The temperature gradient is 3.7-4.4°C/100 m (Kano et al., 2014). The injection was started in April 2016 and ended in Oct. 2019, achieving a cumulative injection of 300,012 tons in the Moebetsu reservoir. The maximum injection rate was approx. 18,000 tons/month, and the maximum pressure increase was 12.6 MPa. Only 98 tons was injected into the Takinoue reservoir by two short tests (METI, NEDO, and JCCS, 2020). Thus, we only focus on the injection in the Moebetsu formation in this study.

Two-Dimensional Radial Flowing Model
As shown in Figure 2, the Moebetsu sandstone aquifer in the CCS site is quite stable and relatively flat over a few tens of kilometers. Thus, it is valuable to make a rough estimate of some key parameters using simple models for which theoretical solutions are available. To this end, a two-dimensional (2D) radial flowing model of an infinite and isotropic horizontal layer was tested. Under the assumption that Darcy's law holds, the theoretical solution of the pore pressure change is given as (Barker, 1988): where Γ represents the Gamma function, t (s) is the time since the start of injection, r (m) is the distance from the injection well, Q (m 3 /s) is the pumping/injection rate, H (m) and k (m 2 ) are the thickness and permeability of the layer, D (m 2 /s) is the hydraulic diffusivity, η (Pa·s) is the dynamic viscosity of water, ϕ is the porosity, S a (Pa −1 ) is the unconstrained specific storage coefficient, and β fl (Pa −1 ) and β pv (Pa −1 ) are the compressibility of the fluid and pores, respectively. This kind of model has been successfully applied to represent fluid pressure and strain observed during pumping tests in soft sedimentary formations (Lei et al., 2019). The thickness of the Moebetsu formation surrounding the injection well is greater than 100 m, and the permeability and porosity are 50 mD and 0.15, respectively (Kano et al., 2014). We assume that H = 100 m, ϕ = 0.15, and k from 1 to 100 mD. In order to perform an overestimation, the injection rate of water is assumed to be Q = 0.0143 m 3 /s, approx. 270,000 tons/year, which is greater than the actual mean injection rate of CO 2 (less than 100,000 tons/year). For the case of K = 10 mD, the predicated pressure buildup surrounding the injection well falls in the range from 22 to 25 MPa after one to 10 years of injection, which is greater than the actual maximum injection pressure (12.6 MPa). Even so, the fluid pressure front (pressure increased by 0.01 MPa) is far from the ITTFZ, even at the end point of the 10-years injection ( Figure 3A). As shown in Figure 3B, the pressure front increases with increasing K. Even for the worst case of K = 100 mD, the front did not reach the fault at 780 days after injection started.

Three-Dimensional Thermal-Hydraulic-Mechanical-Coupled Simulation
The coupled THM simulation was used to predict injectioninduced changes in rock properties, formation deformation, stress redistribution, and fracture/fault stability. In the present study, Tough2 and Flac3D were selected for the coupled THM simulation. Tough2 is a multiphase reservoir simulation program developed by the U.S. Lawrence Berkeley National Laboratory (LBNL). Flac3D is commercially available software for stress analysis (Itasca, 2000). As a promising combination, the "Tough-Flac3D″ approach with couplers and post-processing tools has proven useful in the analysis of deformation accompanied with fluid flow within reservoirs of hard and soft rocks (Itasca, 2000;Todesco et al., 2004;Rutqvist et al., 2008;Rutqvist et al., 2015;Sorai et al., 2015;Lei et al., 2017).
In the coupled simulation, permeability is revised within every time step by built-in functions. Following previous works, permeability has been expressed as a function of volumetric strain (ε v ) (Chin et al., 2000;Cappa and Rutqvist, 2011) as.
Where ϕ and k are porosity and permeability, respectively, with ϕ i and k i being the initial values. A n of~20 results in a 2-order-ofmagnitude permeability increase for a fully reactivated fault, as estimated by in situ testing (Ohtake, 1974) and laboratory experiments (Alam et al., 2014). A numerical model covering an area of 80 km × 80 km centered at the Tomakomai CCS injection site was constructed. The ITTFZ is taken as the thin fault zone. Since the upper limitation of temperature of Tough2 is 300°C, we only focus on the uppermost 6 km ( Figure 4). The model is divided into grids by steps so that the volumes of the fault and injection well are split into small parts while the surrounding matrix becomes coarser as the distance from the fault increases. The total number of grids and elements are 18,954 and 16,796, respectively. The hydraulic and mechanic properties for all nine formations and the fault zone are listed in Table 1. The hydraulic properties are similar to that used in previous studies (e.g., Kano et al., 2014). The upscaled mechanical properties are based on laboratory tests on core samples from the site at 10 MPa confining pressure. The most sensitive parameters are the permeability of porosity of the Moebetsu formation, which were firstly determined based on the permeability profile estimated from injection tests using salt water and logging data (METI, NEDO, and JCCS, 2020) and then optimized by tuning simulation. Aided by the 2D radial flow analysis results, the determination of the optimal parameters required only a few tuning steps. At the same time, we also performed calculations for different permeability (10, 20, 50, 100 mD) and porosity (0.1, 0.15) to examine the influence of the uncertainty of these sensitive parameters on the results. Due to the lack of monitoring data on rock mechanical responses, in-depth studies on the uncertainty of mechanical parameters cannot be carried out. The only water injection pressure data can be well fitted by modifying the reservoir hydraulic parameters (METI, NEDO, and JCCS, 2020).
First, simulation was carried out without any injection in order to obtain a hydraulic and thermal steady state. The top layer (air) has a constant pressure (0.1 MPa) and temperature (10°C), and the temperature gradient is assumed to be 4°C/100 m. Roller boundary conditions were imposed on the four sides and bottom of the model. The fluid pressure at the four sides was assumed to be hydrostatic and was kept constant during injection simulation. Both the flow rate and heat flux were fixed at the lower boundary. The in-site stress values are installed in all zones, and also applied as loads acting on the far-field boundaries. The vertical stress σ V  was calculated from assumed density values. Because the main rock of the model is soft rock, the maximum and minimum horizontal stress values were assumed as σ H = 1.3 σ V , and σ h = 1.2 σ V , respectively. Since the main purpose of the present study is to assess the impact of water injection activities on faults, the background tectonic stress itself is not important. During the coupled simulation, the permeability increases with increasing volumetric strain. The peak well head pressure during injection is 12.6 MPa (METI, NEDO, and JCCS, 2020) and decreases rapidly with increasing distance from the injection zone. Thus, we only consider the elastic behaviors.
In addition to the injection of CO 2 in the Moebetsu formation, a water source was added at the bottom of the fault to simulate natural flow along the ITTFZ. First, fluid flow along the fault zone was carried out at an injection rate of 10,000 tons/year for 1,000 years to reach a relatively stable state of fluid flow and heat flux. The CO 2 injection was then started at an injection rate of~183,600 tons/year for 100 years. This is the base model (Model-1). For comparison, a simulation was also carried out without CO 2 (Model-2) to verify whether CO 2 injection alters deep fluid flow along the fault.
The simulation results show that after three years of injection, by which point more than 550,000 tons of supercritical CO 2 were injected into to the Moebetsu formation, the injected CO 2 is distributed within a few hundred meters around the injection well ( Figure 4D), which is slightly larger than that revealed by the 3D seismic exploration implemented in 2019 (METI, NEDO, and JCCS, 2020). This result is expected since larger injection rate was assumed in the simulation. The pressure buildup is still beyond the fault for both models (Figures 4B,C). As compared with Model-1, the pressure front of Model-2 is farther from the fault zone due to the background flow from deep in the fault. The change of the Coulomb Failure Stress (CFS) acting on the fault planes, which have the same strike, dip, and rake with the seismogenic fault of the 2018 Hokkaido Eastern Iburi earthquake, is governed by fluid pressure. Although there is a redistribution of flow velocity at the intersection of the fault and the reservoir, the resulted fluid pressure (P(A) in Figure 5) by Model-1 of different permeability and porosity for the Moebetsu formation is exactly the same as Model-1, consistently show that CO 2 injection has no effect on fluid pressure at the deep part of the fault event after 100 consecutive years of injection ( Figure 5).

DISCUSSION AND CONCLUSION
The results of the present study show that, even under extreme conditions, when CO 2 is injected at a rate of~180,000 tons/ year (actually less than 100,000 tons/year) for 10 consecutive years (actually less than three years), the 0.01-MPa front of the fluid pressure increment is far from the hypocenter of the Iburi Earthquake. After three years injection, the injected CO 2 is distributed within a few hundred meters around the injection well, which is slightly larger than that revealed by the 3D seismic exploration implemented in 2019 (METI, NEDO, and JCCS, 2020) after injection 207,208.9t CO 2 . If the injected CO 2 , including dissolved CO 2 , could rapidly migrate beyond 10 km as claimed in Sano et al. (2020), 3D seismic exploration would either see nothing or show anomalies on a large area, but this is not the case. Detailed discussion on geochemical issues is interesting but beyond the scope of this study. In the absence of wide-area groundwater flow, the distribution of gaseous and dissolved CO 2 is also limited within a small zone of only a few kilometers at most. Even if there is a wide area of groundwater flow, the distribution range of CO 2 has expanded, compared with The results of different permeability and porosity for the Moebetsu formation consistently show that CO 2 injection has no effect on fluid pressure at the deep part of the fault event after 100 consecutive years of injection.
the stable fluid pressure before water injection, the range where the fluid pressure has increased still does not change significantly. The changes in fluid pressure and pore-elastic effects at the intersection of the ITTFZ and the Moebetsu layer 20 km away are much smaller than the changes caused by daily tidal action, not to mention the 37-km depth of the fault, where the 2018 Hokkaido Eastern Iburi Earthquake initiated. For reference, tidal stress amplitude thresholds required to trigger earthquakes, for static and dynamic triggering, have been estimated to range from 0.01 to 0.03 MPa (pore pressure: 0.015-0.05 MPa) (King et al., 1994;Cochran et al., 2004). Therefore, the implementation of the Tomakomai CO 2 geological storage demonstration project is unlikely to be the triggering factor of the earthquake.
In our simulation, the dissolution of injected CO 2 into the groundwater of the aquifer was counted but other chemical processes, such as mineralization were ignored. In general, minerals, such as silica in the deep source hydrothermal fluid, precipitated to the pressure drop caused by earthquake to heal the fault reduce the permeability of the fault (Saishu et al., 2017). However, no convincing mechanism exists that would allow fluid channels to be healed within a short period of time and block the natural fluid flow along the fault during CO 2 injection, which increases the fluid pressure.
As mentioned above, the geological sequestration of CO 2 is important in mitigating the global warming trend. Demonstration projects of various scales are for the development of technology and the verification of the safety and effectiveness of the corresponding technology. With regard to risk-related injection-induced fault reactivation and earthquakes, both oversized and undersized assessments must be avoided.
The model proposed by Sano et al., (2020) is worth investigating for commercial applications of large-scale and long-term CO 2 storage. Thermally, hydraulically, mechanically, and chemically coupled geomechanics modeling will be a key technology.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
The author confirms being the sole contributor of this work and has approved it for publication.