Time-Lapse Integration at FWU: Fluids, Rock Physics, Numerical Model Integration, and Field Data Comparison

Abstract

We present the current status of time-lapse seismic integration at the Farnsworth (FWU) CO₂ WAG (water-alternating-gas) EOR (Enhanced Oil Recovery) project at Ochiltree County, northwest Texas. As a potential carbon sequestration mechanism, CO₂ WAG projects will be subject to some degree of monitoring and verification, either as a regulatory requirement or to qualify for economic incentives. In order to evaluate the viability of time-lapse seismic as a monitoring method the Southwest Partnership (SWP) has conducted time-lapse seismic monitoring at FWU using the 3D Vertical Seismic Profiling (VSP) method. The efficacy of seismic time-lapse depends on a number of key factors, which vary widely from one application to another. Most important among these are the thermophysical properties of the original fluid in place and the displacing fluid, followed by the petrophysical properties of the rock matrix, which together determine the effective elastic properties of the rock fluid system. We present systematic analysis of fluid thermodynamics and resulting thermophysical properties, petrophysics and rock frame elastic properties, and elastic property modeling through fluid substitution using data collected at FWU. These analyses will be framed in realistic scenarios presented by the FWU CO₂ WAG development. The resulting fluid/rock physics models will be applied to output from the calibrated FWU compositional reservoir simulation model to forward model the time-lapse seismic response. Modeled results are compared with field time-lapse seismic measurements and strategies for numerical model feedback/update are discussed. While mechanical effects are neglected in the work presented here, complementary parallel studies are underway in which laboratory measurements are introduced to introduce stress dependence of matrix elastic moduli.

1. Introduction

1.1. Farnsworth Site Background

The Southwest Regional Partnership on Carbon Sequestration (SWP) is one of seven large-scale CO2 sequestration projects sponsored by the U.S. Department of Energy [1]. The primary objective of the SWP effort is to exhibit and evaluate an active commercial-scale carbon capture, utilization, and storage (CCUS) operation, and demonstrate associated effective site characterization, monitoring, verification, accounting, and risk assessment. The SWP field site is located within the Farnsworth Unit CO₂ WAG (water-after-gas) EOR (Enhanced Oil Recovery) project at Ochiltree County, northwest Texas which is undergoing conversion to a CO₂ flood. All CO₂ utilized by the project is anthropogenic, sourced from a fertilizer and an ethanol plant, and this CO₂ would otherwise be vented to the atmosphere (Figure 1). The CO2 WAG field development scheme being applied at FWU is a is a popular form of tertiary hydrocarbon recovery which also holds promise as a large-scale CO₂ utilization and storage (CCUS) mechanism. The CO₂ WAG process [2] involves cyclic alternation between CO₂ and water injection phases for optimal mobilization and sweep of liquid hydrocarbons remaining in pattern drilling developments after primary recovery and waterflood (secondary recovery). As a potential carbon sequestration mechanism, CO₂ WAG projects will be subject to some degree of monitoring and verification, either as a regulatory requirement or to qualify for economic incentives. Time-lapse (or “4D”) seismic provides a robust method for wide-scale fluid monitoring which has been widely applied in petroleum resource development for decades [3]. More recently time-lapse method has been used to monitor evolution of the CO₂ “plume” on most carbon capture and sequestration (CCS) projects [4]. The Southwest Partnership has conducted time-lapse seismic monitoring at Farnsworth using the 3D Vertical Seismic Profiling (VSP) method.

The efficacy of time-lapse seismic depends on a number of key factors which vary widely from one application to another. Most important among these factors are the thermophysical properties of the original fluid in place and the displacing fluid, followed by the petrophysical properties of the rock matrix which together determine the effective elastic properties of the rock fluid system. Here, is where monitoring of CO₂ WAG systems varies greatly from other oilfield and brine aquifer (CCS) CO₂ storage owing to the thermodynamic conditions dictated by the properties of the original and displacing fluids, reservoir temperature, and pressure. Geologic sequestration of supercritical CO₂ into a brine aquifer, which is the typical case for CCS projects, results in a fluid system with effectively binary fluid properties and relatively simple interface between original and displacing fluids by comparison with miscible CO₂ WAG systems. By contrast, CO₂ WAG operations result in a thermodynamically complex fluid system with multiple fluid contacts and a high degree of ambiguity in thermophysical properties. Further, under miscible conditions an additional phase is introduced. We use the data from the SWP Farnsworth West project in an extensive fluid and rock physics modeling study to understand the unique monitoring challenges presented by miscible CO₂ WAG operations. The rock physics models are applied to compositional reservoir simulator output and the resulting elastic predictions compared to time-lapse 3D VSP surveys acquired at the FWU site.

1.2. Literature Review

There are numerous case histories in the literature documenting the application of time-lapse monitoring for petroleum resource management and carbon sequestration. Here, we focus on the distinct fluid systems presented by the various applications. CCS (brine aquifer storage) projects such as Sleipner [5,6,7], Aquistore [8,9,10], Illinois Basin-Decatur Project (IBDP) [11,12,13], and Ketzin [14] are all non-miscible fluid systems which are thermodynamically stable, have relatively favorable fluid mobility ratios, and in which the single displacing (injected) CO₂ volume is monotonically increasing. The result is a continuous “plume” of CO₂ emanating from the injection well, the interface of which with the aquifer brine may be characterized as a “front”. Similarly, most documented successful (non-EOR) hydrocarbon reservoir surveillance applications of time-lapse seismic have the objective of monitoring an encroaching brine aquifer during field depletion, progression of a water injection front in waterflood operations, or evolution of light hydrocarbon gas upon pressure drawdown, all of which present essentially binary fluid discrimination problems which may be characterized (and modeled) using “plume”, “front”, or “cap” concepts. While not being trivial by any means, these fluid property and displacement scenarios have the added advantage of well-established empirical correlations, analytical models, and reasonable conceptual approximations to simplify and/or constrain interpretation. The binary nature of this problem is further evidenced by the evolution of discrete analytical methods using classification and analytical integration approaches [15,16,17].

By comparison there are relatively few documented case histories of time-lapse seismic monitoring on CO₂ EOR projects. Extensive fluid analysis and rock physics modeling performed by Brown at Weyburn [18] show an estimated maximum of 6.3% variation in compressional velocities between endpoint saturations of 100% brine and 100% CO₂ in high porosity (24%) fractured marly dolostone. While this degree of sensitivity is favorable from a seismic detection perspective, the prediction for intermediate fluid mixtures in a 29% porosity rock with oil and 40% CO₂ are on the order of 1%, which is marginal for seismic detection. Quantitative time-lapse seismic data integration at Weyburn was performed using a novel methodology which optimized a penalty function formed of the distance between simulated CO₂ “front” from ensemble models and that interpreted from thresholded time lapse seismic anomaly maps [19]. At Cranfield Gosh [20] applied a pressure dependent effective media model (PDEM) to invert time lapse seismic data for CO₂ saturation. While success is reported in mapping the extent of CO₂ migration, the author acknowledges the need for additional constraints on gas distribution in order to accurately predict CO₂ movement. Alfi and Housenni [21] compared time-lapse seismic interpretations to simulator predictions at Cranfield. As Cranfield is not a WAG operation the CO₂ “plume” was continuous. Limited success in the comparison was attributed to model deficiencies and uncertainties. No attempt was reported to optimize the simulation model through time-lapse seismic constraints. At the Denbury Bell Creek [22] project researchers performed extensive rock physics and forward elastic modeling from reservoir simulator output. The time-lapse seismic interpretation provided insights into reservoir connectivity which was fed back to simulation modelers for consideration in the subsequent reservoir model update. Incorporated discrepancies were fed back to reservoir engineers for use in model updates. While these case histories show high levels of analytical rigor, all rely either explicitly or implicitly on the concept of a CO₂ “plume” or a distinct CO₂ saturation “front”.

Our literature review also reveals the wide variety of geophysical analyses and resulting attributes which have been utilized for time-lapse seismic data interpretation and integration. These range from computationally simple but robust seismic data transforms which can be extracted without external conditioning, to highly sophisticated and computationally expensive specialized attributes which require conditioning to high quality geophysical logs and/or numerical models. At the Hall-Gurney field in Kansas [23] researchers experimented with use of low cost noninversion attributes for monitoring the effectiveness of EOR in thin, shallow carbonates and were able to identify “overall area effected by injected CO₂”. Time or depth shift (sag, displacement) in events underlying a storage zone is a popular and robust attribute which requires no external constraint for use as a qualitative indicator [24,25]. With suitable constraint on reservoir properties and rock physics, time delay may be used to estimate saturation changes in the injection zone. Time delay has been successfully used at Sleipner [6] and Norne [26]. The next level of analytical rigor is seismic post stack (acoustic inversion), followed by pre-stack (elastic) inversion, both of which have been applied at Sleipner [5,7]. Yet, more sophistication is introduced by implementing wave equation constraints [27]. In rocks displaying anisotropy such as the naturally fractured Weyburn dolostone, multi-component seismic attributes may be used to produce semi-quantitative maps of elastic anisotropy for more representative rock physics modeling [19].

2. Materials and Methods

2.1. General Methodology

We use the extensive body of site characterization data and reservoir modeling performed at FWU as the basis for a comprehensive fluid property and rock physics modeling study, followed by comparison with time-lapse VSP surveys. The calibrated FWU compositional simulation model is used to characterize expected fluid distributions in the injection formation at times corresponding to seismic survey times. These results are used to perform systematic fluid equation of state modeling for representative ranges of anticipated reservoir fluid mixtures. These fluid EOS models are then used with a site-specific petrophysical model for systematic fluid substitution modeling over ranges of representative saturated reservoir rock conditions. These systematic modeling studies reveal important characteristics of the rock-fluid system which are critical for the following comparison of reservoir scale simulations with field time-lapse measurements.

The petrophysical and elastic properties used in our compositonal models and subsequent rock physics computations are correlated through a common 3D porosity distribution which has been developed through analysis of log and core data and interpolated using spatial trends extracted from the available 3D seismic data. First, the 3D permeability distribution required for reservoir simulation was created through poro-perm relationships extracted from the logging/core dataset and interpolated through geostatistical integration with the 3D porosity distribution using the same seismic data for spatial trends. Next, the 3D elastic property distributions required for rock physics computations were developed through correlations with porosity and interpolated through geostatistical integration with the same common 3D porosity distribution, using the same seismic data for spatial trends. In this way we feel that we have maintained correlation across petrophysical and elastic properties throughout our process.

2.2. FWU Geological Model

Site characterization efforts at FWU have produced a rich collection rock and fluid samples, geophysical logs, and multiple time-lapse seismic datasets [28,29]. Core from several wells, including characterization wells drilled specifically for this project, underwent comprehensive petrographic analysis, flowthrough, and mechanical testing. Three wells were drilled as characterization wells (or “science” wells) by the partnership. Approximately 250 ft of core was obtained from each of the new characterization wells. Cored intervals include the entire Morrow B reservoir interval, as well as Morrow shale that underlies and overlies the Morrow B, the B1 sandstone interval, and the Thirteen Finger limestone which forms the remainder of the primary seal. The geophysical logging program for the science wells was designed to support the anticipated geophysical, petrophysical, geomechanical, and geochemical studies as well as coupled process modeling.

Although geotechnical data acquisition activities at the FWU site have transitioned from site-characterization activities to monitoring data collection, ongoing improvements have been made to the site characterization through application of improved processing and analytical methods. Figure 2 shows the type well log for the Farnsworth Unit, the unit boundary, and location of SWP characterization wells. Specialized integration of geophysical logs and mechanical core tests have resulted in creation of detailed wellbore Mechanical Earth Models (MEMs). Pre-stack depth imaging of the 3D seismic dataset has resulted in both higher fidelity structural and stratigraphic imaging and improved image gathers for elastic inversion. Ongoing core flooding experiments and petrographic studies have resulted in definition of distinct hydraulic flow units (HFU’s) [29] within the Morrow B formation as well as corresponding porosity dependencies and porosity-permeability relations. The results of these data analyses have been integrated through a combination of geostatistical and machine learning methods to provide enhanced geologic description, hydrodynamic property estimation (Figure 3), and development of a 3D mechanical property model for support of ongoing numerical elastic and mechanical modeling studies. Although numerical simulations were performed in the Morrow B production interval, porosity and permeability have been interpolated in from the Thirteen Finger to base of the Morrow formation to support ongoing related studies. The property interpolation workflow applied to each formation depended on the data available and the formation characteristics. Integration methods included artificial neural network facies identification from well logs and core, spatial variogram analysis, discrete and continuous distributions, and co-simulation with elastic inversion properties. Due to the limited well log data in all formations except the Morrow B, spatial variograms from seismic impedance were used as proxies for well log data variograms in property interpolation. Such use of variogram proxies, and the use of secondary variables in co-simulation, were justified by observed correlations in available well log data.

The method applied within the Morrow B was distinct from other formations due to availability of legacy well logs and research into HFUs conducted by SWP. A “Winland R35” transform was derived from analysis of core porosity and permeability for 51 wells. Eight different sub populations were identified in poro-perm space and used to create R35 cut-offs defining HUF. Poro-perm relationships were derived for each HFU sub population from core data. The R35 transformation was used to compute R35 logs for the 51 wells with data used in core analysis. For poro-perm interpolation, porosity logs were upscaled into the grid and interpolated through Gaussian co-simulation with seismic acoustic impedance

2.3. FWU Numerical Simulation Model

A field-scale Eclipse 300 compositional reservoir flow model has been developed by SWP researchers for assessing the performance history of the CO₂ flood, and optimizing oil production and CO₂ storage at FWU [30,31,32].

2.3.1. Compositional Fluid Model

The compositional fluid model was constructed from laboratory fluid experiments tuned to an equation of state (EOS) [33,34]. The mixing rules of Pedersen [35] were followed to split C7+ fractions into two pseudo-components using the average molecular weight, average specific gravity and the total mole percent. The 3- parameter Peng Robinson equation of state [36] with Peneloux volume correction [37] was used to perform all the calculations with the resulting composition shown in

Table 1. Fluid composition sampled from the FWU.

Components	Molecular	Molecular	Critical	Critical
	fraction	weight	Temperature	Pressure
	%	gm/mol	°F	Psi
CO2	0	44.01	87.89	1069.8
C1	38.49	16.04	−116.59	667.17
C2	3.86	30.07	90.05	708.36
C3	2.46	44.1	205.97	615.83
C4′s	1.95	58.12	453.65	430.62
C5′s	1.79	72.15	301.12	547.81
C6′s	2.83	86.18	380.71	489.79
HC1 (7–38)	33.48	189.95	802.94	326.19
HC2 (38–70)	15.13	545.65	1077.75	235.69

and phase envelope shown in Figure 4. The viscosity was modeled using the Lohrenz-Bray-Clark correlation [38]. After calibrating the fluid model to equation of state, a slim tube simulation experiment was conducted to obtain the minimum miscible pressure (MMP) for FWU. A one-dimensional 200 cell model was used for the experiment with a CO₂ injection volume of 1.2 pore volume. The MMP of 4009 psia realized from the simulation as compared to an MMP value 4200 psia derived from laboratory experiments provided by the operator represents a less than 5% error.

2.3.2. Numerical Simulation Model

The simulation model for the Morrow B formation was calibrated through primary (depletion), secondary (waterflood), and tertiary (CO₂ WAG) recovery periods using a machine learning assisted methodology and available pressures and injection/production rates and pressures [32]. The original geological model uses a mesh-grid with a dimension of 100 ft by 100 ft. Top and base of the Morrow B formation were modeled as no-flow boundaries. Running simulation using such a fine grid system takes relatively expensive computational cost. Our preliminary investigations indicate one run using the fine-meshed model could take hours, which results in difficulties to the history matching process considering the demands of running hundreds of simulation cases to find the history matching solution. Therefore, the first step of the history matching work is to upscale the model using a coarsened grid of 200 ft by 200 ft mesh without sacrificing the accuracy. Rasmussen et al. [39] presented relative permeability curves based on laboratory experiment which corresponds to each hydraulic flow units from the FWU which is used in the numerical modeling. During the history matching process, the uncertain parameters considered included lateral permeability, vertical permeability anisotropy, and relative permeability curve inputs. A total of 100 simulations were run by randomly combine the uncertainty parameters and neural network based proxies were developed to improve robustness of history matching workflow and save computational time. Particle swarm optimization is employed and coupled with the expert proxy models to minimize the history matching error considering the oil and water rate agreements between the simulated and field observation data. Sun et al. [40] has presented the detailed history matching model and results utilized for this study. The various periods timesteps that coincided with the VSP time-lapse acquisition data were extracted for the rock physics analysis.

2.4. Rock Physics

2.4.1. Properties of Hycrocarbon-CO₂-Brine Mixtures—Analytical Study

Thermophysical properties hydrocarbon-CO₂ fluid mixtures for systematic rock physics sensitivity investigations were computed using NIST REFPROP [41] and SUPERTRAPP [42] databases and FORTRAN subroutines, integrated with python scripts. We use REFPROP for calculating the bulk modulus and density of CO₂ as a function of pressure (P) and temperature (T) because of the established REFPROP’s accuracy for PVT modeling of CO₂ [43]. SUPERTRAPP is used for calculating the bulk moduli and bulk density of hydrocarbon mixture with and without dissolved CO₂. Fluid phase velocities computed with REFPROP and SUPERTRAPP were used to calculate bulk moduli for each phase. CO₂ was in supercritical state at all reservoir temperature and pressure ranges investigated. Depending on modeling pressure and temperature the fluid mixture has as many as 4 phases; original oil, water, CO₂-oil miscible mixture, and free (supercritical) phase CO₂.

The calibrated reservoir simulation model was used to establish representative hydrocarbon-CO₂-Brine mixtures for investigation by extracting fluid composition along injector-producer profiles from WAG simulations. It was determined that, as a result of many decades of primary depletion and waterflood leaving the remaining oil is “dead”, no significant quantities of hydrocarbon gas evolve during WAG production. The native hydrocarbon fluid reaches saturation with CO₂ and additional CO₂ exists as a separate (supercritical) phase. Further, hydrocarbon components do not exhibit selective stripping as a result of the miscible extraction process as has been reported by some experimentalists [44]. Based on this observation we constructed our hydrocarbon-CO₂-brine mixture EOS computations such that as the feed mole fraction of CO₂ and brine in the mixture increase, the feed mole fractions of hydrocarbon components are reduced in proportion with the original composition.

Given the original oil composition (without CO₂);

$${\mathsf{\Sigma }}_1^{ncomp}xH{C}_i=1$$

(1)

CO₂ is introduced to the feed incrementally, such that,

$${\mathsf{\Sigma }}_1^{ncomp}x\prime H{C}_i=1-x{CO}_2,$$

(2)

and,

$${x\prime HC}_i=xH{C}_i\ast \left(1-x{CO}_2\right).$$

(3)

where x’ are the adjusted hydrocarbon mole fractions.

For this analysis a python wrapper script was constructed to systematically sample the relevant EOS and fluid substitution parameter space, execute compiled SUPERTRAPP FORTRAN flash subroutines, calculate thermophysical properties of the hydrocarbon-CO₂-water mixture with additional input from the REFPROP database, and perform fluid substitution. The python script generates a database of fluid properties and fluid substitution results for graphical and statistical analysis. The output from python script was used as input to the rock physics coupled with the numerical simulation model through series of user-defined workflows developed within Schlumberger Petrel platform. All EOS results are converted to volume fractions for conversion to liquid and vapor saturations. Oil and miscible CO₂ Oil mixtures are liquids because the reservoir is above MMP. Free phase CO₂ is in supercritical state. There is no hydrocarbon gas because the reservoir pressure is above the bubble point at the current temperature.

2.4.2. Properties of Hydrocarbon-CO₂-Brine Mixtures—Numerical Model Integration

For integration with the numerical compositional simulation model the fluid properties were obtained from Eclipse 300 EOS calculations. Bulk reservoir fluid density was computed using predicted oil, water, and gas saturations and densities (keywords SOIL, SWAT, SGAS, DENO, DENW, DENG). Bulk fluid modulus was computed as the inverse of the Eclipse 300 solution for total fluid compressibility (keyword TOTCOMP).

2.4.3. Fluid Substitution

The elastic properties of saturated rock were calculated using Gassmann’s relation [45].

$$K_{sat}=K_{dry}+\frac{{\left(1-\frac{K_{dry}}{K_s}\right)}^2}{\frac{\varnothing }{K_f}+\frac{1-\varnothing }{K_s}-\frac{K_{dry}}{K_s^2}}$$

(4)

where;

Φ = Porosity;

Ksat = Bulk modulus of the saturated rock;

Kdry = Bulk modulus of the dry rock;

Ks = Bulk modulus of the mineral constituents;

Kf = Bulk modulus of the pore fluid.

Gassmann’s equation provides a fundamental relationship to relate a fluid saturated formation with an idealized dry frame formation. However, to apply Gassmann’s equation to solve fluid substitution problems the geophysicist must supply dry rock properties. Biot [46] identified systematic relationship between solid matrix moduli (K_s, G_s) and moduli of the material with porosity.

$$K_{dry}=K_s\left(1-\alpha \right).$$

(5)

Introducing the “Biot” coefficient

$$\alpha =1-\frac{K_{dry}}{K_s},$$

(6)

Gassmann’s relation may be rewritten as

$$K_{sat}=K_{dry}+\frac{{\alpha }^2}{\frac{\varnothing }{K_f}+\frac{\left(\alpha -\varnothing \right)}{K_s}}$$

(7)

Applying Wood’s mixing law [47];

$$\frac{1}{K_{fl}}={{\displaystyle \sum }}_{i=1}^n\frac{S_i}{K_i}$$

(8)

where i = fluid phase.

Matrix elastic properties were determined through analysis of the available geophysical logging suite which included magnetic resonance, dipole sonic, spectral gamma ray, pulsed neutron, array induction, formation image, and spontaneous potential. Petrophysical analysis yielded formation intrinsic properties enabling evaluation of elastic moduli for the Morrow B formation. Figure 5 shows the calibration of formation properties to available core data. Figure 6 shows the analysis of shear and bulk moduli used to determine K_dry, G, and α for application of Gassmann’s Equation (3) and elastic modeling (10). The K_dry, G, and K_s obtained from the geophysical log analysis, served as input parameters for the Gassmann fluid substitution model (4). The saturated bulk modulus, expressed as a function of the dry rock frame and pore fluid properties, was computed using the Gassmann Equation (4). A linear regression model was fitted to the data to identify the relationship between the changes in the elastic properties to the rock physical properties. The solid modulus, K_s, was determined by extrapolating porosity to zero on the K_dry vs. porosity plot.

Shear and compressional seismic velocities were computed as;

$$V_s=\sqrt{\frac{\mu }{{\rho }_b}}$$

(9)

$$V_p=\sqrt{\frac{\mu +\frac{4}{3}{K}_{sat}}{{\rho }_b}}$$

(10)

where;

$${\rho }_b=\varphi \ast {\rho }_{fl}+\left(1-\varphi \right)\ast {\rho }_{matrix}$$

(11)

$${\rho }_{fl}={{\displaystyle \sum }}_{i=1}^n\left(S_i\ast {\rho }_i\right)$$

(12)

For analytical studies, these relationships were implemented by python script. For numerical model integration, fluid substitution was performed on the 3D reservoir model and simulator output grids using the PETREL process manager.

2.5. Time-Lapse Seismic Surveys

The project has acquired multiple seismic data sets for site characterization and monitoring. In addition to a full-field surface 3D dataset acquired for site characterization purposes, 3D-VSP surveys were acquired in all three characterization wells [47]. A total of four 3D VSP surveys (pre-CO₂ baseline and 3 monitor surveys), were acquired in well 13-10A, which is the injector in the 5-spot pattern used for this study. The baseline survey was acquired in February 2014 during the well shut-in for conversion to WAG injection. Monitor surveys were acquired January 2015, December 2016, and December 2017. All monitor surveys were acquired after the CO₂ leg of the WAG cycle. Baseline and monitor surveys were processed through an identical three-component (3C) processing workflow. Pre-processing included source/receiver geometry quality control, receiver selection, 3 component orientation, noise attenuation, Surface Consistent Amplitude Compensation (SCAC), 3C wavefield separation, deterministic trace-by-trace wave-shaping deconvolution, and static correction. Images were created from upgoing wavefields using a Generalized Radon Transform (GRT) imaging algorithm. Figure 7 shows maps of depth shift below the injection interval computed for the three p-wave baseline-monitor survey pairs using the method proposed by Nickel and Sonneland [48]. In consideration of acquisition geometry, and image point coverage, images were cropped outside 1000 ft radial distance from the survey well (13–10A). It should be noted here that this is an extremely challenging seismic detection problem due to the relative thinness of the reservoir (~40–45 ft), the estimated wavelength of p waves (~150 ft) in the data at the reservoir interval, and what we will later see to be very subtle fluid effects.

2.6. Numerical Simulaton Model Integration

In addition to aiding in selection of parameter ranges for the analytical fluid EOS and rock physics study, the numerical simulation model was used to model the time-lapse change on elastic properties for comparison with the recorded time-lapse VSP data. Reservoir simulations were configured to output liquid and vapor phase mole fractions for each fluid component, oil, gas, and water densities, and total fluid compressibility as computed by the E300 compositional fluid model. In order to minimize computational effort a sector model was constructed encompassing the 13-10A injector study pattern. Figure 8 shows the full model and sector model domains, the 3D VSP image area, and the well 13-10A/13-12 profile which is used for diagnostic purposes. The reservoir property model was supplemented with elastic moduli based on log correlations as shown in Figure 8. Fluid substitution and elastic modeling was performed on the four-layer Morrow B formation. A shift attribute map was computed for each baseline-monitor pair by vertically integration of the fluid substitution effects within the reservoir. The statistical and spatial characteristics of modeled and measured shifts were compared on absolute value and normalized bases.

3. Results

3.1. WAG Operational Factors

Figure 9 shows the CO₂ and Water injection rates and cumulative volumes injected by well 13-10A during the time-lapse monitoring period. We have adopted the labeling convention that WAG cycles (labeled 1–8) commence with the water leg of the cycle. The baseline survey was acquired just prior to commencement of the cycle 1 water leg, the monitor 1 survey was acquired at the end of cycle 1, just before the cycle 2 water leg, Monitor 2 was acquired at the end of cycle 5, just prior to the cycle 6 water leg, and Monitor 3 was acquired at the end of cycle 8. It can be seen that the durations and injection rates are variable, often due to operational or economic factors. Note that in Figure 9 volumes are reported in STB for liquid (water) and mscf (thousand standard cubic ft) for gas (CO₂). While these are standard conventions for reservoir engineering and production management purposes these units to not adequately represent the relative volumetric proportions within the reservoir which may impact seismic imaging. To normalize the volumes to reservoir conditions we use the gas formation volume factor (FVF) B_g.

$$B_g=\frac{V_{g,r}}{V_{g,sc}}$$

(13)

where;

V_g,r = Volume of gas at reservoir conditions

V_g,sc = Volume of gas at standard conditions

At reservoir conditions of 168 °F and 4500 PSI, B_CO2 equals 0.001455 while B_H2O~1. Volumes will also be adjusted to bulk reservoir volume using the relationship;

$$bulk reservoir volume=\frac{fluid Volume}{\varphi \ast \left(1-S_{wirr}\right)}$$

(14)

where Swirr = Irreducible water saturation.

Table 2. Cumulative volumes of water and gas injected from the time of the baseline survey at surface and reservoir conditions, and in equivalent bulk reservoir volumes using a porosity of 0.15 and S_wirr = 0.3.

		Standard Conditions		Reservoir Volume
Monitor	Date	Water (stb)	CO2 (mscf)	Water (cf)	CO2 (cf)	W/G Ratio
1	1/17/2015	16,550	553,100	885,100	7,665,000	0.12
2	12/3/2016	99,100	1,314,000	5,298,000	18,020,000	0.29
3	1/1/2018	141,000	1,623,000	789,100	22,490,000	0.33

lists cumulative volumes of water and gas injected from the time of the baseline survey at surface and reservoir conditions, and in equivalent bulk reservoir volumes using a porosity of 0.15 and S_wirr = 0.3. After the short initial WAG cycle the Water: CO₂ bulk reservoir volume injection ratio is approximately 1:3 (0.29–0.33). There are several intervening WAG cycles within the M1–M2, and M2–M3 survey times.

In order to achieve a perspective of fluid volume with respect to time-lapse survey timing we invoke a simplistic conceptual model of non-mixing (piston displacement), annular rings representing the sequential fluid phase injection cycles to compute an “equivalent cylindrical radius” for each hypothetical phase front. Figure 10 shows the hypothetical equivalent annular fluid fronts using a reservoir thickness of 45 ft at each time lapse monitor survey time.

While the non-mixing piston front assumption allows us to gain spatial and relative volumetric perspective with regard to injected volumes and pattern pore volume, a more accurate representation of the thermodynamics and hydrodynamics of the system are reflected in numerical simulations. Figure 11 shows streamline (top) and compositional (bottom) simulation results for the monitor survey times. Streamlines clearly show preferential drainage in the pore space between the injector and producing wells. Compositional simulations of gas saturation are consistent with streamlines showing preferential flow of CO₂ toward injectors, albeit exhibiting asymmetry which suggests either differences in well control parameters and/or performance characteristics, or heterogeneity in reservoir properties. Inspection of production data (Figure 12) shows breakthrough in all four producers shortly after acquisition of M2. Figure 12 shows the simulated mole fraction of CO₂ (top) and oil saturation (bottom) along the study profile between wells 13-10A and 13-12. Simulations show no CO₂ breakthrough at the time of M1. Interestingly, oil saturation increases in front of the CO₂ front at the time of the M1 survey (Figure 13). It is hypothesized that this is the result of banking of oil by the leading cycle 1 waterfront before cycle 1 CO₂ moves in and the miscibility process develops.

As a final step in our review of operational factors we inspect the variation of predicted hydrocarbon compositions in time and space as the oil is mobilized and extracted by the miscible process. We normalized hydrocarbon composition by component mole fractions along the study profile (Figure 14) at each survey time and observe that the composition of the remaining oil stays nearly constant, suggesting that there is no significant selective stripping of components by the miscible process.

3.2. Analytical Model Rock Physics Investigations

Based on our data and numerical model review we designed a systematic analytical fluid EOS and rock physics study to investigate the expected elastic response of the reservoir within representative ranges of fluid composition, rock properties, and reservoir thermodynamic state. Parameter ranges for the study are listed in

Table 3. Fluid EOS and rock physics study parameter ranges.

Parameter	Min	Max
Porosity (fraction)	0.075	0.175
Water Saturation (fraction)	0.3	0.75
Oil Saturation (fraction)	0.27	0.7
CO2 Feed Fraction	0	1
Hydrocarbon Fractions	Proportional
Temperature (°F)	163	173
Pressure (psi)	4000	6000

.

Used together, our petrophysical data analysis, the python implementation of NIST subroutines, and databases facilitated a comprehensive investigation of a broad range of realistic fluid substitution scenarios providing valuable insights into the time-lapse integration problem. Here, we show a selection of graphical results which capture many of the important behaviors of the rock-fluid system resulting from WAG implementation at FWU.

Figure 15 shows contour plots of reservoir fluid compositions which represent stiffness endpoints. The top row with the maximum possible water saturation and no CO₂ represents the stiffest fluid in the system. The bottom row with the minimum possible water saturation and maximum CO₂ represents the softest fluid in the system. Contours are relative to the attribute value at user selected reference temperature which is annotated on each plot. For this comparison we used nominal reservoir temperature and pressure (168 °F, 4500 psi) as the reference point. These plots show that the saturated rock properties have little sensitivity to temperature for both fluids. Further, the stiff fluid shows negligible variation with pressure from a seismic detection perspective, but can vary as much as 20% from nominal conditions for highest anticipated pressures. However, this interpretation must be put in the perspective that such high pressures are not optimal for the production scheme and should be considered anomalous.

Figure 16 illustrates another informative visualization tool, again comparing (endpoint) stiff versus soft fluids at a selected nominal reference temperature and pressure. The four plots in Figure 16 show the variation of normalized elastic properties for hydrocarbon-CO₂-water mixtures (left column), saturated reservoir rock (right column), stiff fluid (top row), soft fluid (bottom row) as a function of CO₂ molar (feed) fraction. All curves are normalized to their value at xCO₂ = 0. Additionally, shown are the fluid mixture vapor (red) and liquid (blue) fraction curves which clearly show the saturation point of 0.75 mole fraction of CO₂. Inspection of the top and bottom plots on the left shows the approximately 7% drop in velocity of the reservoir fluid at zero mole fraction CO₂ (owing to the difference in water saturation), and a much-amplified response to the addition of CO₂ for the softer fluid, again owing to replacement of water with potentially miscible hydrocarbon. Comparing left columns to right illustrated the extreme dampening effect of the rock frame on the elastic response. The difference in velocity between stiff and soft fluids (maxim vs. minimum water saturation) at zero CO₂ mole fraction is only ~0.5%. Velocity of the rock saturated with the stiff fluid drops by only ~0.6% at CO₂ fractions up to saturation before the existence of free phase CO₂ begins to take effect on the bulk fluid properties. For the soft fluid the response to additional CO₂ is somewhat stronger but still at just under 1% at CO₂ saturation.

Finally, we looked at elastic property variations in saturated rock at nominal temperature and temperature (168 °F, 4500 psi). Figure 17 illustrates contoured variations of p wave velocity with CO₂ mole fraction and porosity (left), and water saturation (right). The left plot shows a high degree of sensitivity with respect to porosity, which is a source of significant statistical uncertainty in 3D models, while showing little sensitivity to the quantity of CO₂ except for the inflection at saturation (xCO₂ = 0.75). Variation with water saturation and CO₂ mole fraction (right) shows a complex response surface with most significant variations related to the phase transition and at lower water saturations.

3.3. Model Rock Physics from Compositional Simulator

Elastic properties were computed from 3 dimensional arrays of reservoir porosity and reservoir fluid properties at each time-lapse survey time. Computed velocities from the baseline survey and each monitor survey were used to compute synthetic time-lapse depth shifts. Figure 18 shows the computed time-shifts (color maps) and the contour for mole fraction of CO₂ at the saturation level of 0.75. Figure 19 shows distributions of computed monitor survey depth shifts. Monitor 1 survey has the fewest depth shifted samples (simulator blocks) with depth shifts. The number of shifted cells increases over time as expected injection of additional CO₂. Depth shift magnitude for M1 is distributed uniformly by comparison to those for M2 and M3 which show apparent shifting of magnitudes downward. It is possible that this is due to a greater amount of free phase CO₂ in the early stages of miscibility development, with more complete mixing over time yielding lower magnitude depth shifts. In all cases the computed depth shifts are extremely small from a seismic detection perspective.

3.4. Comparison with Time-Lapse Measurements

We compare simulated depth shifts with Schlumberger displacement attribute computations for the three monitor survey times. Figure 20 shows maps of simulated depth shifts (top) and Schlumberger displacement attribute (bottom) in the zone beneath the Morrow B at M1 (left), M2 (center), and M3 (right) survey times. Color scales have been adjusted to normalize the visual comparison between simulated and measured values. Although the visual comparison is not encouraging, we note that the spatial trends in the measured displacement maps are suggestive of preferential fluid saturation changes toward producing wells, consistent with the known hydrodynamics of the 5-spot pattern. Differences from the simulated distribution may be due to local heterogeneities not captured in the reservoir model and not corrected by the calibration which was achieved through optimization at a global parameter scale.

Figure 21 compares the frequency distributions of absolute and normalized, measured and simulated time-lapse shifts for the three monitor surveys. Simulated and measured datasets were normalized independently to their maximum values. Maximum shifts were 0.68 ft for simulations and 3.8 ft for the measured displacement attribute for a ratio of 5.6. Ratios of average distribution values were 7.5, 10.4, and 3.7, respectively for M1, M2, and M3.

4. Discussion and Conclusions

We are tasked with interpretation and integration of time-lapse VSP data optimizing the FWU reservoir simulation model. From previous experience and literature review we know that WAG operations result in fluid properties and seismic detection scenarios that are much more challenging than immiscible systems. Because these processes are difficult to isolate for investigation in a complex compositional reservoir simulation, we first performed an analytical fluid study which enabled us to isolate and investigate complex thermodynamic processes. We reviewed the field operations at FWU in terms of the fluids in place, and the volumes and timing of injected fluids. We discovered that, contrary to the concept of a discrete CO₂ plume, as a conformance management technique the WAG operation is specifically designed to create an optimal CO₂-hydrocarbon mixture with water as an intermediate phase. We used available data and models to perform a systematic EOS analysis of WAG fluid thermodynamics. We found that the FWU WAG operational scenario results in a thermodynamically complex transient distribution of fluid mixtures with largely similar thermophysical properties while the hydrocarbon undersaturated with respect to CO₂. We then used Gassmann’s model to perform a systematic rock physics study which showed that already subtle differences in fluid properties are heavily muted by the rock framework. From these analytical studies we conclude that the most likely seismic detection scenario at FWU is for the existence of free phase CO₂ which exists where the available oil is saturated, or the miscibility process has not yet fully developed.

We applied our rock physics model to compositional simulations using the production data calibrated FWU reservoir model. These simulations verify a continuum of saturations rather than a distinct distribution of high concentration CO₂ which might be characterized as a plume, or as forming a discrete fluid front. These gradational compositional variations are reflected in similarly gradational spatial distributions of simulated depth shifts. Poor spatial correlations between simulated and measured depth shifts may be attributed to either data resolution or model error. Clearly, strict application of seismic imaging fundamentals tells us that these shifts are all well below the seismic resolution. However, it is not yet clear whether or not we are able to achieve the less stringent goal of merely detecting changes in bulk reservoir properties with the given data. As regards to seismic data fidelity, interpretating estimates of effective media properties at seismic wavelengths is a well-studied problem. Worthington [49] reports ratios of up to 10:1 for estimates of bulk reservoir properties from multi-scale measurements of fractured media compliance using cross-well seismic and laboratory ultra-sonic measurements. On the side of the model accuracy, the model used was calibrated at the global (full field) scale. As such, the calibration at the pattern scale may be non-unique with respect to production data only. Recognition of such ambiguities and non-uniqueness in geostatistical model property distributions is one of the main motivators for introduction of time-lapse seismic data as a calibration constraint. If the fidelity of our measured data can be verified, then the spatial discrepancies between simulated and measured shifts can be used to drive model updates. This subject is currently under study.

Other potential sources of simulated versus measured discrepancies owing to numerical simulations but not investigated here include: the accuracy of the Eclipse 300 assumption of first contact miscibility, the effects of simulation grid scale on the modeling of miscible process development, and the correspondence between the Eclipse calculation of total fluid system compressibility and the assumptions used by the Wood’s law mixing equation and Gassmann’s equation. We also assume that the time-lapse response of the reservoir is independent of pore pressure and effective stress. While we are aware of ultra-sonic measurements under stress on Morrow B reservoir samples which indicate potentially significant velocity versus stress sensitivity, inclusion of such effects was not within the scope of funding for this project.

5. Avenues for Future Work

It is clear that successful integration of the time-lapse seismic data at FWU depends on establishing the confidence level for the seismic data. This is not a black and white issue and should not be considered a binary decision. All quantitative integrations provide a mechanism for scaling the contribution of each observation in computation of model updates. Unfortunately, a rigorous estimate of confidence in seismic data is difficult to establish. It is possible that methods based on data signal to noise characteristics may be implemented. However, these are based on seismic amplitude data so a strategy for application to displacement data would be required.

Quantitative integration of the time-lapse seismic data for reservoir model updating presents a number of spatial and temporal sampling challenges. One of these is the discrepancies between spatio-temporal sampling characteristics of the different observation domains. Seismic data are finely and spatially sampled (~50′ × 50′) but coarsely sampled temporally (1–2 years). Production data are coarsely sampled spatially (~1000′ × 1000′) but finely sampled temporally (hourly, daily). Further, the dynamic processes dictating fluid compositional changes and movement in 3D space over time are very complex. This is why that it is essential to use a transient compositional simulation model as the basis for integration. Although integration of simulation models for 3D interpretation is not new, recent advances in distributed computing and machine learning techniques (such as proxy model based optimization) are fuelling progress. Another challenge is the superposition of effects from sub-pattern scale fluvial system heterogeneity and the thermodynamically active fluid system effects. First, we feel that further lateral refinement of the simulation grid scale is needed. An update and optimization scheme such as the one implemented by Ampomah et al. [32] or Sun et al. [40] for production data calibration of the reservoir model and optimization of development schemes at FWU may be used for time-lapse data integration. In these studies, machine-learning technologies such as response surface models (RSM), multi-layer neural networks (MLNN), and support vector machines (SVM) are used to develop proxy models for numerical simulation results against production data. These proxy models are coupled with suitable optimization strategy (such as the evolutionary strategy) to achieve optimization of a penalty function which includes historical data. For time-lapse integration, the proxy model would include both production data from pattern wells and time-lapse seismic measurements. The required proxy for EOS and rock physics response could be developed through integration of NIST databases in a similar manner as was used in our analytical fluid and rock physics study. Although the full field reservoir model property distribution was developed through a pixel based geostatistical method, it may be necessary to implement an object-based property population method using a fluvial system model in order to impose and systematically update representative anisotropy in model porosity and permeability updates. These methods are currently being applied for calibration of coupled hydro-chemical simulation models by SWP researchers [unpublished]. The proxy model process for time-lapse integration may be preconditioned with a geologically realistic and equiprobable realization ensemble strategy such as implemented by Souza et al. [17].