Global warming has been identified as one of the most serious global environmental issues for the last several decades. According to the Climate Change 2014 by IPCC  , the globally averaged combined land and ocean surface temperature anomaly has risen by around 0.85˚C from 1880 to 2012. Moreover, it is likely to rise by 2.6˚C - 4.8˚C by the end of the 21st century. As a result, the temperature rise causes the thaw of the mountain glaciers and snow cover globally, leading to the sea level change. Over the period from 1901 to 2010, the globally averaged sea level has risen by 0.19 m, and will further rise by around 0.07 - 0.36 m by 2050, and around 0.09 - 0.69 m by 2080  .
As the major cause for global warming, anthropogenic CO2 emission into the atmosphere has increased dramatically over the past few decades, and caused negative and irreversible effects on the environment and ecosystems  . In order to mitigate global warming, CO2 capture and storage (CCS) is widely regarded as one of the most effective countermeasures for CO2 emission reduction. The conventional methods for CCS usually refer to CO2 injection and storage into the deep saline aquifers onshore and shallow offshore      , which have been considered as the main stream for CCS. However, there is a social concern about the stability and security of CO2 storage in the deep saline aquifers  . If the cap rocks above the reservoir crack due to the earthquake, the injected CO2 may not be able to keep stable in the deep saline aquifers, and seep out of the cap rocks through the cracks. In the worst-case scenario, the injected CO2 may leak into the ocean, and cause ocean acidification    .
A novel approach of sub-seabed CO2 storage in the form of gas hydrate attracts much attention nowadays  . In this technology, the high-permeability sand layers surrounded by the low-permeability mud layers in the shallow sub-seabed sediments are selected as the target geological strata for CO2 storage, because the sand layers can serve as the reservoir, and the mud layers can serve as the overburden and underburden to restrain CO2 leakage. During the injection, CO2 flows into the pore space of the sand layers, and forms CO2 hydrate gradually with the underground water under low temperature and high pressure conditions. As a result, CO2 can be trapped stably inside the solid hydrate, and the risk of CO2 leakage can also be reduced greatly. However, there is also an obstacle for this new technology. After the injection, the permeability of the sand layers near the injection well will drop sharply due to the hydrate formation, which may cause the CO2 flow blockage, and hinder the further CO2 injection. Therefore, in order to store a large amount of CO2 in the sub-seabed sediments by gas hydrate, it is important and essential to choose a proper CO2 injection rate, and ensure CO2 to spread over a wide area after the injection.
Although the previous researchers have conducted a lot of studies on the topic of sub-seabed CO2 storage by gas hydrate, they mainly focused on the lab-scale experiments and simulations  -  . The behaviors of CO2 hydrate formation in the real sub-seabed sediments with two-phase flow still remain unknown, because few investigations have been dedicated to the large-scale simulations for a relatively long injection period measured by month.
For the practical application of this new technology, a numerical simulator incorporated with an integrated model for CO2 hydrate formation in the sand sediments was developed in our previous study  . Then, this simulator was employed to the experimental cases for the determination of the unknown model parameters. In this study, this simulator is applied to the numerical simulations of CO2 injection and CO2 hydrate formation in a large-scale geological model for an injection period as long as 150 days (nearly five months), to estimate CO2 storage capacity in the real sub-seabed sediments by gas hydrate.
2. Numerical Modeling
2.1. Governing Equations
The numerical simulator used in this study was developed by modifying a gas-liquid two-phase flow code, TOUGH + HYDRATE v1.0  . This improved simulator can describe the mass balance for water, gas, hydrate, heat, and CO2 mass fraction in the aqueous phase using the finite difference method. In this simulator, a series of five primary variables (P, T, SA, SG, and ) are solved iteratively by Newton-Raphson Method using five governing equations (mass balance equations for aqueous, gas, and hydrate phases, heat balance equation, and CO2 mass balance equation in the aqueous phase) as below, respectively  :
where is the volume fraction (i.e. saturation) of phase β º A, G, H (m3/m3), is the density of phase β º A, G, H (kg/m3), and is the mass fraction of the component κ º H2O, CO2, hydrate in phase β º A, G, H (kg/kg). is the flux term of phase β º A, G (kg/m3/s), is the diffusion term of phase β º A, G (kg/m3/s), is the source/sink term of phase β º A, G (kg/m3/s), is the specific enthalpy of phase β º A, G (J/kg), and is the specific internal energy of phase β º A, G, H (J/kg). P is the pressure (Pa), T is the absolute temperature (K), and is the composite thermal conductivity (W/m/K). f is the porosity of the porous medium (−), is the density of the porous medium (kg/m3), and is the specific heat capacity of the porous medium (J/kg/K). is the total hydrate formation rate (kg/m3/s), and is the enthalpy change during hydrate formation/dissociation (J/kg). is CO2 dissolution rate in the aqueous phase (kg/m3/s), is the enthalpy change during CO2 dissolution (J/kg), and is CO2 injection rate (kg/m3/s).
2.2. CO2 Hydrate Formation Model
As mentioned before, in our previous study  , an integrated model for CO2 hydrate formation in the sand sediments was proposed to predict hydrate formation morphologies based on the formation locations in the sand sediments. In this model, CO2 hydrate is assumed to form at three different locations in the sand sediments, and the total hydrate formation rate is described as below:
where δ is a switch to determine whether the gas front exists in a computational cell (δ = 1) or not (δ = 0). , , and are the corresponding hydrate formation rates on the gas front, on the hydrate film, and on the surface of the sand particles behind the gas front (kg/m3/s), respectively, which are given as below:
where is the molar mass of CO2 hydrate (kg/mol), and is the intrinsic rate constant of CO2 hydrate formation (mol/m2/Pa/s). and are the rupture ratios on the gas front and behind the gas front (−), respectively. , , and are the gas-liquid interfacial area on the gas front and behind the gas front, and the sand surface area (m2/m3), respectively. , , , , and are CO2 fugacity in the gas phase, in the aqueous phase, at the three-phase equilibrium point, at the gas-liquid interface on the gas front and behind the gas front (Pa), respectively. is the hydrate formation rate transferred from to (kg/m3/s). For each sub-model in this integrated model, one can find all the details in our previous study  .
3. Model Construction
3.1. Geological Model
For the large-scale geological model simulating the real sub-seabed sediments, an axisymmetric cylinder with a radius of 100 m and a thickness of 160 m is built in this study, as shown in Figure 1. This sediment model is assumed to be located at the depth of 870 - 1030 m from the sea surface (at the water depth of 500 m), and divided into three domains from top to bottom: 1) overburden, 2) CO2 storage reservoir, and 3) underburden. The CO2 storage reservoir is set to be composed of sand layers with the thickness of 100 m, referring to the shallow reservoir (approximately 100 m thick) used for Tomakomai CCS Demonstration Project of Japan  . The overburden and underburden are both set to be composed of mud layers with the thickness of 30 m. According to the previous report
Figure 1. Schematic diagram of the sub-seabed sediment model built in this study.
by Sun et al.  , it may be sufficient for these 30-m-thick overburden and underburden layers to simulate the boundary effects of heat exchange and pressure propagation.
In addition, an injection well is located at the center of the sediment model, which is used for CO2 injection. In order to determine the length of the injection part, the CO2 flow direction in the reservoir has been considered as a main factor in this study. After the injection, CO2 will not only flow horizontally in the reservoir, but also flow upward due to the buoyancy. If the length of the injection part is set to be smaller than the thickness of the CO2 storage reservoir, the injected CO2 may only distribute and form hydrate at the upper part of the reservoir, or in the vicinity of the injection well, which will increase the risk of the CO2 flow blockage. Based on this reason, the injection part is also set to be 100 m, which equals to the thickness of the CO2 storage reservoir, to ensure that CO2 can spread over a wide area after the injection, and flow smoothly in the reservoir without the CO2 flow blockage.
3.2. Computational Conditions
The pore water pressure of the sediment model is assumed to be hydrostatic, and the initial hydrostatic pore water pressure (MPa) can be calculated according to the empirical equation as below  :
where is the standard atmospheric pressure (MPa), is the sea water density (kg/m3), g is the gravitational acceleration (m/s2), h and z are the water depth, and the depth of the sediments from the seafloor (m), respectively.
For the initial temperature condition in the sediment model, the geothermal gradient is taken into account. Other main physical properties of the sediment model refer to the field parameters used for the numerical simulations of gas production behavior from methane hydrate reservoir at the first offshore test site in the eastern Nankai Trough, Japan (2013)  . The model parameters used in this study are listed in Table 1. It is worth mentioning that by the preliminary calculations, it is confirmed that the initial pressure and temperature conditions in the CO2 storage reservoir are located within the hydrate stability zone according to the phase diagram of CO2 hydrate proposed by Kamath  . This indicates that the initial pressure and temperature conditions determined in this study are appropriate and suitable for CO2 hydrate formation, and hydrate can form in the reservoir after the injection. Besides, the intrinsic permeability of the overburden and underburden is set to be 100 times smaller than that of the CO2 storage reservoir, so that the overburden and underburden can both serve as the low-permeability layers to restrain CO2 leakage.
4. Simulation Results and Discussion
In order to estimate CO2 storage capacity in the sediment model built in this study, a proper CO2 injection rate needs to be determined in advance. By the preliminary simulations, it is found that if the injection rate is set to be larger than 100 ton/day, the CO2 flow blockage will occur at the early stage of the injection process. Therefore, three moderate injection rates of 10 ton/day, 50 ton/day, and 100 ton/day are chosen for Case 1 - Case 3. These three cases are also used for the sensitivity analysis to investigate the influence of the injection rate on the behaviors of CO2 reaction and hydrate formation. Besides, the injection period is set to be 150 days (nearly five months), so the total amounts of CO2 injection are 1500 ton, 7500 ton, and 15,000 ton, respectively.
Table 1. The physical properties of the sediment model.
4.1. Evolutions of CO2 Reaction, Free CO2, and Hydrate Formation in the Sediment Model
Figure 2 shows the evolutions of CO2 reaction, free CO2, and hydrate formation in the sediment model during the whole injection period for Case 1 - Case 3. As can be observed in Figure 2(a), the mass rates of CO2 reaction reach the peaks at the beginning of the injection process for all the three cases. The reason is that, as mentioned before, the initial pressure and temperature conditions in the CO2 storage reservoir are suitable for CO2 hydrate formation, and as soon as CO2 is injected into the reservoir, it forms hydrate immediately, leading to the jump in the mass rate of CO2 reaction. On the other hand, during the processes of CO2 dissociation into the aqueous phase and CO2 hydrate formation, a large amount of heat is released in a short time, resulting in the abrupt temperature rise in the sediments, which will have a negative effect on the further hydrate formation. Therefore, the mass rate of CO2 reaction drops sharply after the peaks, and the heat release decreases accordingly. This phenomenon repeats again and again, and it is the reason that there are many fluctuations on the curves of the mass rates of CO2 reaction for all the three cases, especially in Case 3. In addition, it can also be seen that the larger the injection rate is, the more intense the fluctuations become, because the amount of heat release also becomes larger. During the whole injection period, the average mass rates of CO2 reaction in the reservoir are approximately 0.20 ton/day, 0.90 ton/day, and 1.86 ton/day, respectively, for Case 1 - Case 3.
By the integral of the mass rate of CO2 reaction, the amount of CO2 reaction over time can be obtained accordingly, as shown in Figure 2(b). Since the curves of the mass rates of CO2 reaction are fluctuant in Figure 2(a), the curves of the amount of CO2 reaction are also fluctuant, especially in Case 3, and present linear behaviors by appearance for 150 days. If given a much longer injection period, the status in the sediments will get close to the CO2 flow blockage gradually. As a result, the amount of CO2 reaction may approach to a maximum value, and the curves may present exponential behaviors by appearance. By the end of the whole injection period of 150 days, the total amounts of CO2 reaction in the reservoir reach approximately 32.9 ton, 138.0 ton, and 286.5 ton, respectively, for Case 1 - Case 3.
After the injection, a part of CO2 neither forms hydrate nor dissolves into the aqueous phase. Instead, it just remains in the reservoir as free CO2. Figure 2(c) shows the amount of free CO2 in the reservoir over time. As can be seen in the figure, there are also some fluctuations on the curves for all the three cases. However, the fluctuations in Figure 2(c) are just opposite to those in Figure 2(b), because the more CO2 participates in the reaction and forms hydrate, the less it will remain in the reservoir as free CO2. With the increase of the injection rate, the amount of free CO2 also rises accordingly. By the end of the whole injection period, the total amount of free CO2 in the reservoir reach approximately 104.2 ton, 646.3 ton, and 1310.2 ton, respectively, for Case 1 - Case 3.
Figure 2. Evolutions of CO2 reaction, free CO2, and hydrate formation in the sediment model during 150 days for Case 1 - Case 3. (a) The mass rate of CO2 reaction [ton/day]; (b) The amount of CO2 reaction [ton]; (c) The amount of free CO2 [ton]; (d) The amount of CO2 hydrate formation [ton].
For the amount of CO2 hydrate formation in the reservoir, the curves in Figure 2(d) are actually in the same shapes as those in Figure 2(b), because the amount of CO2 reaction means the part of CO2 which forms hydrate, and it is closely related to the amount of CO2 hydrate formation. Besides, it is clearly observed that the more CO2 is injected into the reservoir, the more CO2 hydrate forms. By the end of the whole injection period, the total amount of CO2 hydrate formation in the reservoir reach approximately 110.5 ton, 462.7 ton, and 960.7 ton, respectively, for Case 1 - Case 3.
4.2. Spatial Distributions of Physical Properties in the Sediment Model
Since the total amounts of CO2 reaction and CO2 hydrate formation are the largest in Case 3, this case has been extracted as a best case to investigate the behaviors of CO2 hydrate formation in the sediments with two-phase flow. Figure 3 shows the spatial distributions of the physical properties in the sediment model after 30 days, 90 days, and 150 days, respectively, for Case 3.
As can be seen in Figure 3(a), after the injection, the isopiestic lines in the sediment model have been disturbed, and showed an upward shift in the vicinity
Figure 3. Spatial distributions of the physical properties in the sediment model during 150 days for Case 3. (a) Pressure [MPa]; (b) Temperature [˚C]; (c) CO2 mass fraction in the aqueous phase [kg/kg]; (d) CO2 saturation [m3/m3]; (e) Hydrate saturation [m3/m3]; (f) Water saturation [m3/m3].
of the injection well, especially at the depth of 895 - 910 m. This is because that after CO2 hydrate forms in the reservoir, the solid hydrate occupies the pore space of the sediments, and causes the reduction in the effective permeability, which will hinder the CO2 flow to a certain extent. As a result, the injected CO2 accumulates in the vicinity of the injection well, leading to the upward shift of the isopiestic lines as mentioned before. On the other hand, during the injection, the temperature jumps significantly in the reservoir, and a high temperature zone forms in the sediments, as shown in Figure 3(b). For the heat source of the high temperature zone, some of the heat comes from CO2 hydrate formation heat, but most of the heat actually comes from CO2 dissociation heat into the aqueous phase. This can be directly proved by Figure 3(c), in which it is clearly observed that during the process of the two-phase flow, a narrow zone of high CO2 mass fraction due to CO2 dissociation into the aqueous phase appears in the reservoir. This narrow zone moves forward with the CO2 flow, and expands gradually in the radial direction, which indicates that the amount of CO2 dissociation becomes larger during the process of the two-phase flow, and releases a great deal of heat in the reservoir. Meanwhile, due to the low thermal conductivities of the sand layers and mud layers, the heat generated by CO2 hydrate formation and CO2 dissociation increases and accumulates in the reservoir instead of being discharged from the sediments, resulting in the high temperature zone as shown in Figure 3(b).
In order to obtain a better understanding of the behaviors of CO2 hydrate formation in the sediments with two-phase flow, the evolutions of CO2 saturation, hydrate saturation, and water saturation over time are presented in Figures 3(d)-(f), respectively. As can been seen in Figure 3(d), after the injection, CO2 nearly moves horizontally in the reservoir, leading to the graded distribution of the CO2 saturation in the radial direction. However, since the density of CO2 is smaller than that of water, the injected CO2 also flows upward gradually. Although the intrinsic permeability of the overburden and underburden is set to be 100 times smaller than that of the CO2 storage reservoir, a small part of CO2 still seeps out of the reservoir, and leaks into the overburden due to the buoyancy. On the other hand, as shown in Figure 3(e), CO2 hydrate forms gradually in the reservoir as a result of the suitable temperature and pressure conditions, and generates a lot of heat as mentioned before. This part of CO2 hydrate formation heat, especially along with the large amount of CO2 dissociation heat, has a negative effect on the further hydrate formation, so the hydrate saturation in the reservoir shows a trend of heterogeneous distribution. Meanwhile, since the temperature and pressure conditions in the overburden are also within the hydrate stability zone, most of the leaked CO2 forms a thin layer of hydrate in the overburden, especially just above the injection well, where a high hydrate saturation spot can be observed clearly. This thin layer of hydrate may serve as a self-sealing cap, and restrain the further CO2 leakage, so that the injected CO2 can be stored safely in the sub-seabed sediments without leaking into the ocean. Besides, after the injection, CO2 pushes the aqueous phase forward in the radial direction, and a part of pore water is consumed to form hydrate, so the water saturation also displays a graded distribution, creating a low water saturation zone in the vicinity of the injection well, as shown in Figure 3(f).
For the estimation of CO2 storage capacity in the real sub-seabed sediments by gas hydrate, a large-scale geological model simulating the real sub-seabed sediments in the ocean was built in this study, and numerical simulations of CO2 injection and CO2 hydrate formation in the sediments with two-phase flow were conducted at three different injection rates of 10 ton/day, 50 ton/day, and 100 ton/day, respectively. It is found that, at the injection rate of 100 ton/day, a total amount of 15,000-ton CO2 can be injected into the sediments for an injection period of 150 days. After the injection, a part of CO2 can be stored in the sediments in the form of gas hydrate, and the rest part remains in the reservoir as free CO2 or dissolves into the aqueous phase. For a CO2 storage reservoir with the thickness of 100 m as built in this study, at the injection rate of 100 ton/day, i.e., averagely 1 ton/day/m, a maximum amount of 36,500-ton CO2 could be injected and stored in the sub-seabed sediments per year. For the practical scenario, this average value of 1 ton/day/m could also be used to determine the actual injection rate based on the thickness of the real sub-seabed sediments.
Moreover, in order to investigate the behaviors of CO2 hydrate formation in the sediments with two-phase flow, the spatial distributions of the physical properties in the sediments over time were presented for the case of the injection rate of 100 ton/day. The simulation results indicate that during the injection process, a large amount of heat is released due to CO2 hydrate formation heat and CO2 dissociation heat into the aqueous phase, leading to a high temperature zone in the reservoir which has a negative effect on the hydrate formation. After the injection, CO2 not only flows horizontally in the reservoir, but also flows upward due to the buoyancy. As a result, a small part of CO2 permeates into the overburden, forms hydrate, and serves as a self-sealing cap to restrain the further CO2 leakage. Although the long-term injection and monitoring are still needed to fully evaluate the potential and feasibility of the technology of sub-seabed CO2 storage in the form of gas hydrate, it is reasonable to believe that this novel technology can be expected to be applied in the field demonstration in the future.
This work was supported by the Hirosaki University Grant for Exploratory Research by Young Scientists and Newly-appointed Scientists, and the Open Fund of Key Laboratory of Ocean Energy Utilization and Energy Conservation of Ministry of Education at Dalian University of Technology.
*A part of this paper was presented at 14th International Conference Fluid Dynamics (ICFD14) in Sendai, 1-3 Nov. 2017.