Deep saline aquifers have strong heterogeneity under natural conditions, which affects the migration of carbon dioxide (CO2) injection into the reservoir. How to characterize the heterogeneity of rock mass is of great significance to research the CO2 migration law during CO2 storage. A method is proposed to construct different heterogeneous models from the point of view of whether the amount of data is sufficient or not, the wholly heterogeneous model with sufficient data, the deterministic multifacies heterogeneous model which is simplified by lithofacies classification, and the random multifacies heterogeneous model which is derived from known formation based on transfer probability theory are established, respectively. Numerical simulation is carried out to study the migration law of CO2 injected into the above three heterogeneous models. The results show that the migration of CO2 in heterogeneous deep saline aquifers shows a significant fingering flow phenomenon and reflect the physical process in CO2 storage; the migration law of CO2 in the deterministic multifacies heterogeneous model is similar to that in the wholly heterogeneous model and indicates that the numerical simulation of simplifying the wholly heterogeneous structure to the lithofacies classification structure is suitable for simulating the CO2 storage process. The random multifacies heterogeneous model based on the transfer probability theory accords with the development law of sedimentary formation and can be used to evaluate the CO2 migration law in unknown heterogeneous formations. On the other hand, by comparing the dry-out effect of CO2 in different heterogeneous models, it is pointed out that the multifacies characterization method will weaken the influence due to the local homogenization of the model in small-scale research; it is necessary to refine the grid and subdivide the lithofacies of the local key area elements to eliminate the research error. The research results provide feasible references and suggestions for the heterogeneous modeling of the missing data area and the simplification of large-scale heterogeneous models.
Carbon dioxide (CO2) sequestering in deep saline aquifers is considered to be an effective way to mitigate climate change, with great storage potential [1–5]. Before the injection of CO2 into deep saline aquifers, it is quite important to screen potential reservoirs and validate the injection strategy through large-scale numerical simulation. Supercritical CO2 (scCO2) injection into deep saline aquifers is a multiphase (gaseous, liquid, solid) and multicomponent (CO2, water, salt) seepage problem. Commonly used numerical software to simulate this process includes STOMP CO2 , ECLIPSE 300 , and TOUGH . Normally, heat exchange, phase change, and transport of CO2, water, and salt mixtures  must be considered in the numerical model. In addition, it is important to ensure the accuracy of the simulation and to ensure that the parameters of the model match with the real site.
To improve the computational efficiency of large-scale model simulation and the scarcity of real site parameters, the homogeneous model [9–12] has been used by many scholars to simulate the multiphase fluid migration process of CO2 sequestration in deep saline aquifers. Nordbotten et al.  established an analytical model related to fluid mobility and gas pressure distribution and solved the time-space evolution process of CO2. Oruganti and Mishra  used the analysis model based on fractal theory and numerical analysis method to predict the change of CO2 saturation (SCO2). Lagneau et al.  used data from CO2 sequestration experiments in the Paris Basin and the North Sea to simulate the transport, dissolution, and mineralization of CO2 in gas-liquid two-phase and phase state. The simplified layered model is used to describe the vertical heterogeneity of the formation, but it is still a simplified homogeneous model in a certain layer [15–17].
With the development of more and more site demonstration projects and field-scale projects on CO2 storage [18, 19], subsurface monitoring data from the site are widely collected (such as seismic methods [20, 21], electromagnetic and electric methods , and geophysical logs ), providing a reference to verify the accuracy of numerical simulation. The monitoring information not only confirmed some of the estimates from the numerical simulations but also found some differences from the numerical predictions. For example, cross-well seismic at Frio CO2 storage site finds that the front of the CO2 plume moved more quickly than simulation results conducted with a homogeneous model , reservoir modeling obtains unsatisfactory CO2 migration results, and saturation profile matches .
Due to the discontinuity of the actual sedimentary process of the rock mass, deep saline aquifers have a strong objective heterogeneity [26–28]. How to characterize the heterogeneity of rock mass is of great significance to the study of the transport law of CO2 in deep saline aquifers. At present, the methods of characterizing formation heterogeneity can be counted in two categories. One is to assign the hydrogeological parameters to each unit of the simulated area directly from the field data. For example, Doughty  simulated the distribution of CO2 by distributing the porosity and permeability acquired from the field to the simulated area and found that the distribution of CO2 preferred the area with higher permeability. Another method to characterize the heterogeneity of formation is to generate the distribution field of global hydrogeological parameters based on field data and stochastic theory. Ide et al.  simulated CO2 migration in heterogeneous formations with different correlation lengths. Flett et al.  used the standard geophysics method to generate a random distribution of porosity and to simulate the effect of heterogeneity on CO2 transport. A series of three-dimensional (3D) permeability random fields are generated based on a large numerical model by Obi and Blunt  to simulate the transport and capture of CO2, and it is found that strong heterogeneity can reduce the storage efficiency of underground porous media. Liu et al.  used the randomly distributed porosity of the log-normal distribution to characterize the heterogeneity of the deep saline aquifers and simulated the pressure and distribution characteristics of CO2. The above methods only consider how to construct a heterogeneous model and study the effect of heterogeneity on the regularity of the simulation results. These models are not supposed to represent real site heterogeneous features; how to build a model consistent with the in situ heterogeneity as far as possible from the data available is a problem that needs to be further solved at present.
In summary, how to characterize the heterogeneous structure of deep saline aquifers is the key to studying the CO2 migration law. Because the buried depth of the potentially suitable CO2 storage formation site is more than 800 m, it is difficult to obtain detailed hydrogeological parameters. This paper proposes that according to the detailed degree of field data available, the characterization of heterogeneous deep saline aquifers can be simplified in varying degrees, as shown in Figure 1. Specific ideas are as follows: (I) when the measured data are more detailed, a wholly heterogeneous model can be established, and different hydrogeological parameters can be assigned to each element of the grid in the model (Figure 1(a)); (II) when there are few measured data, the lithology of the deep saline aquifers can be roughly classified according to the measured hydrogeological parameters, and the rock masses with the same hydrogeological characteristics can be classified into the same kind of lithofacies, the structure of multifacies combination is called multifacies heterogeneous model, which simplifies the simulation modeling (Figure 1(b)); (III) when there are little or no data in the target area, according to the adjacent known formation data, the spatial distribution of multifacies permeability can be randomly generated based on transfer probability theory to characterize the heterogeneity of research target formation, establish a random multifacies heterogeneity model (Figure 1(c)). In addition, a large-scale numerical simulation case analysis is used to verify the feasibility and applicability of the method based on simplified multifacies heterogeneity and random multifacies characterization based on adjacent area data.
2.1. Establishment of Heterogeneous Model Based on Sufficient Field Data
In general, the mesh size of numerical simulation should be set as which has no mesh dependence on the problem solved. Under the condition of satisfying this premise, if each grid has corresponding hydrogeological parameters, then it is said that the data is sufficient. Based on a set of comprehensive geostatistics data of permeability provided by GSLIB , these data are assumed to be the permeability of a real reservoir. A two-dimensional (2D) heterogeneous model of is established, in which each central node of the grid is assigned different permeability parameter values; the distribution of permeability values is shown in Figure 2. We named a wholly heterogeneous model in this research.
2.2. Establishment of a Heterogeneous Model Based on Lithofacies Classification
Based on the measured hydrogeological parameters (such as permeability and porosity ), these hydrogeological parameters are classified in a certain way to represent different stratigraphic and lithological characteristics, which are called different types of lithofacies, and the structure of the multifacies association is called a multifacies heterogeneous model.
Based on the geostatistics data of the same permeability as in Section 2.1, the data are classified into four groups according to the numerical value of the permeability. The permeability of the first 30% is called lithofacies 2, the permeability of the second 40% is called lithofacies 1, the permeability of 20% is called lithofacies 3, and the minimum permeability of 10% is called lithofacies 4. A 2D model of is established. The distribution of lithofacies is shown in Figure 3. It is called the deterministic heterogeneous model of multifacies. We named the multifacies heterogeneous model in this research.
2.3. Establishment of a Heterogeneous Model Based on Transition Probability Theory
However, the usual situation is that we usually only get partial field data in research. When there are no data or the data is not enough, how to get partial field data can roughly reflect the sedimentary law of sedimentary rocks? The problem to be solved in this paper is to extend permeability distribution of the known formation to the unknown formation and to obtain close numerical simulation results for large-scale simulation.
Based on the transfer probability of correlation points, the transfer probability of different lithofacies under certain correlation length is calculated by the Markov chain; then, the conditional probability of a lithofacies in the center of the model grid is simulated by the sequential indicator simulation (SIS) method, and the optimal indicator lithofacies in the center of the model grid is determined by the objective function.
2.3.1. Transition Probability Theory
2.3.2. Markov Chain
Markov chain is a mathematical model used to describe a certain random phenomenon and to study the natural science process based on mathematical analysis, which assumes that what will happen in the future is related to what happens now, but it has nothing to do with what happened. Taking a simple one-dimensional (1D) Markov chain as an example, Figure 5 shows the distribution characteristics of three types of lithofacies in a 1D embedded system.
2.3.3. Sequential Indicator Simulation (SIS) Method
2.3.4. Characterization Process of Random Multifacies Heterogeneous Structure
For each type of lithofacies that may occur in each grid center node, the objective function is calculated, and when the objective function is minimum, the corresponding lithofacies is the optimal indicator lithofacies of the grid center node; thus, the indicator lithofacies of each grid center node is obtained. When the randomly generated random paths through all the central nodes of the grid are different, the generated multifacies heterogeneous model is different.
Specifically, the 2D deterministic multifacies heterogeneous model is extended to a 3D random multifacies heterogeneous model. First, the transition probability is calculated according to the existing data in Section 2.1, and a continuous Markov model is generated, as shown in Figure 6.
In Figure 6, according to the existing permeability data, the transfer probability of a certain type of lithofacies under the premise of the occurrence of these lithofacies or other lithofacies is calculated under the influence of different correlation lengths. For example, the sequence of the first 4 diagrams shows the transfer probability sequence of lithofacies 1 under the influence of different correlation lengths of 0 ~ 20 m and the premise of the occurrence of lithofacies 1, lithofacies 2, lithofacies 3, or lithofacies 4. The 4 diagrams in the first column show the transfer probability sequence of lithofacies 1, lithofacies 2, lithofacies 3, or lithofacies 4 under the influence of 20 different correlation lengths. The continuous graph in Figure 6 shows that based on the calculated discrete transition probability, the continuous transition probability value is generated by the Markov model; thus, the transition probability of any length within 20 m can be obtained. According to the transfer probability matrix in different directions, the transfer probability in 3D space can be obtained. The SIS method is used to generate a random path through all grid nodes. Combined with the determination of the objective function, the optimal indicative lithofacies distribution of all grids under the random path is generated; that is, a spatially random multifacies heterogeneous model is obtained, as shown in Figure 7.
The 2D random distribution characteristics of permeability classification in one of the cross sections (different from the deterministic multifacies heterogeneous model section) are assigned to the 2D model of , which is the random multifacies heterogeneous model in this study, as shown in Figure 8.
The randomly generated multifacies heterogeneous model considers the spatial cross-correlation of geological bodies and accords with the development law and evolution process of sedimentary formation.
3. Simulation of CO2 Migration in Different Heterogeneous Models
The above three heterogeneous models are built based on different data volumes. In general, the more data quantity, the constructed model will match the real site situation better. In the case of little or no data, the multifacies model and the random model are established, respectively. Several questions need to be studied: (1) why use the heterogeneous model to study CO2 sequestration? (2) Can flow characteristics be simulated by simplified multifacies models and random models? (3) How does the simplified model differ from the wholly heterogeneous model? (4) What should be paid attention to when using this different heterogeneous model?
A large-scale simulation of CO2 migration in different heterogeneous models is conducted in this section. The geometric size of the numerical model is the same as that of the heterogeneous model; the upper and lower boundary of the model is set as the no-flow boundary. Considering the injection well buried depth is 1,200 m underground, the formation pressure is 12 MPa, and the formation temperature is 40°C, CO2 is injected into a horizontal well at a distance of 22 m from the lower boundary, and the injection rate is 0.25 kg/s; at this point, CO2 is in a supercritical state. The rest of the left boundary is the no-flow boundary; the right boundary is the constant pressure boundary. The horizontal grid of the model is divided into 58 elements, the dimensions of the first 8 elements are 1 m, 1 m, 2 m, 2 m, 4 m, 5 m, 10 m, and 15 m, respectively, and the dimensions of the next 49 elements are 40 m. The last element size is 0.001 m (as a constant pressure boundary), and the buried depth grid is divided into 50 elements with a unit size of 4 m.
The migration law of CO2 with three heterogeneous models is simulated by using the ECO2N module in the TOUGH2 program. The total simulation time is 2 years, and the plume of gas saturation () distribution of CO2 and the isoline map of gas pressure is drawn. The parameters of the deterministic multifacies heterogeneous model, random multifacies heterogeneous model, and wholly heterogeneous model are shown in Table 1.
In Table 1, represents rock permeability, and φ represents rock porosity. The relative permeability equation is based on the van Genuchten-Mualem model [36, 37]; capillary pressure equation is based on the van Genuchten model . In relative permeability equation, is the empirical parameter related to core size distribution, is the residual liquid saturation, is the residual gas saturation, and is the liquid saturation when saturated. In capillary pressure equation, is the air entry pressure, and is the maximum capillary pressure. For the initial condition, the formation pressure () and temperature () of the injection well are constant, is initial gas saturation, is brine salinity, and is the CO2 injection rate.
4. Simulation Results and Discussion
4.1. Simulation Results in Different Heterogeneous Models
Figure 9 shows the distribution and the isoline map of gas pressure after 2-year CO2 injection into horizontal wells in the wholly heterogeneous model.
Due to the low water content in the gas, the distribution shown in Figure 9(a) can be considered as a CO2 plume. When the formation pressure is 12 MPa and the temperature is 40°C, the density of scCO2 is about 715 kg/s, which is lower than that of brine. After CO2 is injected into the wholly heterogeneous aquifer, the gravity difference causes CO2 to suffer upward buoyancy, so the distribution characteristics of CO2 show a funnel shape as a whole. Because of the heterogeneity of the formation, when CO2 moves upward driven by buoyancy, a small part of the rock mass with low permeability is infiltrated slowly, and most of it bypasses the low permeability rock mass to find the rock mass with higher permeability for further migration. In comparison to Figures 2 and 9(a), the permeability of lithofacies 2 is the highest, and the distribution of is the largest in the area representing lithofacies 2. CO2 is primarily distributed in high-permeability rock mass after a 2-year injection, and its distribution decreases as permeability decreases; a substantial dominating flow phenomenon is shown.
Figure 9(b) shows that in the wholly heterogeneous model, the gas pressure is maximum near the injection well and gradually transfer to the saline aquifers as CO2 migrates. The lower the gas pressure, the further away from the injection well is.
Figure 10 displays the isoline map of gas pressure and the contour map of the distribution following a 2-year CO2 injection into horizontal wells in a deterministic multifacies heterogeneous model. The contour map of and the pressure transfer law in Figure 10 are comparable to those in the wholly heterogeneous model. As can be shown, the CO2 migration rule in the wholly heterogeneous model can be researched by conducting a numerical simulation by reducing the model to a lithofacies categorization model.
Figure 11 shows the contour map of and isoline map of gas pressure in a random multifacies heterogeneous model following a 2-year CO2 injection into horizontal wells. The distribution of and pressure transfer is generally consistent with the distribution of the contour map in the nearby area (Figure 10). The gas pressure near the injection well in Figure 11(b) builds higher, and the transmission speed is slower due to a modest change in the heterogeneity of the rock formations.
4.2. Physical Processes in CO2 Storage Based on a Heterogeneous Model
In the physical process of CO2 storage, the role of the different permeabilities of rock can be expressed as follows: (1) viscous fingering of CO2 along high-permeability fast flow pathways and (2) dynamic intrusion of CO2 from high-permeability zones into low-permeability zones by capillarity, as well as buoyancy . Taking the simulation results of the deterministic multifacies heterogeneous model as an example, we output contour maps of gas distribution (as shown in Figure 12) after 30 days, 100 days, 1 year, and 2 years of CO2 injection; the purpose is to observe whether the simulation process could reflect the physical process caused by different permeability of rocks in the real CO2 geological storage process.
When determining the classification of lithofacies according to the permeability value, the order of permeability is lithofacies 2, lithofacies 1, lithofacies 3, and lithofacies 4. Select the area near the well () and mark the area with high permeability, as shown in Figure 13.
In Figure 12, the orange arrows indicate physical processes (1) in CO2 storage, the green arrows indicate physical processes (2) in CO2 storage, and high-permeability areas are outlined by the magenta dotted line. CO2 preferentially occupies the high-permeability region (Figure 12(a)); it shows the physical processes (1) in CO2 storage: viscous fingering of CO2 along high permeability fast flow pathways; as the injection time goes on, gas CO2 gradually full of the high-permeability area and expands to the low-permeability area; the upward expansion trend is larger than the downward expansion; indicating the superimposed effect of buoyancy (Figures 12(b)–12(d)), it shows the physical processes (2) in CO2 storage: the dynamic intrusion of CO2 from high-permeability zones into low-permeability zones by capillarity, as well as buoyancy. Such physical processes cannot be captured in the simulation of the homogeneous model  or layer model but exit on a real storage site .
4.3. Comparison of CO2 Migration Laws in Three Heterogeneous Models
The distribution characteristics of gas CO2 are similar in Figures 9(a) and 10(a), but the high area in Figure 10(a) is more concentrated than in Figure 9(a), and the high area in Figure 9(a) is speckled (as shown in Figure 14). This is because the permeability values allocated to each grid in the wholly heterogeneous construction are spread, ranging from 0.1 D to 22 D, causing the distribution to be dispersed and the gas concentration to occur only in individual grids (Figure 9(a)). The deterministic multifacies heterogeneous model, on the other hand, is a simplification of the wholly heterogeneous model, which is classified into four values based on the permeability value, causing adjacent grid permeability values of different lithofacies to jump and adjacent grid permeability values of the same lithofacies to be the same value, resulting in the local gas concentration phenomenon (Figure 10(a)). This discrepancy, however, does not affect reservoir pressure (the comparison between Figures 9(b) and 10(b)). In terms of injection research, it is sufficient to reduce the wholly heterogeneous model to a lithofacies classification structure for numerical simulation to analyze the CO2 migration law.
When comparing Figure 10(a) with Figure 11(a), the distribution characteristics of CO2 are mostly the same, with high values in the two regions. The first area shows that CO2 is distributed at the top of the formation layer affectedby buoyancy; the second area shows that there is a local low permeability rock mass within the range of 100~ 150 m vertically and 500~ 800 m horizontally from the origin, which blocks CO2 migration and results in CO2 accumulation within 500 m of the injection well, forming another area of high (as shown in Figure 15). The unknown geological profile’s CO2 distribution characteristics are similar to those of the known geological section, indicating that CO2 migrated from the known geological section to the unknown geological section. As can be shown, the multifacies heterogeneous model constructed at random using the transfer probability theory could be used to simulate CO2migration in the lithofacies classification structure’s saline aquifers and to evaluate CO2 distribution features in unknown formations.
Based on the foregoing, it is sufficient to simplify the wholly heterogeneous model to lithofacies classification model for numerical simulation to research the CO2 migration rule when the hydrogeological data collected in the field are not comprehensive. When the hydrogeological parameters of the formations are unknown, the spatial distribution of multifacies permeability can be randomly produced using the transfer probability theory to characterize the formation’s heterogeneous properties based on neighboring known hydrogeological data. A random multifacies heterogeneous model reflecting the development of sedimentary rock masses is established to evaluate the CO2 migration in heterogeneous deep saline aquifers.
4.4. Comparison of Dry-Out Effect in Three Heterogeneous Models
The features of the dry zone induced by the dry-out effect following dry gas CO2 injection were investigated for the three types of heterogeneous models mentioned above. Draw the solid saturation () curve as a function of horizontal distance from the well, as shown in Figure 16.
Figure 16 shows that the precipitation of salt crystals decreases with increasing distance from the injection well, and the dry-out effect occurs mainly in the formation of the near-well area, and in the wholly heterogeneous model, the value of the wellbore is extremely high, which will reduce the efficiency of CO2 injection .
The formation hydrogeological characteristics, the distribution range of , and the value of sidewall in the near well area are listed in Table 2.
Table 2 demonstrates that the formation heterogeneity influences the dry-out effect in the near-well region, and that the more the heterogeneity, the greater the salt precipitation, and the shorter the CO2 drying impact influence distance. When comparing the results of the multifacies heterogeneous model with the results of the wholly heterogeneous model, it can be seen that when the classification of the wholly heterogeneous model is simplified, the original strong heterogeneity is reduced to local homogeneity, especially near injection wells, and the dry-out effect is eliminated (as shown in Figure 17). It is therefore required to adjust the grid and partition the lithofacies of the elements around the well when utilizing the random multifacies heterogeneous model to analyze the CO2 dry-out impact.
The deterministic multifacies heterogeneous model, the random multifacies heterogeneous model, and the wholly heterogeneous model are created to define the heterogeneous properties of deep saline aquifers based on the data available. The distribution and pressure transfer properties of CO2 following a 2-year injection were investigated. The following are the primary conclusions:
We offer the following strategy for constructing alternative heterogeneous models based on whether the amount of data is sufficient: firstly, when the in situ hydrogeological parameter data are more detailed, a wholly heterogeneous model can be established by assigning different hydrogeological parameters to each grid unit in the model; under limited data conditions, rock masses with the same hydrogeological characteristics can be classified into the same lithofacies and establish a multifacies heterogeneous model. When the hydrogeological parameters of the target formation area are unknown, the heterogeneity of the formation can be characterized by randomly generating the spatial distribution of permeability of multifacies based on the transfer probability theory and establishing a random multifacies heterogeneous model as close as possible to the field distribution characteristics of the formation
The migration of CO2 in heterogeneous deep saline aquifers shows a significant fingering flow phenomenon and reflects the physical process in CO2 storage. The completely heterogeneous model is simplified by the deterministic multifacies heterogeneous model. The distribution of gas CO2 produced using the two models is identical; it shows that simplifying the wholly heterogeneous model as a lithofacies classification structure is sufficient for numerical simulation of the CO2 storage. The CO2 distribution in the unknown formation segment is close to the known formation section; the multifacies heterogeneous model based on the transfer probability theory can be used to evaluate the distribution of CO2 in the unknown formations and is an effective way to explore the CO2 migration law of the unknown formation
In the simplified classification of heterogeneous models, the strong heterogeneity is simplified to local homogeneity, and the dry-out effect is weakened. It is necessary to refine the grid and subdivide the lithofacies of the local key area units when the spatial scale of the research problem is small
The method provides a feasible reference for heterogeneous modeling in the data-deficient region and the simplification of large-scale heterogeneous modeling works
All data generated or analyzed during this study are included in this published article and references.
Conflicts of Interest
The authors declare that there is no conflict of interest regarding the publication of this paper.
The research was supported by the Jiangsu Planned Projects for Postdoctoral Research Funds (No. 2020Z006).