Reservoir modeling for Wabamun lake sequestration project

A study, entitled “Wabamun Lake Sequestration Project” or “WASP,” was performed to evaluate large‐scale CO2 storage opportunities in the Wabamun area including potential risks. The project examined the feasibility of storing 20 megatons (Mt) of CO2 per year over 50 years. This scale is one order of magnitude larger than the typical benchmark (1 Mt/year) used in academic research and commercial projects that are currently in place or under review. The study was conducted by a group of researchers from several universities as well as industry consultants. This study presents an overview of the reservoir modeling part of this study, which the authors were responsible for. The main objectives of the reservoir modeling for WASP project were as follows: (1) estimation of storage capacity (traditionally, this value is projected based on the available pore space, but we have an additional practical consideration: the maximum amount one can inject within short period of time [~50 years] and within a localized injection area [~30 km × 90 km]); (2) investigation of CO2 plume movement and pressure distribution during and after injection including the effect of formation dip angle on the plume shape and its migration; (3) investigation of the long‐term fate of injection associated with free phase CO2 (risks of leakage) and aquifer pressurization (possible geomechanical changes and related phenomena); (4) investigation of the phase behavior of H2S initially available and dissolved in brine during CO2 sequestration process. The WASP reservoir modeling study mentioned above led to a few important findings. The most important one is that, when CO2 is being injected into a sour aquifer, initially dissolved H2S will release into the expanding CO2 plume and accumulate at the leading edge of the plume. Also, the large‐scale injection scheme (20 Mt/year), which requires multi well injectors, provides very different pressure response compared to a one well (1 Mt/year) scenario.


Introduction
Carbon dioxide emissions from use of fossil fuels are likely to be the dominant drivers of climate change over the coming decades [1]. The use of carbon capture and storage (CCS) technology offers the possibility of maintaining access to fossil energy while reducing emissions of carbon dioxide to the atmosphere. However, for this technology to play a significant role in managing global emis-sions, it should be carried out at very large scales. Although aquifers offer a large storage potential for CCS, the injection of large volumes to fill up this capacity can be challenging. In this study, we performed a comprehensive numerical simulation of injection into the Nisku aquifer located in Wabamun Lake Area, Alberta, Canada, where large CO 2 emitters are present (Fig. 1).
Together the four coal-fired power plants in the Wabamun study area emit~30 Mt CO 2 per year. The Alberta Geological Survey has identified this area as a potential storage site with different underground storage targets. Identification of the targets was based on evaluation of a variety of parameters such as stratigraphy and lithology, fluid composition, rock properties, geothermal, geomechanical, and pressure regimes. One of these targets, the Nisku Formation in the Devonian Winterburn Group, was selected for WASP study. The Nisku was deposited at the edge of a carbonate shelf. From southeast to northwest, relatively pure platform carbonates change into interbedded limestone and shale of ramp and ultimately basin slope characteristics. The depth to the top of the Nisku formation, and for its whole areal extension, ranges between 1600 m in the northeast and gradually increases to 2150 m in the southeast at a dip angle of 0.5° [2]. This amount of change in the formation depth can consequently bring about a wide range of variation in fluid properties such as temperature and pressure. However, in the limited area of the injection site considered for this project, these variations are minimal and therefore neglected. Table 1 summarizes the average aquifer and fluid properties used in the study.
Both compositional and black-oil models have been reported in the literature to capture the interactions between CO 2 and brine in storage processes. The compositional modeling relies on equation of state (EOS) flash calculations, and therefore they are usually computationally very expensive and impractical for large simulation models. On the other hand, in the black-oil approach, PVT (pressure-volume-temperature) calculations are facilitated through simple table-look up of properties (e.g., fluids density, viscosity, and compressibility) which change solely as a function of pressure. Therefore, the "black-oil" approach is much faster and less vulnerable to convergence issues, which both are crucial for simulation of large models. Accurate black-oil PVT model for a twocomponent CO 2 -brine system have been developed and validated against the experimental data and were adopted in this project. The details of the PVT model and its implications can be found in Hassanzadeh et al. [3,4].   Figure A1(A-C) in Appendix A1 shows the properties of the fluid system pertinent to this study. The characteristic relative permeability curves of the Nisku carbonate, which have been measured at in situ conditions of pressure, temperature, and brine salinity are shown in Figure 2 [5].
Although the presented curves belong to a drainage process, the same curve is assumed for the imbibition process. Nevertheless, one can investigate the effect of relative permeability hysteresis on the residual gas trapping using different curves for imbibition and drainage processes. In addition, the effect of capillary pressure in all simulation models was neglected.
An important rock property, which is crucial to any injection scenario, is the fracture pressure gradient and hence the fracture pressure along the injectors wellbore. Practically, the injection pressure should be limited to 90% of the local fracture pressure to avoid the formation and propagation of fractures from the aquifer into the cap rock [2]. According to Bell et al. [6] the fracture pressure gradient for Wabamun area could be as high as 24 kPa/m; however, Karsten et al. [2] have used a more conservative value of 20 kPa/m. Based on these values, an average base case value of 22 kPa/m was assumed in study which corresponds to 40 MPa of fracture pressure at the average depth considered. The maximum bottomhole pressure of the injectors in the simulation models are constraint to this value. This constraint will ascertain that the pressure in the aquifer will be remaining in the safe operating window.
In the following sections, the results of reservoir simulation of saturation (plume size) and pressure distributions in a simplified conceptual homogeneous model, as well as full size homogeneous and heterogeneous reservoir models, are presented. Secondly, the effect of reservoir dip on plume evolution is discussed. Thirdly, the longterm fate of average reservoir pressure and CO 2 in the free gas phase is investigated. Finally, the phenomena of H 2 S exsolution during CO 2 injection into a sour aquifer are briefly addressed.

Sensitivity of Simulation Results to Grid Size
The numerical simulations in the current project were conducted using IMEX, the CMG commercial black-oil simulator. Since the horizontal extension of the presented models can reach several hundred kilometers, larger grids were used to construct the overall parent grid and local grid refinement (LGR) was applied in the vicinity of the injectors where the maximum changes in flow parameters happen. The LGR option in IMEX is a very practical way to quickly examine the effect of grid block size, especially when the dimensions of the model are very large. In Appendix B, more details about the grid size and grid refinement applied in simulation models are provided.

Preliminary Conceptual Model
To develop a benchmark at the beginning of the project, a simplified conceptual model was developed assuming homogeneous properties (provided in Table 1) and an infinite acting aquifer. A square (200 km 9 200 km) simulation domain was chosen to represent the aquifer (the results were not sensitive to an increase in model size to 250 km 9 250 km). By setting the model to these dimensions, the aquifer behaves as though it were infinite acting for the injection of the target volume of CO 2 . Figure 3 shows the model configuration for different numbers of vertical injector wells (1-25).
All wells are perforated from the top to the bottom of the aquifer. The number of wells (n) and placement configuration were designed to allow the use of an element of symmetry and hence reduce the total number of grid cells by a factor of four. The distance between the wells in x and y directions are the same and equal to k. The total cumulative amounts of injected CO 2 (denoted as V1, V4, V9, V16, and V25), increase with the number of wells and these amounts are split equally between injectors in each case. For example, in the case of nine wells, the yearly flow rate per each well is V9/(nine wells 9 50 years). In all cases, the injection time is chosen as 50 years and injection capacity is defined as the maximum amount, which is possible to inject without exceeding the average reservoir fracture pressure of 40 MPa.
Simulation results demonstrate that, for a single (n = 1) injector completed in the whole formation interval and at the center of the reservoir, the plume radius reaches 4.6 km after 100 years. Most plume propagation  takes place during the injection period and virtually stops after 100 years from the start of injection (50 years after injection stops). The injection rate capacity was found to be~1 Mt/year for this vertical injector (constrained by fracture pressure), while horizontal well 1 km in length and completed open-hole, can improve the injectivity tõ 1.5 Mt/year. For the multiple wells (n > 1) injection scenario, CO 2 saturation plumes do not interfere and n individual plumes each having a radius 4-5 km are observed. However, the pressure field behavior is completely different from the saturation field. By the end of CO 2 injection period, there are no individual pressure plumes. Instead, different pressure plumes merge into a single large (scale of hundred km) pressure envelope. Injection capacity increases with the number of wells but this effect is rapidly vanishing after n~10. These phenomena will be discussed in detail later, when a full Nisku aquifer geomodel is considered.

Detailed Modeling with Full Aquifer Extent
In this section, simulation results for CO 2 injection into the Nisku formation are presented. First, a homogenous model is used to investigate injectivity into a semiinfinite formation. Then, a heterogeneous model is populated with realistic permeability and porosity fields in order to demonstrate the effect of heterogeneity and reservoir dip angle on the evolution of a CO 2 plume as well as the impact on reservoir pressure. Figure 4 shows a top view of the Nisku aquifer. The illustrated area is approximately 450 km 9 640 km while the area outlined in red corresponds to the injection area  (1500 km 2 ). The injection site is chosen to be near wells with available core and log data. The base properties used for homogeneous model are the same as for the conceptual model (base case).

Dependency of Plume Size and Pressure
Field on the Number of Injectors Figure 5 presents the plume extension and pressure distribution after 50 years of injection using base case properties.
For the case with one well, the plume radius in the top layer is about 4.6 km, which is consistent with the conceptual model as well as analytical solution radius [7]. It is noticeable that the size of "pressure plume" is much larger (about 65 km), even for one well. In the cases of n number of injectors, one can see the existence of n individual plumes for CO 2 saturation. As the number of wells increases, the individual injector flow rate decreases (fracture pressure constraint) and consequently the plume radius decreases. For pressure, however, strong interference between injectors occurs leading to pressurization of the total area covered by injectors and far beyond that.
It is very important to mention that the dynamics of the pressure field is very different for one well injection and for multiple injection wells, as depicted in Figure 6. One can see that the pressure fall-off after ending CO 2 injection is much quicker for a single well than multiple wells. The pressure fall-off for multiple wells is much slower because of the larger area of pressure influence.

Sensitivity study of injection capacity
Starting with one well, the maximum achievable rate was determined to be as high as 1.1 Mt CO 2 /year (matches the results of conceptual modeling), which is equivalent to 0.055 Gt after 50 years. This flow rate causes the bottom-hole pressure to reach 40 MPa (fracture pressure) at the end of the injection period. This value for the formation fracture pressure was assumed by WASP team at the beginning of the project based on some literature data for Alberta reservoirs. Later, midway through the project, the Geomechanical Simulation Group (part of the WASP team) estimated this value to be around 35-37 MPa. Because our original assumption was very close to this estimate, the 40 MPa estimate was kept in our reservoir numerical modeling (a sensitivity of capacity to different fracture pressures within the range of 30-40 MPa is presented later in this section). When the next five wells are placed in the zone, the corresponding flow rate for each well reduces to 0.625 Mt/year per well with cumulative injection of 0.15 Gt. Increasing the number of wells to ten, brings the flow rate to 0.418 Mt/year per well with total injected CO 2 of 0.209 Gt. Finally, the values for 20 wells are equal to 0.238 Mt/year and 0.238 Gt, respectively. These results are shown in Figure 7 (solid-line curve). In addition, we estimated the reservoir properties required to achieve the target storage of 1.0 Gt. The dashed-line curve on Figure 7 presents the injection capacity of the focus area with the following aquifer properties: porosity of 20%, horizontal permeability of 90 mD. Although these values cause significant difference in the outcome, the limitation in injectivity improvement for more than 10 wells still exists.

Some options to increase capacity (horizontal injection, fracturing)
In the storage process, the term "capacity" could have two meanings. The apparent capacity is the available and accessible pore volume of the aquifer. However, the injection capacity is the amount of CO 2 that can be injected realistically into the formation and is a function of the number of wells and the fracture pressure of the formation. As discussed earlier, for a restricted injection area, such as one in the Nisku study, increasing the number of wells beyond a certain limit (which is controlled by formation properties and injection site area) has a minor effect on the injection capacity. A few methods that could increase injection capacity are discussed below.
The first method is using horizontal wells instead of vertical wells. For vertical wells, it is preferable to use fully penetrated wells over the whole thickness of the aquifer. To find the minimum length of the horizontal well that would increase the injection capacity over vertical wells, the effective radius of pressure disturbance around the vertical injection well, which is again a function of formation properties, should be determined. After the start of injection through a vertical well, the pressure around the wellbore increases rapidly and causes the development of a narrow width pressure peak in the vicinity of the well (Fig. 8).
Using horizontal wells with a total length greater than the scale of the vertical injector's pressure peak "L 1 " (Fig. 8) will diminish these high-pressure peaks in the vicinity of the wellbore and cause the injectivity to increase. For the Nisku formation, this minimum required well length was estimated to be 3000 m. The middle points of the horizontal laterals are placed at the same locations of the vertical injectors. They are all in N-S direction and parallel to each other and located in the middle of the model thickness. In Appendix A2, some details of the gridding for simulation of horizontal wells are presented. Figure 9 compares the effect of different well configurations on the achievable CO 2 injection capacity. While the vertical injectors provide the minimum capacity, their replacement with the horizontal wells will increase storage capacity by~50%. Application of stimulation techniques such as hydraulic fracturing (HF) can also improve injectivity, with the results depicted as the last two bars in Figure 9. Nevertheless, the technical feasibility of implementing these techniques requires careful geomechanical characterization of the formation. For example, maintenance of the caprock integrity is a major concern in using such techniques for secure underground storage of CO 2 . The details of simulation for HF can be found somewhere else [8].

Heterogeneity study
Two different heterogeneity models suggested and provided by geo-statistics group were considered in this study: stochastic and object based modeling. The stochastic modeling was based on existing quantitative data (wireline log, acoustic impedance) and geo-statistical tools. It relies on resistivity-derived porosities and permeabilities (nearly 60 wells). For this study, five equiprobable realizations of properties were generated (Fig. 10). Porosity ranges from 1.3% to 28.6% (mean 4.9%) and permeability from 3.1 to 393 mD (mean 22.37 mD). Details of the development of the geo-statistical realizations, generated by a geo-statistics group for the WASP project, are presented in the geomodeling section [9] of WASP report. It should be noted that the presented geomodels only cover the study area (see Fig. 4); however, the pore volume of the boundary grids (far from the well locations) have been modified such that the total pore volume of the complete model is still honored. This is an important consideration for running such a huge model without sacrificing the accuracy of the results.
Object-based models use "objects" (of same size) with higher porosity and permeability zones. The geometry and distribution of the objects in this study were constrained by dimensions of existing modern carbonate analogs, conceptual understanding of the Nisku carbonate in the Wabamun area, wireline log data, and seismic data.
Two kinds of objects were used: (1) dark blue (minor width = 500 m, major/minor ratio = 5, and thickness = 5 m) and (2) light blue (minor width = 300 m, major/minor ratio = 5, and thickness = 2 m). These are oriented along dip and distributed in each zone (upper, middle, and lower, Figure 11 left) of the Nisku. Although one can see some differences on a small scale for saturation and pressure fields of these two models, the injection capacities for all cases are almost identical   The above results represent the storage capacity when the fracture pressure was set to 40 MPa. To investigate the effect of fracture pressure limit on the storage capacity, simulations were run with sequential Gaussian simulation (SGS) model 1 in Figure 13. Sensitivity of capacity to fracture pressure is shown in Figure 14; a 25% reduction in fracture pressure will reduce the aquifer capacity by 40%, which is significant.
For object modeling, the saturation, and pressure fields are shown in Figure 15.
Injection capacity, as in the case of stochastic modeling, is very close to the homogeneous model (Fig. 16).
Based on this limited number of realizations, it is tentatively concluded that the small-scale heterogeneity (compared to plume size) considered in this study does not exert a strong control on pressure and saturation fields or on overall capacity of the injection site. However, it is expected that for layered systems, and for objects comparable to plume size, the effect of heterogeneity will be substantial. This heterogeneity will especially affect selective placement of injectors. Such a study would require more detailed knowledge of the property distribution within the aquifer.

Long-term Fate of CO 2 Injection
In this section, the long-term fate of CO 2 injection is discussed in association with the following phenomena: • Increased aquifer pressure during and after injection • Migration of CO 2 beyond injection area due to dip • Buoyancy of CO 2 over long period of time (leakage risk) • H 2 S issues

Pressure field evolution during and after injection
As discussed in the previous sections, the pressure in the aquifer (within and around injection area) increases rapidly during the injection period and then gradually decreases toward the initial pressure distribution after injection stops, due to the very large volume of the Nisku aquifer. It is important to know how long it will take for a substantial pressure disturbance to dissipate. The simulation results of a 10-well scenario are presented in Figure 17.
One can see that pressure does not reach initial reservoir pressure (P i = 16 MPa) even 650 years  after injection stops, although the difference in DP = P À P initial is small compared to the maximum difference (DP max = P final À P initial = 24 MPa at the end of injection). The graph of DP versus time is presented in Figure 18, which allows estimation of the time scale of pressure decay. From this graph, it can be inferred that DP falls e (~2.7) times at~120 years, when the pressure fall-off curve starts to level off and thereby providing the time scale for the fate of pressure under injection design.

Effect of aquifer dip on plume movement and its size
The effect of aquifer dip was evaluated with simulation runs using the base Nisku properties (Table 1) and a single well injection rate of 1 Mt/year for 50 years. Simulations were run up to 1000 years after injection started, and two cases were considered for comparison: (1) dip angle = 0°a nd (2) dip angle = 0.5°. The second case corresponds to     Figure 19. One can see that at base conditions the effect of dip angle on the plume movement is marginal (Fig. 19A), although when permeability was increased (while all other parameters remained the same) noticeable plume migration along dip was observed (Fig. 19B).
The above results suggest that if the plume reaches regions with higher permeability, it will migrate upwards, and this should be taken into consideration for management of leakage risk.
Estimation of onset of convection for free phase CO 2 The CO 2 injected into a deep aquifer is typically 10-40% less dense than the resident brine. Driven by density contrasts, CO 2 will flow vertically first and then horizontally spreading under the caprock. If there are breaches in the caprock, leakage could occur through these high permeability zones or through artificial penetrations such as abandoned wells. It is very important to know how long the CO 2 remains as a free phase beneath the caprock and how long complete dissolution of CO 2 into the brine phase will take because this determines the time scale over which free phase CO 2 has a chance to leak out. After injection, freephase CO 2 (gas or supercritical fluid) will be partially trapped as residual saturation and the rest will slowly dissolve in the brine [10]. Depending on reservoir properties, different mechanisms may be responsible for dissolution. Here, we estimate the dissolution mechanisms at Nisku conditions and the associated time scales of dissolution. In the short term (Fig. 20A), during and after injection, some amount of CO 2 is residually trapped and the remainder may be dissolved by natural convection (Fig. 20B). The amount of trapped CO 2 depends on rock fluid properties (relative permeability curves), and in the Nisku case, is about 30%. The amount of dissolved CO 2 depends on reservoir pressure and temperature as well as brine salinity, and in the Nisku case, is about 3-5%. The remaining amount of free CO 2 (~65%) will slowly dissolving into the brine by a combination of diffusion or natural convection.
The important parameter to describe convection is the porous medium Rayleigh number, defined as Ra = kghΔq/ (lφD), where k is the permeability, φ is the porosity, g is the acceleration due to gravity, h is the aquifer thickness, l is the viscosity, Δq is the density difference (CO 2 saturated and fresh brine), and finally D is the molecular diffusion coefficient. If Ra > 49 then natural convection occurs. At Nisku conditions, Ra~400, from which we estimated the onset of convection at these conditions (t onset~8 0 years) and time scale of convective dissolution (t dis~3 000 years). Estimations are made based on Hassanzadeh's model [11]. One can see that the time scale for possible leakage is very long and this can be a serious obstacle for implementation of an injection project due to regulatory and monitoring issues. Active reservoir engineering might be needed (for  example, see ref. [12]) to reduce the dissolution time to the scales of hundreds years.

H 2 S issue
While working on the WASP project, we discovered that when CO 2 is injected into a sour saline aquifer, the H 2 S initially dissolved in the brine will be exsolved and released into an expanding CO 2 plume. A comprehensive study of this phenomenon has been made and all detailed results are presented elsewhere [13].

Summary
In this study we performed numerical modeling of injection of large volumes of CO 2 (1 Gt target over 50 years) into the Nisku Formation located in Alberta. Injection modeling was performed within localized injection area of 30 km 9 60 km. The main objectives of the study were to: • Estimate the injection capacity and CO 2 plume movement and pressure distribution during and after injection; • estimate the time scale of pressure relief upon injection, time scale of CO 2 dissolution, effect of dip angle on plume shape, and its migration; and • assess the possible H 2 S evolution in the CO 2 plume over time and space.
The following observations were made in this study: • The capacity of injection is limited not by available pore space but by ability to inject without exceeding the fracture pressure of the formation. The capacity increases with the number of injectors but increasing the number of wells has a limit even considering a large well spacing. Very strong interference between pressure plumes was observed with no substantial benefits beyond 10 wells. Horizontal injection wells and aquifer fracturing may be considered as alternative options to increase the capacity. Sensitivity of capacity to reservoir permeability, rock compressibility, and well placement was investigated.
• For multiple injection scenarios (n wells), CO 2 saturation plumes have no interference; we see n individual plumes of radius 4-5 km for each injector. The pressure field behavior is totally different from the saturation field. There are no individual pressure plumes but a single large (scale of hundred km) pressure envelope.
• The dip in the Nisku formation does not affect the results (i.e., no substantial plume movement) although free-phase CO 2 may migrate along the dip if it reaches a zone with higher permeability above 100 mD. We estimated that at Nisku conditions, the time scale for pressure relief is 120 years and time scale for free phase CO 2 dissolution is~3000 years where the mechanism of dissolution is natural convection.
• H 2 S dissolved in aquifer brine will be released into the CO 2 plume during injection and will reach a high mole fraction at the outer edge of the CO 2 front.

Appendix 2: Details of Grid Size and Associated Refinement in Simulation Models
Sensitivity of the result to the grid size is a common approach to validate a simulation model. In this project, since the models are very large in horizontal directions, it was inevitable to use coarse grid size to create the global model. However, LGR was used to obtain a suitable level of refinement around the wells. Figure A2 displays the variation of grid resolution for the global model of the Nisku formation (with homogenous properties). The model uses smaller grids in "Region 1," where it is enclosing the injection site and becomes coarser in "Region 2" and away from the injectors (see the map scale on the figure). The horizontal dimensions of the parent grids in Region 1 are 3220 m 9 3220 m. In the following the   refinement procedure and the way that different levels of refinement was applied to each well is demonstrated. A square-shape area around each vertical well with three blocks on each side (area~8000 m 9 8000 m) is considered as the area for the first level of refinement (Fig. A3).
At this stage, each parent block is divided into 121 blocks with 11 divisions on each side. The central block in the refined area is then selected and further refined into 121 blocks (called nested refinement) to achieve the level 2 of refinement. Finally, the very central block from the previous refinement task is divided into 25 blocks to complete the refinement at level 3 (Fig. A4). Table A1 shows the dimensions of the blocks at different levels of refinement.
It should be noted that in all the conducted simulations with LGR, the same properties as those of the parent grid block are assigned to all the refined grid blocks. Also, the proper level of refinement is determined by studying the convergence behavior of the solution. For example, in the case of CO 2 sequestration, the amount of free gas at the end of the injection scenario could serve as a worthy criterion to consider. As such, to obtain the required number of divisions at level 3, we plot the average free CO 2 gas saturation in the whole refined area ver-   sus the reciprocal of degree of refinement. At the same time, proper refinement in the vertical direction is also required to capture the buoyancy effect of free CO 2 movement in the upward direction. Figure A5 shows the free gas saturation in the refined area at the end of a scenario in which the CO 2 injection happens continuously at constant injection rate of 1.0 Mt/year and for 50 years.
According to this figure, the parent block of refined area should be divided into five blocks in each x, and y directions with equal length and width of 5.3 m. This degree of refinement is corresponding to the refinement level at which the gas saturations starts to level off. It should be noted that the level of refinement in the hori-zontal direction should be considered combined with grid size gradient in the z direction. Figure A6 shows the final results of the applied grid size in the vertical direction.
Similar steps were performed to obtain a refined model for simulation of horizontal wells. The first level of refinement was applied exactly the same way it was performed for vertical wells. At the second level, the central column of the refined blocks in center block is divided into 121 blocks each. Finally, the central columns of the refined blocks are divvied into five blocks in the perpendicular direction to the horizontal well trajectory (here E-W). Figures A7 and A8 show the refinements details for levels 2 and 3, respectively.