The heterogeneities and depositional trends in point bars have stimulated several scholarly contributions among researchers and modelers (e.g., Allen, 1964, Allen, 1965, Allen, 1970; Visher, 1964; Thomaset al., 1987). The heterogeneities are formed by episodic migration of sinuous channels, which leads to erosion and deposition of channel sediments. Several periodic lateral accretions are prominent patterns that are observed in point bars (Ghazi & Mountney, 2009; Miall, 1988). The lateral accretions define the aerial dimensions of the heterogeneities while the inclined heterolithic stratifications (IHS) define the vertical dimensions (see Figure 1A). Within point bars, a fining upward trend is common. Sand-prone sediments dominate the bottom of the sequence and transition to silt-prone and mud-prone sediments at the top (Thomas et al., 1987; Labrecque et al., 2011; Fustic et al., 2012). In the past, research studies (e.g., Li and White, 2003; Yue and Shiyue, 2016; Sun et al., 2017) have shown that the lateral accretions (and IHS) constitute the most important heterogeneities that influence fluid flow in point bar systems. This is because of the shale drapes that occur along the surfaces of the lateral accretions and the IHS. The shale drapes are potential flow baffles (Richardson et al., 1978; Hartkamp-Bakker and Donselaar, 1993); they compartmentalize the point bar reservoir and greatly reduce CO₂ storage capacity of the point bar (Issautier et al., 2013; Issautier et al., 2014). In this study, we will make an attempt to represent these heterogeneities in the geological model of the point bar.

The geological modeling process basically involves recreating the channel flow path and its subsequent migration. The heterogeneities are gridded separately and combined to form a 3D point bar grid. Finally, the point bar property distribution is modeled within the gridded geometry using geostatistical simulation, by incorporating a grid transformation scheme. Figure 1B summarizes the workflow for modeling the point bar reservoir.

Grid generation constitutes a crucial part of the workflow for preserving the point bar heterogeneities, and internal architecture. The scheme implemented in this study generates grids that adequately capture the curvilinear geometry of point bars. The gridding scheme was implemented separately for the lateral accretions and IHS. A domain of interest is initially defined, which for the lateral accretions is the region bounded by the initial and current channel path (see Figure 5A). If the number of grid blocks along each channel path is n x , then the cumulative distance at a grid node i along each channel, denoted l i , can be computed as l i   =   L n x   ⋅ ( i − 1 ) . l i is substituted into Eqs 2a, 2b, to generate the coordinates of the grid nodes ( X i j , Y i j ) along the channel paths (see Figure 5B). In Figure 5B, all the grid nodes on any of the gridded paths have their corresponding pairs on the other gridded path. Specifying the number of grid blocks between these corresponding pairs of grid nodes as n y , we can compute the coordinates of the grid nodes between the gridded channel paths using Eqs 3a, 3b, to complete the gridding of the lateral accretions (see Figure 5C), with the equivalent curvilinear grid displayed in Figure 5D.

The same procedure is repeated to generate the grid for the IHS (see Figures 5E,F).

Combining the grids for the IHS and the lateral accretions generates a 3D grid for the entire point bar. x i j = X i j ± ( j − 1 ) ⋅ D i j n y ⋅ cos ⁡ α i j (3a) y i j = Y i j ± ( j − 1 ) ⋅ D i j n y ⋅ sin ⁡ α i j (3b) where D i j   and     α i j are respectively, the distance and angle between a grid node at position i , j on the current meander path and its pair on initial meander path; n y is number of grid blocks between meanders; X i j and Y i j are grid node coordinates along meanders.

As indicated earlier, traditional geostatistical methods when applied on regular Cartesian gids cannot properly preserve the curvilinear characteristics and the heterogeneity of point bars. To overcome this challenge, a grid transformation scheme was incorporated. Grid transformation is a way of standardizing the position parameters with respect to the channel boundaries, where the layers are unraveled and flattened onto an orthogonal grid. Geostatistical simulation is then conducted in the transformed space, after which the properties are mapped back into the original curvilinear grid. A summary of this workflow is shown in Figure 6. To account for the pinch out of point bars, the pinch-out grid blocks are treated as net-effective zero thickness blocks, using volume modifiers; therefore, the pinch-out arrays will not contribute to fluid flow to affect flow simulation results.

Identifying the direction of spatial continuity is also necessary for preserving the heterogeneities and depositional trends. For point bars, maximum continuity is in the downstream direction while least continuity is in the perpendicular direction between successive IHS sets (i.e., accretion direction). Between these limits is the medium continuity, which occurs in the dip direction. Accordingly, the maximum range ( h m a x ) , medium range ( h m e d ) and minimum range ( h m i n ) were chosen to be in the direction of continuity in the downstream flow, dip, and accretion directions, respectively. Figure 7 shows the directions of spatial continuity as depicted in curvilinear and orthogonal space. The point bar is modeled with 5 inclined layers (i.e., IHS sets) and for each set, h m a x , h m e d and h m i n were specified as the major, medium, and minor axis of anisotropy, respectively, for the variogram model (Table 1). Using the SGSIM algorithm (Deutsch & Journel, 1998; Remy, 2004), the IHS sets were modeled separately on an orthogonal grid with dimensions 150 × 30 × 25 . The complete 3D point bar model was obtained by stacking the IHS sets orderly from set 1 to 5 and mapping their properties into the point bar curvilinear grid of dimensions 150 × 150 × 25 (see Figure 8A). Typical trends seen in point bar systems like the Cranfield reservoir include overall fining up trend in the sediments. High porosity sediments like conglomerates occupy the bottom of the sequence followed by relatively lower porosity sediments like sandstones, muddy-sandstones, and mudstones. These trends were imposed in the simulation with locally varying prior probabilities. It was ensured that the probability of encountering high porosity is maximum at the bottom of each IHs set and minimum at the top. Horizontal slices taken at different depths across the point bar illustrate a fining upward trend as captured in the geostatistical simulation (see Figures 8B–D).

The point bar property model displayed here is one of the 50 realizations that were generated using geostatistical simulation. The reservoir response computed over the suite of realizations would provide an assessment of uncertainty associated with response predictions. For a reliable assessment of the CO₂ sequestration potential of the point bar reservoir, the uncertainty in the point bar geologic model must be reduced by constraining the model to available flow-related data. In the next section, we will discuss how to achieve this by further conditioning the realizations of the point bar model to observed CO₂ injection data.

Model calibration (history matching) is an essential step in characterizing reservoirs and forecasting future reservoir performance. It is grounded on the premise that the static or primary reservoir variables (e.g., porosity, permeability, seismic behavior etc.), influence the dynamic or secondary response of the reservoir (e.g., bottom-hole pressure, CO₂ saturation, CO₂ injectivity etc.). Accordingly, in this study, the static reservoir variables will be systematically adjusted to ensure that simulated dynamic variables acceptably agree with the observed historic data.

In the model calibration procedure, we consider the fact that: 1) reservoir geometry influences fluid flow response variables (e.g., pressure and transmissibility), and therefore should be accounted for in the history matching, 2) the Cranfield reservoir exhibits non-gaussian characteristics (due to the continuity of the depositional structure); and furthermore, the joint relationship between the primary reservoir variable (like permeability), and the secondary variable (like bottom-hole pressure) is non-linear. We use a two-step ensemble-based data assimilation procedure where: step 1 uses ensemble Kalman filter (EnkF) to update ensembles of point bar reservoir model geometries, to select the geometry that yields the closest match to observed data; and step 2 applies indicator-based data assimilation (InDA) to update ensembles of permeability models within the optimal reservoir geometry determined in step1.

With the limited coverage of wells in the CO₂ injection area, there is likely to be significant uncertainty associated with the prediction of reservoir geometry. Calibrating models for point bar geometry using injection data is necessary to constrain the uncertainty. Every angular displacement ( ω ) in the SGF defines a unique point bar geometry; therefore, 50 realizations of ω were drawn from a normal distribution with mean 35 and standard deviation 10, to generate an ensemble of point bar reservoir geometries. Figure 9 shows some of the realizations of point bar geometry ensembles that were generated. The ensemble of geometries are the initial uncertainties that would be updated upon the availability of secondary variables (e.g., bottom-hole pressure, reservoir pressure, injection rate). Using CMG-GEM simulator (CMG-GEM, 2019), flow simulation was run on the ensemble for a period of 450  days. The geometries were described using curvilinear grids of dimension 150 × 150 × 25 , making 562500 grid blocks in total. Reservoir properties like porosity and permeability were kept the same over the ensemble. The injection schedule used for simulation is shown in Figure 10.

Flow simulation was conducted on the ensemble to obtain the corresponding simulated dynamic bottom-hole pressures that can be used to infer the required covariance between state parameters, and subsequently update the initial geometry ensembles. The error variance representing the measurement error is made proportional to the variance of bottom-hole pressure observed over the ensemble. Updates were performed on the reservoir geometries using the EnKF formulation in Eq. 4 (Kumar and Srinivasan, 2019). z o a =   z o f + [ C p z ( h j o ) ] 1 × n [ [ C p ( h l j ) ] n × n + d i a g [ C ∈ ∈ ( h l j ) ] n × n ] − 1 [ P o b s − P s i m ] n × 1 (4)

In Equation 4, C p z ( h j o ) terms are the covariance between the primary and secondary variables, the matrix containing the C p ( h i j ) terms are the covariance between the dynamic response variables. Equation 4 yields the updates z o for the ensemble locations o = 1,2,3... n ; z o f is the initial value of the primary variable, which is angular displacement in our case; The superscript a represents the updated value. P o b s and P s i m are respectively, the observed and simulated secondary variables, which are bottom-hole pressures in our case. ∈ is a vector of uncorrelated observation errors (which is made proportional to the variance of pressures observed over the ensemble), C ∈ ∈ therefore becomes a diagonal matrix. To account for possible errors in the observations, samples (number of samples equals the size of ensemble being used) drawn from the distribution N ( 0 , C ∈ ∈ ) are added to the observed data. (Please note that in this study C ∈ ∈ = 10   k P a ).

The update (or error) term constitutes all the terms after the first term in Eq. 4. Figure 11B illustrates the updates realized after EnKF implementation. It could be observed that the updates at higher angles are higher than those at lower angles. After performing EnKF updates and running flow simulation on the updated ensembles, a reduction in the spread (i.e., uncertainty) in the simulated bottom-hole pressures is observed (see Figures 11A,C). The reservoir geometry with the least uncertainty (i.e., least error) has an approximate angular displacement of 19.4 ° . For this reservoir geometry, the update made to the point bar geometry is illustrated as the black point in Figure 11B. The smaller update indicates that the updates using the CO₂ injection data have stabilized. The corresponding simulated bottom-hole pressure after EnKF updates is the red plot in Figure 11C, and the updated geometry is shown in Figure 11D. This reservoir geometry was selected for the next phase of the history matching process.

Figure 11 also suggests that the geometry alone does not fully explain the dynamic response characteristics of the point bar reservoir under study. A further confirmation is seen after running flow simulation on the ENKF-updated reservoir geometry ensembles. As shown in Figure 11C, after ENKF updates, the simulated bottom-hole pressures still exhibit appreciable uncertainty even though, overall, the uncertainty is reduced. Yet, this calibration step is important as it allows us to select a reservoir geometry that will ensure a more successful history match of static reservoir properties.

Please note that the standardized errors on Figure 11B, denoted S E ∈ [ a , b ] were obtained by using Eq. 5, to scale the updates to a desired range [ a , b ] . S E = E − min ( E ) max ( E ) − min ( E ) ⋅ ( b − a ) + a (5) where E is the absolute error or update and S E is the standardized absolute error for the ensembles.

The applied thresholds and indicator definitions are used in the InDA formulation to perform updates. Updates are performed by updating the cumulative density function ( c d f ) of the primary variable. In the InDA formulation, the update c d f equation, expressed as a conditional expectation of the indicator variables is written as: F ( Z α f ( u ) ≤ z k | Y p )   =   E ( ( I ( Z α f ( u ) , z k ) ) +   ∑ i = 1 i = n λ i ( Y i − E ( Y i ) ) (9)

Like the EnKF update equations discussed previously, λ i are the Kalman gain computed on the basis of the indicator cross-covariance between the primary indicator data I and the secondary indicator data Y. InDA is immune to non-linear transformations because of the binary nature of the indicator variable. In addition, it does not suffer from Gaussian assumptions as applying Eq. 9 directly generates the conditional c d f without making any parametric assumptions.

Using the point bar geometry corresponding to the lowest error, the next task is to perturb the spatial variation of rock properties to match the observed data. Like the traditional geostatistical methods, InDA runs optimally on orthogonal grids; therefore, we modified its implementation by incorporating the grid transformation scheme within the update process.

At the time of modeling, hard data for permeability was not available, therefore, from the porosity ( ∅ ) ensemble of models that have been generated earlier, the permeability ( k ) ensemble of models were generated from Eq. 10, which comes from micro-modeling for enhanced small scale porosity-permeability relationship (Boisvert et al., 2012). ln ( k ) = 4121.2 ⋅ ∅ 4 + 3963.6 ⋅ ∅ 3 − 1353.3 ⋅ ∅ 2 + 202.05 ⋅ ∅ − 4.3571 (10)

The use of Eq. 9 for performing updates is contingent upon defining appropriate indicator thresholds for the primary and secondary variables. Figure 12A represents the c d f of the reference permeability model under consideration, and the thresholds that were retrieved. The secondary data thresholds were defined based on the mismatch between the simulated and observed bottom-hole pressure data from the Cranfield injection well (see Figure 12B)

Using the defined thresholds and the indicator definitions, CO₂ flow simulation was run on the entire initial ensemble of permeability models. Updates were performed by assimilating the mismatch between the observed and simulated bottom-hole pressure at the CO₂ injection well location, using Eq. 9. Figure 12C shows the initial and InDA updated c d f at each location.

To determine the updated permeability whose initial value is say, 800 md, it is read on the horizontal axis (position 1) and the equivalent probability is read on the initial c d f (position 2). This same probability is read on the updated c d f (position 3). The updated permeability value is therefore the permeability that corresponds to the initial probability drawn from the updated   c d f , which is position 4 on Figure 12C.