Numerical Simulations of Carbon Dioxide Storage Efficiency in Heterogeneous Reservoir Models

. The U.S. Department of Energy ’ s National Energy Technology Laboratory (DOE-NETL) has been developing methods and tools (the online Carbon Dioxide Storage prospeCtive Resource Estimation Excel aNalysis (CO 2 -SCREEN) tool) to estimate carbon dioxide (CO 2 ) storage potential in subsurface reservoirs. The CO 2 storage e ﬃ ciency terms are input in the tool to calculate storage potential in targeted reservoirs. In this e ﬀ ort, two CO 2 storage e ﬃ ciency terms were evaluated: volumetric displacement ( E V ) and microscopic displacement ( E d ). The ﬁ rst term deals with e ﬃ ciency of CO 2 propagation into an accessible reservoir volume, while the second term evaluates e ﬀ ectiveness of native ﬂ uid displacement with CO 2 . The interpreted well logs and core sample measurements were applied to create the heterogeneous reservoir models including geostatistical realizations of porosity and intrinsic permeability ﬁ elds. Supercritical CO 2 was injected over the course of 30 years into brine-saturated reservoir models for clastics, limestone, and dolomite lithologies and deltaic ﬂ uvial, aeolian, shallow marine, and reef depositional environments by means of varying reservoir parameters and injection scenarios. The reservoir models providing vertically heterogeneous petrophysical properties and designated as “ layered reservoir models ” (with homogeneous parameters along each layer of the model) were not determined to be a transition between the homogeneous and heterogeneous models in respect to storage e ﬃ ciency. Another ﬁ nding shows that high-e ﬃ ciency factors do not necessarily mean increased CO 2 storage; they rather indicate that the available volume and pore space are more fully utilized. The CO 2 storage e ﬃ ciency factors were evaluated dynamically at the select time points using P 10 ‐ P 50 ‐ P 90 percentiles. The results of this study show that the P 10 ‐ P 90 distribution for volumetric e ﬃ ciency is wider when compared to the microscopic e ﬃ ciency. It was found that where dominant buoyancy forces drive the plume to the top of a target formation, the volumetric e ﬃ ciency is low. Tighter sandstone and carbonate formations show prevalence of capillary forces and better utilization of reservoir volume.


Introduction
Subsurface storage of CO 2 is gaining interest as a means for mitigating greenhouse gas emissions. There is a need for development of methods that estimate how much CO 2 can be stored in the subsurface to be widely available and easily accessible. Recently, a standard guide for estimating storage of CO 2 in the subsurface with consistent terminology was developed by the Society of Petroleum Engineers (SPE) and the American Association of Petroleum Geologists (AAPG) called the Storage Resources Management System (SRMS) [1]. This system focuses on developing a framework for large-scale commercial storage of CO 2 in the subsurface where a project is clearly defined with active injection wells and the amount of storage space available is quantified. There are three storage classifications under the SRMS: storage capacity, contingent storage resources, and prospective storage resources. Capacity quantifies storable quantities for developed projects actively injecting CO 2 , while both contingent and prospective storage resources cover projects that are not proven but include various degrees of regional assessments. The U.S. Department of Energy National Energy Technology Laboratory (DOE-NETL) has worked to develop storage methods at the contingent and prospective storage levels for saline formations, unmineable coal seams, depleted unconventional shale formations, and residual oil zones [2][3][4][5]. An online version of these methods is available as the CO 2 Storage prospeCtive Resource Estimation Excel aNalysis (CO 2 -SCREEN) tool [3]. These methods and tools have been used worldwide to estimate contingent storage and prospective storage and incorporated into best practice field manuals compiled from CO 2 injection tests at various pilot-scale sites [6][7][8][9][10]. The introduction of CO 2 storage efficiency factors utilized in the CO 2 -SCREEN tool is provided below.
The U.S. DOE is considering geologic carbon storage as a means for reducing anthropogenic emission of CO 2 that contributes to climate change and formed a nationwide network of regional partnerships to help determine the best approaches for capturing and permanently storing CO 2 [10]. The U.S. DOE Carbon Storage Program also developed a CO 2 storage method in conjunction with the Regional Carbon Sequestration Partnerships (RCSP) for the estimation of volumetric CO 2 storage in deep saline formations, oil and gas reservoir, and unmineable coal seams in the U.S. and portions of Canada [6][7][8][9]. The RCSP were government/ industry efforts tasked with determining the most suitable technologies, regulations, and infrastructure needs for carbon capture, storage, and sequestration in different areas of the country [10]. Saline aquifers show high potential in terms of storage potential and safety as evidenced by development phase field projects for CO 2 storage [10].
As reported in Ref. [11], the U.S. DOE-NETL methodology introduces a volumetric method to estimate prospective mass CO 2 storage resource (G) using where A d , h s , ϕ s , and ρ s are the areal size of the formation, the thickness of the formation, porosity, and CO 2 density (estimated at average pressure and temperature of the storage formation), respectively. Superscripts d and s indicate that a corresponding parameter is estimated deterministically and stochastically (the superscripts are further omitted for clarity). The storage efficiency (E saline ) term reduces the estimation of stored CO 2 mass at a specific site to accommodate the complexities of the fundamental processes associated with production, injection, and storage within saline aquifers.
where E A , E h , and E ϕ are the fraction of the geologic area, thickness, and porosity accessible for CO 2 storage, respectively. Since the individual efficiency terms in Equation (2) are determined stochastically, Equation (2) is required to be written in log odds notation and expanded to Equation (3) (as described in Ref. [9]) to determine storage efficiency (E saline ) using Monte Carlo sampling. E V , the volumetric displacement efficiency, represents the fraction of reservoir volume accessed by CO 2 plume. It is calculated as the ratio of CO 2 injected to the pore volume of the reservoir accessed by the CO 2 plume.
where V i is the volume of injected CO 2 ; A, h, and ϕ are the area, thickness, and porosity of the accessed volume, respectively; and S w irr is the irreducible water saturation. The volume of the injected CO 2 can also be expressed through mass flow rate ðQ i Þ, injection duration ðtÞ, and averaged CO 2 density (ρ) at in situ conditions (Equation (4)). To determine the area, A, the minimum-area-circle approach is used as depicted in Figure S1 (Supplemental Materials). The accessible volume is the multiplication of the area and thickness of the formation. E d , the microscopic displacement efficiency, is the fraction of water displaced by CO 2 , which is calculated at the plume level by the average water saturation (S w ave ) as shown in Equation (4) [12].
where S CO 2ave is the average CO 2 saturation within the CO 2 plume. It is assumed that CO 2 would not be able to occupy all effective porosity defined by ð1 − S w irr Þ [12]. The volumetric and microscopic displacement efficiencies are dynamic terms in Equation (2) and Equation (3) since their values are time-dependent and change with the CO 2 plume propagating through a target formation. As injection progresses with time, plume shape changes and, therefore, E V changes. In addition, E d changes as average saturation values change due to buoyancy forces, capillary trapping, and other processes. The time dependency of storage efficiencies emphasizes the importance of the heterogeneous representation of porosity and intrinsic permeability in reservoir models mimicking actual pore network in a saline aquifer and coupled permeability controlling mobile phase flow.
Shao et al. [13] performed simulations on various stochastic model realizations under different conditions and showed that nested geological heterogeneity is a major player in controlling plume propagation. Rasheed et al. [14] found that an aquifer with a porosity higher than 20% and a low to medium level of heterogeneity is more suitable for CO 2 storage. Stochastically generated and spatially correlated permeability distributions have shown that heterogeneous permeability significantly affects both the accumulation and distribution of CO 2 [15][16][17]. Heterogeneity also promotes the local capillary trapping that may deviate the path of rising CO 2 [18,19]. Both experimental and numerical studies have identified the great impact of 2 Geofluids heterogeneous capillary trapping on the performance of CO 2 storage [20,21]. Several studies have explored CO 2 storage in saline formations using numerical simulation [22]. Storage potential and efficiency factors have been investigated by developing heterogeneous models and performing sensitivity analysis using geologic (e.g., ratio of net to gross thickness, ratio of net to total area, and ratio of effective to total porosity), macroscopic (e.g., gravity), and microscopic (e.g., saturation) terms [23]. To estimate storage efficiency for specific lithologies and depositional environments, heterogeneous geocellular models and geological architecture common to each of the environments had to be developed [24][25][26]. However, prior storage efficiency values were mostly based on limited geologically nonspecific relative permeability data.
Recently, CO 2 storage capacity in deep saline aquifers and residual oil zones (ROZs) was investigated using machine learning algorithms [27,28]. The developed smart tools could be used to provide fast and accurate storage efficiency for aquifers that have similar parameters falling within the range of the database [27]. The proposed artificial neural network models can predict the CO 2 -enhanced oil recovery (EOR) and storage performance with high accuracy in ROZs [28].
In our recent work, storage efficiency parameters were estimated based on homogeneous reservoirs [29]. Homogeneous CO 2 storage efficiency values for supercritical CO 2 injection into brine-saturated reservoirs for clastics, limestone, dolomite lithologies and marginal marine, strand plain, deltaic complex fluvial, aeolian, shallow marine, and reef depositional environments were determined for 30 years of injection [29]. In this work, we add complexity to the storage efficiency factors by incorporating heterogeneity into the reservoir models. The numerical simulations of fluid dynamic models to estimate volumetric displacement (E V ) and microscopic displacement (E d ) efficiency factors for saline formations were performed using the TOUGH3 code [30]. In this effort, heterogeneity was accounted for in the reservoir models designed for CO 2 injection using downhole well log information retrieved at pilot sites of different lithologies and depositional environments. The models are developed using the corresponding experimental relative permeability data available in the CO 2 -Brine Relative permeability Accessible (CO 2 BRA) database [31]. The heterogeneous storage efficiency values are also available in DOE-NETL's CO 2 -SCREEN tool v4.1 [3].

Methodology
The workflow of the manuscript includes creating heterogeneous reservoir models that represent five combinations of lithology and depositional environment. Next, simulation cases are run with varying pressure, temperature, permeability anisotropy, injection rate, reservoir thickness, porosity, and permeability distributions. Then, the results of CO 2 injection over 30 years are processed to calculate timedependent storage efficiency terms using Equations (4) and (5). The computed efficiency factors are then analyzed and compared with factors estimated using homogeneous reser-voirs. A summary table listing the storage efficiency factors as P 10 -P 50 -P 90 percentiles is provided at the end of the manuscript.

DOE Carbon Storage Partnership Sites Used as Proxy for
Heterogeneous Reservoir Models. Five saline reservoir models and five depositional environments were considered in this study based on the data collected by the RCSP. The heterogeneous reservoir models were based on the corresponding sections of saline aquifers showing the best quality in terms of total porosity and intrinsic permeability for continuous CO 2 injection. Below is the list of the formation/depositional environments utilized for making the models and their brief descriptions. A suite of wireline logs was collected at each site and interpreted with core sample measurements to provide total porosity and intrinsic permeability distributions within the reservoir thickness. Supplemental material contains derived porosity and permeability logs used to create the heterogeneous reservoir models ( Figures S2, S5, S7, S10, and S13).  Besides heterogeneity of the petrophysical properties in a target formation, another strong factor influencing the efficiency terms is relative permeability (k r ) controlling multiphase flow in pore space. The impact of different k r relationships and their parameters is critical to accurately predict the propagation supercritical CO 2 (scCO 2 ) and saturation distribution in a reservoir [38,39]. In this study, CO 2 BRA, a freely available database containing experimental data for rock types of different depositional environments [31,40], was utilized to correlate lithology and depositional environment of the select heterogeneous reservoir models with k r experimental data measured for corresponding core samples. The k r data for the corresponding lithology/deposition were made available in the TOUGH3 code [30] used to conduct simulations. Details of the CO 2 BRA database implementation into the code and selection of experimental k r curves for a lithology/deposition environment of interest can be found elsewhere [29]. Table 1 lists the five target formations and their dominant deposition environments, described in the previous section, against relevant CO 2 BRA sample names of the same lithology and deposition. The table also collects porosity, permeability, and parameters of residual saturations and maximum k r values derived from the unsteady-state experiment of CO 2 injection into the brine-saturated core samples [40]. Noticeably, there are low endpoint k r values of CO 2 (<0.3) and high residual brine saturations (>0.4) for the k r curves that were typical for the unsteady-state method of scCO 2 injection [41,42]. The method relies on the injection of scCO 2 to displace brine in a fully saturated sample and never truly settles on a steady-state condition due to the incremental displacements of the primary fluid [43,44]. The low endpoint numbers for scCO 2 and high residual saturation values for brine are attributed to capillary end effects, sample heterogeneity, and low mobility ratio by Jeong et al. [45].
Supplemental material contains descriptions of sidewall core samples (where available) extracted at the geological wells and the CO 2 BRA samples corresponding to similar lithology and depositional environment ( Figures S3, S4, S6, S8, S9, S11, S12, S14, and S15 and Tables S1-S2). The images of the samples are accompanied with general  To populate the mesh with porosity and permeability data, the well log information depicted in Figures S2, S5, S7, S10, and S13 and core measurement correlations (Supplemental Materials) were used. Two types of heterogeneity were introduced into the reservoir models: (1) layered. Porosity and permeability values were varied with depth following the mesh discretization in the vertical direction. To do so, the corresponding values were upscaled (red lines in the above-mentioned figures in Supplemental Materials) to fit layers of the meshes. In the lateral direction, the properties were assigned uniformly assuming they are laterally homogeneous. This model implies ideal geological control of the properties in the horizontal directions. (2) Heterogeneous. For those models, the density porosity logs across the reservoir thickness were utilized to estimate empirical variograms in the vertical direction. Since only one well log was adopted for each formation, to quantify spatial correlation in the lateral direction, the anisotropy factors were applied following the guidance from Gorecki et al. [23] for horizontal variogram ranges or information about geologic control of bedding planes at nearby wells if available (like for Danielson and Wallewein wells drilled into the Middle Duperow formation) [46]. Since the reservoir models utilize a 2D approximation, no azimuthal control for model variogram was used. That implies the same geospatial property correlation in the horizontal direction. Using the model variogram as weighting functions, ordinary kriging along with Sequential Gaussian Simulation was performed using the normal score-transformed porosity as the input data on a rectangular reservoir grid. The output from the stochastic simulation, which is a normal score porosity, was transformed back to actual porosity values. A geostatistical software, WinGSLib v.2015 [47], was used to calculate porosity distributions for each heterogeneous reservoir model. To generate coupled intrinsic permeability data, the correlation derived from permeability logs, or core measurements, was used to develop permeability values. Up to nine realizations of porosity and coupled permeability distributions were generated for use in simulations to account for uncertainty in property distributions. Figure 2 shows the layered and heterogeneous representation of the sandstone/fluvial reservoir models generated using the Cranfield Sandstone formation as a proxy.
The initial pressure and temperature conditions were selected based on data from Gas Information System (GASIS) [48], which is a public database containing approximately 20,000 oil and gas reservoir records. The initial P-T values correspond to CO 2 in supercritical condition. CO 2 injection was executed for 30 years using an open interval across the entire reservoir thickness. Table 2 collects pertinent reservoir model parameters and injection scenarios. The simulation cases designed to estimate statistical distributions of the CO 2 storage efficiency terms are considered in the next section. Table 1, simulation cases were designed by varying five parameters for sensitivity analysis, as shown in Table 3.

Simulation Cases. For each combination of lithologies and depositional environments listed in
They included pressure, temperature, permeability anisotropy (i.e., the ratio of vertical permeability to horizontal permeability, k v /k h ), reservoir thickness, and injection rate. A min/max value was selected for each parameter. The min/max pressure was coupled with min/ max temperature representing shallow/deep formations, respectively. The minimum pressure and temperature were selected in such a way to assure the existence of scCO 2 . The min/max for permeability anisotropy was set at 0.1 and 0.5 [52,53]. The min/max injection rates were selected in such a way to be consistent with field data [32][33][34][35][36][37] and the maximum pressure buildup around the wellbore and formation fracture pressure [29]. The thickness variability was based on consideration of typical perforation intervals  [34].
In this study, the numerical simulations were carried out using the ECO2M module of TOUGH3 [30] designed for simulations of CO 2 storage in saline aquifers. The module represents a comprehensive description of thermodynamic properties of a H 2 O-NaCl-CO 2 mixture [54] within broad ranges of pressure, temperature, and salt concentrations. The module accounts for phase changes that can arise during CO 2 injection involving gaseous, liquid, and supercritical phases of CO 2 together with brine and has improved convergence during the appearance and disappearance of a nonaqueous phase. PetraSim [55] was used to postprocess the results of the simulations. Table 3 lists 16 total cases created using various injection scenarios and reservoir parameters, which together with nine geostatistical realizations of porosity and coupled permeability distributions for each case provide 144 simulations carried out for each reservoir model type: layered and heterogeneous models. In total, 288 simulation runs were conducted. The results of the simulations were postprocessed following the approach described in Section 2 to calculate CO 2 efficiency terms and corresponding percentile values (Equations (4) and (5)). Figure 3 compares volumetric displacement efficiency (E V ) and microscopic displacement efficiency (E d ) computed using layered and heterogeneous reservoir models with homogeneous values calculated in the previous study [29] for select lithology/depositional environment. The efficiency terms are represented as horizontal bars (probability distribution ranges) where the minimum and maximum values correspond to P 10 and P 90 values (10 th and 90 th percentile), respectively. The table also contains minimum and maximum porosity and permeability values extracted from the GASIS database [48] for homogeneous models and average porosity/permeability values evaluated using layered and heterogeneous (one realization) reservoir models.

Results and Discussion
The analysis of Figure 3 indicates that introducing geologic heterogeneity into petrophysical properties strongly   Reservoir thickness * For sandstone reservoir models, the Brooks-Corey capillary pressure function was used using the parameters from Refs. [49,50], and for limestone and dolomite reservoir models, the van Genuchten capillary pressure function was utilized using the parameters from Ref. [51]. Leverett scaling of capillary pressure was invoked since heterogeneous porosity and permeability fields were used. 6 Geofluids influences flow characteristics determining the plume size and shape. This translates into broad variations of the volumetric efficiency that range from 11.5% in sandstone/aeolian to 90.5% in carbonate/shallow marine depositional environments. In some cases, probability distribution ranges of volumetric efficiency become narrower by adding heterogeneity (e.g., sandstone aeolian), and in others, ranges become wider (e.g., carbonate reef limestone). Given that the volumetric efficiency is based on a CO 2 plume shape and maximum propagation of CO 2 into pore space of a saline aquifer, it is expected that the volumetric efficiency is sensitive to porosity and permeability distributions in a reservoir. The results indicate that there is no trend in probability range changes moving from homogeneous to layered and then to heterogeneous models. We also concluded that layered models utilizing idealistic property distribution in the lateral direction cannot be considered as "transition" models from homo-to heterogeneous models in terms of volumetric efficiency evaluation. Since layered models use only vertical variability in porosity and permeability, locations of layers with high values of porosity and permeability in the horizontal direction ( Figures S2, S5, S7, S10, and S13) facilitate plume extension into those layers and determine an efficiency value. Microscopic efficiency ranges are generally less impacted and demonstrate narrow ranges compared to the distribution ranges of the volumetric efficiency. Only the carbonate/reef data show wide distribution ranges for both layered and heterogeneous models. For that combination of lithology and deposition, large shifts are observed in both microscopic and volumetric efficiency ranges, moving from homogeneous to layered and then to heterogeneous models. These variations of CO 2 plume shape and storage efficiencies are further explained in Figure 4. The presence of high permeable layers ( Figure S13) in the middle of the layered model (Figure 4(b)) allows for CO 2 to channel through the formation, resulting in a different plume shape compared to the homogeneous or heterogeneous models (Figures 4(a) and 4(c)). As a result, volumetric efficiency is reduced because of the extended plume shape for the layered model. Reduction of microscopic efficiency of the layered model is associated with diminishing of CO 2 saturation in the areas behind the advancing CO 2 front. Such behavior indicates dominance of capillary trapping over transverse buoyancy crossflow and is typical for layered systems of contrasting porosity/permeability [56]. The heterogeneous model provides a larger E V value when compared to the homogeneous model (Figures 4(a) and 4(c)). In the heterogeneous case, the efficiency term is calculated based on the minimum-area-circle approach ( Figure S1) and near vertical plume interface which results in more effective reservoir volume utilization.
Thus, the balance between viscous, capillary, and buoyancy forces affects the efficiency terms. A dimensional analysis has been widely used to characterize multiphase fluid flow in porous media to quantify the ratio of a certain force over another one [57][58][59][60]. The residually trapped nonwetting phase is directly affected by a tradeoff between viscous, capillary, and buoyancy forces and was shown to be correlated with dimensionless number magnitude, like capillary and bond numbers [61][62][63][64][65]. To analyze dominant forces acting during CO 2 flow in porous media of the reservoir, the dimensionless bond (B o ) number was engaged which assesses the ratio of gravitational (negative to buoyancy) and capillary forces (Equation (6)) from Shook et al. [66], Kuo and Benson [67], Trevisan et al. [68], and Ben [69].    Figure 3: Impact of heterogeneity on probability distribution of the efficiencies at 30 years of CO 2 injection. The values of the homogeneous models were obtained from Haeri et al. [29].  where Δρ, σ, and θ are the density difference of two fluids, interfacial tension, and contact angle between supercritical CO 2 and brine at storage conditions, respectively, g is the gravitational constant, and k V is the average vertical permeability of the formation.   9 Geofluids pressure, and temperature conditions is depicted in Figure 5. To calculate the B o number, the density difference between CO 2 and brine at storage conditions and permeability are averaged over the grid cells within the plume area. The average values of 27 mN/m [70][71][72] and average 22° [73,74] were used for the interfacial tension and contact angle, respectively. As a general trend, the larger the B o number, the more tendency of CO 2 to move upward reflecting the dominance of buoyancy force. In the case of the sandstone/aeolian formation (B o = 2:6), the high permeability facilitates migration of CO 2 upward under the influence of density difference and its accumulation at the boundary with the seal acting as a structural trapping. On the other hand, the B o numbers smaller than 1 indicate the dominance of capillarity resulting in retarded vertical movement and more capillary trapping in horizontal direction, as seen in Figure 4. Consequently, the B o number serves as a criterion to determine the dominant force for flow that correlates with the resultant plume shape (Figure 4).
Both E V and E d are impacted by variation in the parameters listed in Table 3. However, E V was influenced the most. Figure 6 illustrates the change in the probability distribution of E V (P 10 -P 90 range) with variation of the modeling parameters, including pressure and temperature (associated with depth), permeability anisotropy, injection rate, and reservoir thickness in different depositional environments. The light and dark bars of each color correspond to minimum and maximum values of a certain parameter, respectively. Raising the pressure and temperature (i.e., going deeper in the formation) consistently shifts the efficiency distribution range upward, showing an increase in the volumetric efficiencies ( Figure 6, the first plot from the top). The deep formations typically are more favorable for CO 2 storage as higher temperature associated with viscosity reduction makes the CO 2 more mobile [60]. Gorecki et al. [23] demonstrated that with increase of temperature, there is a small increase in the microscopic efficiency and a small decrease in volumetric efficiency. The impact of temperature is not studied independently in this work. However, pressure is shown to be the leading factor compared to temperature, resulting in a consistent increase in the volumetric efficiency ranges. The impact of injection rate ( Figure 6, second plot from the top) is similar to pressure and temperature influence. Using a high injection rate, a larger CO 2 mass enters the reservoir for the same injection time compared to cases with a low injection rate. Consequently, capillary entry pressure is likely overcome, leading to diminishing capillary trapping and accessing more pore volume and thus higher volumetric efficiency.
The influence of permeability anisotropy is the opposite (Figure 6, third plot from the top). Increasing the permeability anisotropy (changing K v /K h from 0.5 to 0.1) shifts the 10 Geofluids volumetric efficiency ranges downward (i.e., reduction in volumetric efficiency). With more contrast permeability anisotropy, CO 2 moves predominantly in the lateral direction that leads to a decrease in volumetric efficiency. In other words, the mechanism of CO 2 fluid flow is shifted from buoyancy-driven to permeability-controlled by the predominant flow in the horizontal direction. An exception is the sandstone/aeolian case, in which the average permeability is high enough to compensate the lat-eral flow advantage with transverse buoyancy flow, leading to volumetric efficiency being insensitive to the variation of permeability anisotropy. To illustrate the influence of permeability anisotropy, Figure 7 depicts the plume shapes for the sandstone/shallow marine case with a relatively more drastic E V reduction under the influence of an increase in permeability anisotropy. The reduction of the B o number correlates with more extensive plume propagation in the lateral direction under the dominance of capillary forces, when

Geofluids
K v /K h = 0:1. The capillary trapping makes the propagation front uneven with multiple areas of elevated CO 2 saturations.
Increasing the reservoir thickness leads to decreasing both P 10 and P 90 values for sandstone/shallow marine and carbonate/reef limestone. For the remaining lithology/deposition, the changes of P 10 and P 90 with thickness increase occur in different directions ( Figure 6, bottom plot). That means that thickness may be strongly coupled with other reservoir parameters and implemented heterogeneity. It should be noted that in most cases, the trends of the parameter impact on P 10 and P 90 values remain preserved over the course of 30 years of injection.
As an example, Figure 8 displays dynamic changes of the efficiency values with time for carbonate/reef limestone. The P 10 and P 90 values of the higher thickness remain under corresponding storage efficiencies of the lower thickness over the course of 30 years of injection.
The effect of the modeling parameters and injection rate on E d is less pronounced compared to E V . Figure 9 shows the microscopic efficiency terms affected by those parameters and the rate. There is a general increasing trend in efficiencies with increasing depth (higher pressure and temperature). Sim-ilar trend is found for the volumetric efficiency that agrees with the previous work using the homogeneous reservoirs [23,29]. The permeability anisotropy contrast (changing K v / K h from 0.5 to 0.1) undermines the ability of CO 2 to displace brine and, thus, decreases the microscopic efficiency (similar trend with the volumetric efficiency) in all the selected lithologies and depositional environments. This finding is different than what was found in the previous work [29] that showed little sensitivity of the efficiency terms on anisotropy after a few years of injection using homogeneous reservoir models. There is a general trend of increase of the P 90 value for E d with an increase of injection rate while effect of thickness shows no certain trend across all cases. That implies that for a reservoir with dominance of one lithology and depositional environment, the thickness of a formation plays little role in microscopic efficiency. Table 4 collects dynamic E V and E d efficiency terms and their product calculated using P 10 , P 50 , and P 90 values for select time points: 1, 5, 10, 20, and 30 years. It is noticeable that efficiencies are time-dependent. After 1 year of injection, the P 10 , P 50 , and P 90 values are the lowest, and then, they increase toward 5 years and get relatively flat after 10  Table 4 are transferred to the CO 2 -SCREEN tool [3], and updated values are available in the online version of the tool.
In should be noted that since the efficiency terms are estimated up to 30 years of injection, reactive transport of CO 2 in the formations and the effects of geochemical reactions were not included. Although it is known that CO 2 can significantly impact pore structure, especially in carbonate reservoirs [41,75], the data show that mineral dissolution is most likely to develop during hundreds of years (except carbonates) [76]. Similarly, since postinjection periods and possible imbibition of the displaced fluid were not accounted for, the relative permeability hysteresis effects were not included into the simulations. The estimated heterogeneous efficiency terms are designed to be applied for formations with dominant lithologies and depositional environments listed in Table 1. For a target reservoir with large variability in depositional environment within stratigraphic layering, those factors should be applied with caution. Distinct structural features like dome, anticline, and incline were not considered in this work. It should be kept in mind that the reported efficiency factors can be influenced by such structural traps and acting gravitational forces. The 2D approximation used to create the reservoir models implies that the CO 2 plume propagates uniformly from a wellbore. However, geospatial property correlation may vary in the azimuthal directions that together with possible presence of sealing faults deviate CO 2 plume expansion from its symmetrical shape. Consequently, the uncertainty in efficiency terms due to an asymmetric plume shape will increase with injection time. For our next effort, we consider studying CO 2 injection using 3D heterogeneous reservoirs and developing a methodology of efficiency terms estimation for multiple well injection scenarios.

Conclusions
The results of this work indicate that there is a strong effect of the reservoir parameters and injection rate on the volumetric efficiency, while microscopic efficiency is less influenced as it is evident by narrow ranges of P 10 -P 90 values for E d compared to the corresponding ranges for E V . The increase of pressure and temperature and injection rates results in a general trend of increasing efficiencies. Increasing permeability anisotropy decreases both efficiency terms contrary to prior work utilizing homogeneous reservoir models that showed little dependency. The data suggest that there is no clear trend in the efficiency changes because of using homogeneous, layered, and heterogeneous reservoir models. Layered models providing laterally homogeneous petrophysical properties cannot be considered as transition between the homogeneous and heterogeneous models in respect of storage efficiency change.
There are no consistent trends in P 10 -P 90 ranges comparing homogeneous and heterogeneous model results. The volumetric and microscopic efficiencies are sensitive to injection time. There is a trend of increasing efficiencies with time indicating better utilization of reservoir volumes.
The heterogeneous reservoir models using high porosity and permeability had the lowest volumetric efficiency. That was attributed to dominance of buoyancy forces leading to poor utilization of reservoir volume by the plume. For other heterogeneous reservoirs, the higher contribution of capillary forces results in better E V values (expressed through P 10 and P 90 values). It is shown that tighter reservoirs with low permeability and porosity demonstrate higher E V and E d , implying more efficient volume and pore network utilization. In other words, higher efficiency factors do not mean that more CO 2 can be placed in a formation, rather that the available volume and pore space will be more fully filled. A low porosity reservoir with a high-efficiency factor might hold less CO 2 than a high porosity reservoir with a lowefficiency factor.

Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.

Disclosure
This research was supported in part by an appointment to the National Energy Technology Laboratory Research Participation Program, sponsored by the U.S. Department of Energy. Neither the United States Government nor any agency thereof, nor any of their employees, nor the support contractor, nor any of their employees makes any warranty, expressed or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or usefulness of any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights. Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise does not necessarily constitute or imply its endorsement, recommendation, or favoring by the United States Government or any agency thereof. The views and opinions of authors expressed herein do not necessarily state or reflect those of the United States Government or any agency thereof.

Conflicts of Interest
The authors declare that they have no conflicts of interests.

Supplementary Materials
Supplementary Material includes graphics showing an approach to determine the accessible rock volume, porosity and intrinsic permeability logs and their average value distributions for reservoir models of select lithologies and depositional environments, and images of the samples with geological descriptions. (Supplementary Materials)