Abstract
Reservoir development forecasts depend on accurate descriptions of the spatial distribution of rock properties that impact subsurface fluid-flow pathways and volume connectivity. Reservoir models constructed using geostatistical methods combine analogous facies dimension data with sparse subsurface data to predict spatial variations in rock properties. This study uses a physics-based depositional process model to define realistic facies variations within a river-dominated delta deposit formed during multiple shoreline regressions and transgressions. Geostatistical models are conditioned to varying amounts of information extracted from the depositional model to examine how well they reproduce the facies patterns. Reservoir simulation is used to examine the impact of analogous dimension data and varying conditioning constraints on reservoir performance predictions of water displacing oil. The dimensions of surface depositional features underestimate the continuity of preserved facies patterns, proportional grids following major flooding surfaces allow significantly better predictions than uniform rectangular grids, and trend constraints are more important when defined facies correlation length is significantly less than well spacing. When geostatistical model parameters are poorly chosen, reservoir simulation of the resulting weakly-structured facies patterns overpredict recovery and water breakthrough time. It is demonstrated that process-based depositional models can be used to optimize geostatistical model construction methods and input parameters to reduce uncertainty of reservoir development assessments.
Forecasts of reservoir production behaviour require three-dimensional models of geological heterogeneities that influence subsurface fluid-flow pathways. Accurate characterization of reservoir variations is difficult because data to constrain predictions tend to be either local and widely spaced (e.g. core and well logs) or lack details of fine-scale rock property patterns (e.g. seismic images of strata at typical reservoir depths). The details needed for reservoir characterization are generally added by applying geological concepts to correlate larger-scale reservoir units between wells and incorporating analogous information to infer finer-scale heterogeneity patterns appropriate for the interpreted depositional environment of the reservoir internal rock facies. Geostatistical simulation is the most common method to reproduce interwell-scale heterogeneity patterns within correlated reservoir units because these methods are well developed in standard reservoir modelling software, and results can be conditioned to field data and used for uncertainty analysis by defining multiple geostatistical realizations.
Analogy information is used in geostatistical simulation to provide finer-scale spatial correlation structure, typically in the form of semi-variograms for sequential Gaussian simulation (SGS: Journel and Deutsch 1997; Deutsch 2002) or training images for multipoint statistical methods (MPS: Deutsch 1992; Guardiano and Srivastava 1993; Strebelle 2002). Semi-variograms summarizing lithic variation patterns at interwell scale can be calculated from high-resolution subsurface analogous examples or inferred from the size, shape and organization of lithic bodies observed in analogous modern depositional systems or outcrops. Multipoint statistical methods can provide greater control over predicted lithic spatial continuity patterns by sampling a training image (Strebelle 2002; Hu and Chugunova 2008) but developing suitable 3D training image descriptions for the range of possible analogous depositional environments has proven difficult.
Training images can be defined by predicting facies patterns within high-resolution subsurface models of extensively drilled mature fields, closely sampled modern deposits, large three-dimensional outcrop exposures, or using geological object-based or process-based depositional forward models. In most realistic cases, however, rock properties vary over hierarchical nested scales, and non-stationary trends across the deposits complicate geostatistical quantification of the regionalized (stationary) patterns needed in MPS applications. The critical sensitivity of MPS simulations to the training image has led to development of specialized programs to generate stationary training images (Deutsch and Tran 2002; Maharaja 2008; Pyrcz et al. 2009; Boucher et al. 2010; Fadlelmula et al. 2016) and databases summarizing the dimensions of the nested hierarchy of lithic bodies formed in different depositional settings (Reynolds 1999; Gibling 2006; Colombera et al. 2012, 2016; and many others). Most of these applications define geobodies as geometrical shapes with relatively simple spatial organization relationships because it is easier to generate stationary images using these methods than by using real analogous examples. Training images are typically inferred to be reasonable if they look qualitatively like conceptual facies models or surface images of depositional bodies within modern environments (Fig. 1a and b), and the scaling of the reservoir and analogous bodies are broadly similar.
Comparison between modern deltaic surface geomorphic features and variations within the depositional model. (a) Surface features on the modern Mississippi River delta show multiple scales of distributary channel bifurcation and associated depositional units. Major delta lobes and sub-lobes form as the main river channel avulses (deltaic complexes of Fisk and McFarlan 1955 and Frazier 1967). Deltaic elements are defined by a local cluster of mouth-bar storeys at the end of a terminal distributary channel (USGS EarthExplorer, Image Landsat). (b) Details of mouth-bar storeys deposited within the Wax Lake deltaic depositional element (USGS EarthExplorer, Image Landsat, overlain by bathymetry from Shaw and Mohrig 2014). (c) Surface topography of the delta model used in this study shows a similar pattern of terminal channel and lobe deposits as in Wax Lake but these elements are smaller scale than the Mississippi system. (d) Detailed topography of one deltaic element in the model shows radiating mouth bars at the terminal end of a distributary channel comparable to the Wax Lake example. (e) Proximal to distal cross-section of the river-dominated delta model used (along the dip axis). Lithic predictions defined in this study are for the area within the dashed polygon (see also the limits of the map view area show in c). Basin subsidence rate was constant through time but increased from proximal to distal areas within the basin. (f)–(i) Horizon slices through the delta model, parallel to the mean dip of the preserved regressive–transgressive sandstone layer (f) & (g) or along a horizontal plane (h) & (i). The dashed line crossing (f)–(i) shows the location of the cross-section in (e).
Although standard training images may capture the key rock property spatial relationships in many cases, it is difficult to demonstrate the geological realism of simple geometrical models for application to specific types of reservoir deposits. In some cases, when the modelled lithic bodies have reasonable scaling, their spatial organization and connectivity relationships may be erroneous. In other cases, lithic variation scales not defined well in a training image may have critical control on reservoir performance (e.g. narrow threads of coarse-grained lag deposits or shale beds that have low volume but are laterally extensive). It is also possible that complicated realistic-looking training images that are not generated by statistical methods may not contain enough volume or the repetitive stationary 3D facies patterns needed to define MPS summations with distinct spatial organization. Given that MPS has become a standard method for heterogeneity prediction within reservoirs, there are surprisingly few studies that attempt to demonstrate that the geometrical training images typically used for modelling reservoirs are adequate to predict realistic subsurface fluid-flow patterns.
This study starts with a process-based depositional model that defines multiple nested scales of lithic variation formed during regressions and transgressions of a river-dominated delta. Geostatistical simulations are constructed to replicate these lithic patterns using standard reservoir modelling methods and various amounts of conditioning data extracted from the process model. Subsurface flow simulation of a line-drive flood displacement across reservoir models constructed from the depositional model and the different geostatistical simulations are compared to measure the impact of modelling assumptions and the value of conditioning information on reservoir behaviour forecasts. The goal of the study is to determine if simple variography with SGS or geometrical model training images with MPS are adequate to define complex reservoir variations and the value of different types of conditioning data in reservoir modelling.
River-dominated delta deposit
A physics-based depositional process modelling program (CompStrat) developed within Chevron uses depth-parameterized Navier–Stokes equations to predict flow and sediment transport of mixed grain sizes (following computations described in Parker et al. 1986 and Parker 2008). For this study, CompStrat was used to predict preserved lithic variations within a river-dominated delta deposit formed during a succession of regressions and transgressions driven by variations in sea level. The model was constructed assuming a constant mix of sand and mud entered a basin through a fixed inlet and expanded in area over the initial 5 km into a 5 × 8 km rectangular area (Fig. 1c). Subsidence rate within the basin increased from proximal to distal areas like that expected on a passive continental margin. Although cyclic sea-level variations were modelled to force the shoreline to regress and transgress the basin, sediment input and subsidence rates were chosen such that the shoreline never reached the distal end of the model. The resulting deposits are a relatively complicated succession of sandstone and mudstone layers with internal channel and mouth-bar deposits.
Details of these modelled deposits are described elsewhere (Amaru et al. 2017; Willis and Sun 2019), and so will be summarized only briefly here. The model predicts channels that pass sediment across the exposed delta top and deposit laterally-stacked mouth bars along the shoreline. When deposition was located on a drowned delta top area (following flooding by rising sea level), distributary channels frequently bifurcate and avulse as they pass into relatively thin mouth-bar deposits. During sea-level lowstands, a main distributary channel is entrenched across the delta top and thicker mouth-bar deposits formed along the prograding delta front edge. The delta front edge continued to build basinwards over successive sea level cycles. Lithic variations within the deposits showed repetitive patterns defined by sand-filled channels that pass down-dip into a localized cluster of terminal mouth bar sands. The lithic patterns are not stationary, in that the channel deposits changed in shape (width and thickness) both from proximal to distal areas of the model, and, at a given location, during regressions and transgressions of the shoreline. Similarly mouth bar sands tended to be smaller and less clustered when deposited on drowned delta top areas and larger and more clustered where fed by a lowstand channel entrenched to the delta front edge.
Depositional features observed on the CompStrat model surface at any one time can be compared directly with those observed on the surface of modern analogue systems. The main channel entering the model bifurcates down basin. Individual distributaries deposit a radiating cluster of mouth bars at the shoreline that together comprise a lobe-shaped shoreline protrusion (Fig. 1c and d). These features are comparable to Wax Lake delta within the larger Mississippi River delta system (Fig. 1a and b), which is similarly composed of a central distributary channel attached to a branching lobate set of mouth-bar storeys. Such lobe deposits define a reservoir element, in that the internal sands tend to be connected over these scale depositional units. For the Wax Lake example, the element is about 12 km across and individual mouth-bar storeys that comprise the element are about 2 × 3 km. In the CompStrat model used here, a similar depositional element is about 1 km in diameter and individual storeys are about 100 × 500 m (Fig. 1d). A few isolated finger-like mouth-bar storeys along the submerged lowstand delta are slightly longer (1.2 km). Geostatistical reservoir models typically incorporate these types of depositional dimensions to define expected distances for the spatial correlation limits of reservoir properties within a field.
This study examples lithic variations within a 2 × 3 km by about 15 m-thick block cut out of the larger CompStrat model (Fig. 1c and e: see the dashed polygon showing the map view and cross-section position, respectively). Grain-size variations across this sub-volume are shown in sandstone layer-parallel (Fig. 1f and 1g) and horizontal slices (Fig. 1h and 1i). See the inset in Figure 1e for the locations of these deposit slices. These slices are shown at the same scale as the detailed topographical surface map in Figure 1d. On average, older beds dip more basinwards because subsidence rate increased offshore. Many of the lithic bodies observed in the deposit are larger and more dip-elongate than the depositional features observed on the surface of the model. This reflects that lithic bodies are emergent features formed as channels migrate and mouth bars prograde over time. The size differences between the observed depositional surface features and the preserved lithic bodies suggest caution in using the dimensions of surface geomorphic bedforms directly as inputs for reservoir property modelling.
Predicting lithic variations based on sparse subsurface data
The CompStrat model defines lithic variations at the centre of 20 × 20 m plan view grid blocks for 2800 depositional surfaces. Grid block thickness, which reflects local sediment accumulated during an increment of time, varies from essentially zero to about 1 dm (Fig. 2a). Because standard geostatistical simulation is done on grids with blocks of more uniform volume, grain-size values were copied by near-point assignment to a ‘pixelated’ grid (Fig. 2b), with blocks uniformly 20 × 20 m by 0.05 m thick. This model contains about 70% net sand (Fig. 2c and d), such that the sandbodies are completely connected across the volume (Fig. 2d). Predicted grain size varied between 1 and 300 μm (mud to medium sandstone: Fig. 2e). Permeability (0.001 mD–1 D: Fig. 2f) and porosity (0.001–0.27) were calculated from model grain size based on empirical functions fitted to core-plug data from a conventional deltaic reservoir operated by Chevron. A vertical semi-variogram of the model grain size shows a cyclic structure with repetition of about 5 m, which defines the average thickness of successive regressive–transgression layers (Fig. 2g). Horizontal semi-variograms of grain size reflect variations across laterally stacked channel and mouth-bar deposits (cf. Figs 1f–i and 2h). In the strike direction, the horizontal semi-variogram shows a range between 0.5 and 1 km, reflecting the average width of the channel and mouth-bar lithic bodies (Fig. 2h). The lack of a semi-variogram sill in the dip direction reflects an increase in the mud/sand ratio across the whole model from proximal to distal areas (a non-stationary trend).
Details of the river-dominated delta model used in this study. (a) Depositional grid that follows bedding surfaces (every 10th grid line is shown by a black line). (b) Uniform pixelated grid with rectangular grid blocks. The pixelated grid sampled grain size from the depositional grid by grid-block centre near neighbour assignment. (c) Model gross thickness. (d) Model net thickness (sand thickness). (e) Model cumulative frequency of grain size. (f) Model cumulative frequency of permeability. (g) Vertical semi-variogram of log base 2 grain size. (h) Horizontal semi-variogram of log base 2 grain size in the dip (red) and strike (blue) directions.
Facies are defined in the model based on grain-size classes (facies 1 is mud, facies 2 is very fine sand, facies 3 is fine sand and facies 4 is medium sand). Maps of the vertically averaged facies proportion in the model (Fig. 3) show a proximal to distal decrease in facies 3 (fine sand) relative to facies 1 (mud). Medium sand, restricted to the basal part of larger channels, also decreases in abundance basinwards. There is a slight increase in better quality facies upsection, and more pronounced local variations across individual regressive–transgressive layers.
Facies proportions maps show proximal (left) to distal (right) changes in the grain-size-defined facies.
Five injector wells are placed along the distal end of the model and five producer wells along the proximal end (Fig. 4a). The wells penetrate better reservoir facies than average (cf. Fig. 4a and b), reflecting that a geologist might bias well placement to improve field development. Even if there were no a priori subsurface information, so that the wells were placed randomly within the model, there would be uncertainty in the facies proportion estimates from wells (Fig. 4b1 and b2). Although potential variations in estimates of net sand decrease with the number of wells, a Monte Carlo simulation of random wells placed in this model indicated that a 10% uncertainty in net sand estimates remained after 25 wells are placed (average of 60 acre well spacing). Varying facies proportions could have significant impact on subsurface flow predictions.
(a) Distribution of injector and production wells. (b) Vertical changes in the facies proportion within the model (b1) and of facies contacted by the wells (b2). The changes in the proportion of facies 1 (mud), facies 2 (very fine sand), facies 3 (fine sand) and facies 4 (medium sand) reflect regressions and transgression of the deltaic shoreline and compensational stacking of larger-scale deltaic lobes within an area larger than the model considered here. (c) Variation of estimates of net/gross with a varying number of randomly placed wells based on thousands of Monte Carlo well-placement realizations.
Geostatistical simulations were constructed based on two sets of scaling relationships: (1) the dimensions of the channels and mouth bars observed on topographical surface maps of the modelled deposits (Fig. 1d: short-correlation-length models); and (2) the dimensions of preserved lithic bodies observed in horizontal slices of the modelled deposit (Fig. 1f–i: long-correlation-length models). For SGSs, observed dimensions were used to define strike, dip and vertical semi-variogram correlation lengths. For MPS, a Chevron proprietary training image generator was used to build simple geometrical models that looked like the deposits (Fig. 5). This program constructs a stationary representation of the deposits from a succession of nested ellipses with concave, convex, or flat top and base. Statistical ‘rules’ define local variations in facies-body size and nesting relationships (see the caption to Fig. 5 for an example of the rules applied). Similar MPS training image generators are now available within most commercial reservoir modelling software. Although many different long- correlation-length and short-correlation-length MPS models were constructed using this method, only two representative long and one short are discussed here. For comparison, the pixelated copy of the CompStrat model was also used directly as a training image.
Training images constructed using a geometrical modelling system were selected to be reasonable approximations of variations observed with the CompStrat delta depositional model used in this study. (a) Conceptual strike-orientated cross-section of facies geobodies defined within the training image shown at about ×100 vertical exaggeration (V.E.). Facies 1 (blue) is the non-reservoir background facies. Facies 2 are upward-concave elliptical bodies resembling mouth bars (green), whereas facies 3 and 4 are downward-concave bodies with more channel-like geometry (yellow and red, respectively). (b) Short-range training image. Body dimensions were estimated by looking at the surface depositional features of the CompStrat model. The dimensions of facies 2 are defined as 0.5–1 km long, 0.25–0.5 km wide and 0.5–1 m thick bodies elongate along a depositional dip trend (randomly assigned to be between ±20° from east to west). Facies 3 are defined to vary from 2 to 4 km long, 25–50 m wide and 2–5 m thick, and are generally enclosed within facies 2 (minimum 50% overlap). Facies 4 are defined to be downward-convex elliptical bodies varying from 0.5 to 1 km long, 10–20 m wide and 0.5–1 m thick that are enclosed within (minimum 95% overlap), and along the base of, facies 3. (c) Long-range training image dimensions were defined by examining variations with the deposits preserved within the CompStrat depositional model. The same facies relationship rules as in (b) were used but dimensions for facies 2 (3–6 km long, 0.5–1 km wide and 1–2 m thick), facies 3 (2–4 km long, 150–500 m wide and 2–5 m thick) and facies 4 (500–750 m long, 25–75 m wide and 1–2 m thick) were different. (d) A second, long-range training image is as in (c) but facies 2 bodies are wider and thinner, and facies 4 bodies are wider and thicker. The training images were four times larger than the volumes shown, which are displayed here cropped to match the dimensions of the final MPS models. The global facies proportions within the MPS models are not defined by the training images but, rather, match that in the CompStrat model. (b)–(d) are shown at shown at ×100 vertical exaggeration.
Geostatistical simulations were performed within two different grids, both with 20 × 20 m map view blocks but different block thicknesses (Fig. 6). The ‘uniform’ grid had rectangular blocks, 0.05 m thick, like the pixelated copy of the CompStrat model (Fig. 6c). For the second ‘stratigraphic grid’, five surfaces were correlated between the wells based on recognized maximum flooding surfaces. The correlated surfaces were simpler than the actual bedding surfaces within the model, which were complicated by local incision of younger channels (cf. Fig. 6a and b) but these surfaces generally aligned similar age stratigraphic layers. The stratigraphic grid had the same number of grid blocks as the uniform grid but block thicknesses varied proportionally between the correlated surfaces (Fig. 6d).
(a) Sea-level-defined stratigraphic units within the delta depositional model delineate complex lithic variation trends within the deposits. (b) Surfaces correlated between the wells show the stratigraphic variations that could reasonably be corrected across the model. (c) Stratigraphic variations shown on a uniform grid cross-section (rectangular gird blocks). (d) Stratigraphic variations shown on a stratigraphic grid defined by layers spaced proportional to the separation of the correlated surfaces (same square map view grid-block geometry as for the uniform grid). Note that it is only the grid geometry than has changed between (c) and (d) (and slight differences in the viewed orientation of the models) but no differences in the physical geometry of the stratigraphic units.
MPS simulations were conditioned to three different levels of data extracted from the depositional model: (1) wells only; (2) wells with global vertical proportion curves; and (3) wells with global facies proportion volumes constructed by combining the vertical proportion curve with facies proportion maps (cf. Figs 3 and 4). In all cases, the global facies proportions were held constant to values in the full-resolution CompStrat model, so that impacts of different levels of conditioning reflected variations in facies organization rather than simply differences in relative facies abundance. In addition, a post-processing method run on some models used local block reassignment to increase the horizontal continuity of facies bodies or, alternatively, to increase both vertical and horizontal continuity of facies bodies (Strebelle and Remy 2005). This post-processing did not change the global facies proportion within the models.
The SGS models defined grain-size variations directly, without first defining facies classes. The MPS models were converted from facies to grain-size volumes by randomly adding the volume-corrected grain-size distributions defined by grain-size mean and variance of each facies class within the full-resolution CompStrat model (using very short range, 60 m, SGS). Grain size was then used to define permeability and porosity variations. A compilation of 60 models was used in a designed reservoir flow simulation study. The base case for comparisons to the geostatistical simulations is the full-resolution CompStrat model gridded along depositional bedding. The pixelated version of the CompStrat model is included to determine the impact of the fine-scale re-gridding like that used in the geostatistical models. SGS models used short- correlation-length and long-correlation-length variograms applied variably to the uniform and stratigraphic grids (Fig. 7). MPS models were defined by image training on the pixelated copy of the depositional model: the two long-correlation-length geometrical training images or the short-correlation-length geometrical training image (Fig. 5). Each MPS training model result was applied with reference to one of three levels of data conditioning and different MPS post-processing on both the uniform grid and stratigraphic grids (Figs 8 and 9).
Examples of grain-size pattern predictions based on sequential Gaussian simulation (SGS) across stratigraphic and uniform grids. (a) Pixelated CompStrat model shown for comparison. (b) SGS with longer variogram range, uniform grid. (c) SGS with shorter variogram range, uniform grid. (d) SGS with longer variogram range, stratigraphic grid. (e) SGS with shorter variogram range, stratigraphic grid.
Examples of multipoint statistics (MPS) grain-size pattern predictions on the uniform grid with varying amounts of conditioning data. Correlation structure defined by: (a) the short-range training image (Fig. 5b); (b) the long-range training image (Fig. 5c); and (c) the CompStrat depositional model used as the training image. Cases were conditioned using: (1) wells only; (2) wells and the vertical proportion curve; and (3) wells and the facies proportion volume.
Examples of multipoint statistics (MPS) facies predictions on the proportional stratigraphic grid with varying amounts of conditioning data. See the caption to Figure 8 for an explanation of the cases.
Comparing lithic variation predictions
The full-resolution and pixelated versions of the CompStrat depositional model look very similar (Fig. 2a and b). The pixelated grid (Fig. 7a) has fine enough resolution (6 × 106 grid blocks) to preserve all the storey-scale lithic bodies and most of the finer-scale internal bedding variations.
The SGS models were conditioned to the well data, and used variogram lengths like those observed for the delta surface features (Fig. 7b and d) or based on the observed lithic variations within the CompStrat deposit (Fig. 7c and e). Although the results matched at the well locations, as the method ensured, grain sizes become more diffuse away from the wells. The range of different-sized channel and mouth-bar interfingering was not replicated well by these models. The pronounced layering within the CompStrat model is defined better in the stratigraphic proportional grid than in the uniform grid. Because the wells penetrated more coarse-grained deposits than average, and the final facies proportions were specified as input conditions, the dominate pattern predicted is a gradual lateral fining away from the wells.
The short-range MPS training images (Figs 8a and 9a) are based on the observed dimensions of the depositional bodies on the final depositional surface of the model, like that which could be observed in an air photograph (Fig. 5b). The resulting MPS models are dominated by local variance, giving models a weakly structured salt and pepper texture. Despite the diffuse short-range variations, there are distinct coarser-grain-size bodies that are elongate along dip (reflecting the preferential correlation direction imposed by the training image). Adding the vertical facies proportion curve produced a weakly defined layering in the stratigraphic grid cases but made little obvious difference in the uniform grid models. This contrast reflects that the vertical proportion curve averages horizontally across stratigraphic units lower in the volume that are slightly included due to long-term downdip increases in subsidence rate. For the stratigraphic grid, the proportional grid layering between the main flooding surfaces keeps the vertical proportion curve constraints better aligned with the stratigraphic layers. Models with horizontal proportion maps as added constraints correctly predict the gradual decline in grain size downdip. Post-processing the MPS results makes the boundaries of coarse-grained bodies slightly less diffuse but the results are visually subtle (and, thus, are not shown here).
The models constructed with long-range MPS training images (Figs 8b and 9b) show more distinct layering with discrete coarser- and finer-grained bodies. In both dip and strike, the imposed spatial correlation length was greater than the well spacing, and thus the models are locally more layered than the CompStrat model in intervals where lithic variations in adjacent wells fortuitously align. Adding a vertical facies proportion curve (Figs 8a2 and 9a2) strengthens delineation of the larger-scale stratigraphic layering, and adding a facies portion map improves the reproduction of the proximal to distal grain-size fining. As for the other models, the proportional stratigraphic model better defines lithic trends across the slightly inclined stratigraphic intervals lower in the model.
The results using the CompStrat model directly as a training image were broadly like that generated using the long-range MPS training images. The defined variations are more organized than the models produced using short-range training images but slightly more diffuse than those using the long-range training images (Figs 8c and 9c). Global MPS geobody definitions are presumably weakened by conflicting patterns in non-stationary areas of the depositional model. The example on a uniform grid with full facies volume constraints is distinctly more diffuse. This may indicate that applying too many constraints on the regionalized summations can also weaken MPS-defined patterns. For example, the MPS geobodies in the uniform grid defined by the full facies proportion volume (Fig. 8c3) show significantly more diffuse patterns than the example conditioned only to wells (Fig. 8c1).
Reservoir flow simulations
A commercial black oil flow simulator (Schlumberger Intersect™) was used to model the displacement of oil by a water flood between the line of injector and producer wells (Fig. 4a) for the depositional and geostatistical models. Production rates, injector pressures, well completions, relative permeability, pressure–volume–temperature properties, irreducible saturation correlations and other dynamic parameters are kept constant in the flow simulations. Simulations are simple two-phase water displacing oil with no gas coming out of solution at reservoir conditions. Since all dynamic parameters are similar, contrasts in subsurface fluid flow can be attributed directly to static rock property model differences.
Models were compared based on two production metrics: initial water breakthrough time of fluids from injector to producer wells (Fig. 10); and oil recovery after 0.3 pore volumes of water were injected (Fig. 11). Other metrics, like pressure at wells and flow patterns through the grid, were also monitored for model comparisons. The full-resolution CompStrat depositional model is used as the base case for comparisons. The pixelated version of the CompStrat model resulted in similar predictions of water breakthrough time (Fig. 10) and only slightly lower recovery (Fig. 11), indicating that variations in the distribution of rock properties within the different models had substantially greater impact than contrasts in the grid structure.
Simulation results defined by differences in water breakthrough time. Results for models defined on the uniform and stratigraphic proportional grid are displayed separately. The dot–dash line shows the results of the full complexity depositional model grid (CompStrat), which is used as a reference for results based on the pixelated CompStrat model: and (a) MPS models trained on the CompStrat Model (Fig. 7a); (b) long-range MPS models trained on the training image in Figure 5c; (c) long-range MPS 2 models trained on the training image in Figure 5d; (d) short-range MPS model defined by channel and mouth-bar dimensions observed on the surface of the depositional model (Fig. 5b); and (e) sequential Gaussian simulations (SGS). For each geostatistical model specification, cases are compared conditioned on: (1) wells only; (2) wells and vertical proportion curves; and (3) wells and vertical proportion volume. For the MSP models, results are shown for those without post-processing (yellow), with post-processing to increase horizontal facies-body continuity (green), and with post-processing to increase both vertical and horizontal facies-body continuity (blue). For the SGS cases, results are based on variograms of different correlation length: long (yellow), medium (green) and short (blue).
Simulation results of models defined by differences in recovery after 0.3 pore volumes injected. See the caption to Figure 10 for an explanation of labels and symbols.
The SGS and short-range MPS models, both with visually more diffuse lithic variation patterns than the others, had significantly longer breakthrough times and higher predicted recovery than the base case. Adding trend data (vertical proportion curve or the full facies proportion volume) improved predictions somewhat but results remained significantly worse than the longer correlation-length models that referenced the preserved geobody dimensions within the depositional model. Although post-processing MPS models to locally increase geobody continuity typically improved estimates of recovery, water breakthrough time predictions generally became worse. Differences between the pixelated uniform grid and the more geologically informed stratigraphic grid were subtle for these cases, with correlation range shorter than the well spacing.
Predictions based on the stratigraphic grid were significantly better than those based on the uniform grid for the longer correlation-length MPS models (where predicted correlation lengths were similar to the well spacing). This was not surprising, as simulations on the stratigraphic grid visually better preserved the larger-scale regressive–transgressive defined layering (cf. Figs 8 and 9). Adding trend data generally improved predictions but less dramatically than for the simulations based on the short-range geostatistical models. This suggests that trend data become less important as lithic correlation distances increase relative to well spacing. Post-processing MPS models to locally increase geobody continuity generally improved simulation predictions based on the uniform grid but in many cases made predications based on the stratigraphic grid worse. In cases where the MPS simulation predicted lithic patterns well, post-processing to enhance the continuity of lithic patterns predicted breakthrough times that were too early and recovery that was too pessimistic.
Using the CompStrat depositional model directly as a training image for MPS simulation did not result in significantly better subsurface flow predictions than a geometrical training image that was constructed to reproduce observed variations in the depositional model (cf. CompStrat training with long-range models in Figs 10 and 11). Although the depositional model defines the ‘true’ lithic variations, non-stationary patterns were averaged during MPS training image sampling, which in turn probably degraded local pattern replication. If the simulated area was smaller relative to the depositional model, such that lithic patterns remained more stationary across the area of interest, then MPS training directly on variations within the depositional models might be better than training on a stationary geometrical representation of the lithic patterns.
Differences between the geostatistical models can also be observed when comparing water-displacement patterns within the simulation models (Fig. 12). Map view slices through the middle of simulation models showing water saturations after 2 years of water injection were used to compare production patterns. Short-range models are more homogeneous at the scale of well separation, resulting in relatively uniform displacement that produced slow water breakthrough and high recovery (Fig. 12a). Longer-range correlation models result in a more irregular displacement, and heterogeneities allow flows to progress faster along some pathways relative to others (Fig. 12b and c). Flows through the full complexity lithic patterns defined by the depositional CompStrat model are even more irregular and have more pronounced 3D lithic variations that cause flow to pass vertically between the simulated grid layers (Fig. 12d). With added trend data (Fig. 12e and f), the short-range models better represent displacement patterns within the non-stationary deposit but the flood displacement front remains too regular and, thus, recovery predictions are too high.
Displacement patterns away from wells through different models at 0.3 pore volumes injected. Each example shows a horizontal slice in the same vertical position through the respective model. (a) MPS model based on the short-range training image and conditioned to well data. (b) MPS model based on the long-range training image and conditioned to well data. (c) MPS model trained on the pixelated CompStrat model and conditioned to well data. (d) Full complexity grid CompStrat model. (e) MPS model based on the short-range training image, and conditioned to both well data and the vertical (v) facies proportion curve. (f) MPS model based on the short-range training image, and conditioned to both well data and the full facies proportion volume. All models have the same global facies proportions, and MPS simulation was applied to the stratigraphic grid.
Discussion
This designed reservoir simulation study suggests that geostatistical methods can define reasonable reservoir models of river-dominated delta deposits despite complex stratigraphic layering with non-stationary internal lithic patterns. Models that define correlation patterns based on the size and geometry of surface geomorphic features (e.g. in this case, the static surface expression of channels and mouth bars) tend to underpredict the spatial continuity of lithic patterns. As discussed earlier in this paper, the true geometry of lithic patterns can only be defined by considering how surface features grow, migrate and avulse over time, or from direct observations of the preserved deposits. Applying supplemental conditioning data to MPS simulations might not compensate for inaccurate initial estimates of correlation length. For example, Figure 13a shows the predicted water production rate estimated using the short-correlation-length training image designed to match surface depositional features. Although supplemental conditioning data improve these estimates, even simulations supported by a complete facies proportion volume estimate delayed water breakthrough relative to the full complexity depositional model base case. The reason for the delayed water invasion of the production wells is clear when comparing differences in the predicted lithic patterns. With more conditioning data the large-scale facies patterns in the MPS simulation become more like the depositional model but the fine-scale variations remain more diffuse. Compare MPS simulation predictions with progressively more conditioning data (Fig. 9a1, a2 and a3, respectively) with grain-size variation patterns in the full complexity depositional model (Fig. 7a).
Change in predicted reservoir development behaviour defined by the water production rate for (a) a varying level of conditioning data and (b) varying estimates of net/gross. In all cases, MPS models used the short-correlation-length training image and the stratigraphic grid.
All of the geostatistical simulations presented here had the same facies proportions as the CompStrat depositional model (i.e. and thus the same global grain size, permeability and porosity distributions). For real field assessments, uncertainty of reservoir volume (e.g. net/gross estimates, structural closure and original fluid contacts) can have a significant impact on subsurface flow patterns (cf. Figs 4c and 13b). If the range in lithic pattern variations does not have a first-order impact on reservoir behaviour predictions of a given field, then details of how these variations are assigned become less important. The range of results for the different models discussed here is rather low (e.g. recovery varies only between 10.3 and 11.3% after 3 years). This reflects, in part, that despite complex stratigraphic layering, the models are not particularly heterolithic: 87% of the deposit is sand, the sands are 100% connected and bed-scale heterogeneity is not included (constant river discharge over time). Subsurface flow predictions through deposits with lower net/gross, more internal isolated compartments and more pronounced nested scales of heterogeneity may be more sensitive to the details of the spatial patterns of lithic variations predicted in the static reservoir model. More studies like the one presented here are clearly needed to determine the sensitivity of reservoir production forecasts to facies patterns characteristic of a broader range of depositional environments and to determine how well these different types of rock property variations are replicated by different reservoir modelling techniques.
Simulation of river-dominated delta deposits for subsurface prediction is arguably simpler than facies patterns observed within deposits formed in some other types of depositional environments. Although channel and mouth-bar sands changed in both size and proportion across the simulated volume, and within the stacked regressive–transgressive stratigraphic units, the general pattern of coarser-grained channel deposits expanding into basinwards fining lobes remained. Therefore, in this case it is arguably reasonable to define the deposits as stationary for geostatistical summation, despite observed spatial variations in delta deposit character. This might not always be the case. For example, stacked wave-dominated shoreface deposits might be characterized by distinct belts of shoreline sands adjacent to parallel belts of back-barrier heterolithic deposits with completely different lithic character. In such cases, facies belts would need to be defined correctly first, before modelling of the finer-scale internal lithic patterns within each of the facies belts.
Although results from this study cannot be directly applied to geostatistical modelling of rock property patterns formed in all depositional settings, some learnings have broader application. Process-based depositional models can provide 3D analogues to better define preserved lithic patterns in different depositional and stratigraphic settings. Variations predicted by depositional models can help to constrain correlation distances and nested facies relationships for geostatistical simulations. Reservoir simulation of production scenarios using depositional models can be compared with results of geostatistical simulation to characterize uncertainty of results and likely ranges of production forecasts. Such comparisons can also validate that geostatistical methods and input parameters are appropriate for specific reservoir modelling applications. For example, in cases where simple stationary models fail to match depositional model cases, the added work of applying more complex modelling techniques would be earlier to justify (e.g. nested MPS models or object-MPS hybrid models).
Conclusions
Physics-based depositional process models provide 3D digital analogies to define spatial patterns of subsurface rock property variations for reservoir characterization and subsurface flow simulation.
Rock property correlation dimensions and training images for geostatistical simulation need to reference the preserved deposits rather than observed static land-surface features. This reflects that subsurface lithic patterns gradually emerge as landforms grow, migrate and avulse over time.
Multipoint geostatistical models with key depositional model input and conditioning did a reasonable job in defining heterogeneity patterns in a high net/gross river-dominated delta deposit.
Less-structured geostatistical models tend to overestimate recovery and water breakthrough time.
Trend constraints on MPS simulations are most important when shorter correlation lengths are used, and the grid layering does not follow stratigraphic surfaces. If the variogram range spans the wells and the geological guide surfaces incorporated into grid construction correctly align correlations, then additional constraints become less important for predicting global facies patterns.
Using a depositional model directly as a training image is not clearly better than using a geometrical stationary training image generator with reasonable input dimensions. The averaging of non-stationary patterns in the depositional model can degrade global MPS pattern definitions.
Post-processing to improve geobody continuity makes modest improvements when initial geostatistical models predictions are weakly structured but can degrade results when geostatistical simulation already adequately predicts the general lithic patterns.
Process-based depositional models can be used to optimize geostatistical model methods and input parameters to reduce uncertainly of reservoir simulation results for development planning. Deltaic facies patterns defined by a river with constant discharge can be replicated well by MPS simulations. Additional studies like the one presented here are needed optimize geostatistical model methods for deposits formed under more variable depositional conditions and for other depositional environments.
Acknowledgements
We thank Chevron Energy Technology Company for permission to publish this work. Sebastien B. Strebelle and Maisha Amaru provided constructive reviews of the initial manuscript.
Funding
This research was funded in toto by Chevron.
Author contributions
BJW: conceptualization (lead), formal analysis (equal), methodology (equal), writing – original draft (lead); SK: formal analysis (equal), methodology (equal), writing – review & editing (equal); TS: methodology (equal), project administration (lead), software (lead), visualization (equal).
- © 2020 Chevron U.S.A, Inc.. Published by The Geological Society of London for GSL and EAGE
This is an Open Access article distributed under the terms of the Creative Commons Attribution 4.0 License (http://creativecommons.org/licenses/by/4.0/)
Defining lithic patterns within river-dominated delta deposits for geostatistical simulation
