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.

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).

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.

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.

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.

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).

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).

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.

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.


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).

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).


  • 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.


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.


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).

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/)