# Reproduction of channel stacking patterns in geomodeling: Metrics and impact of the modeling strategy on reservoir flow behavior

Enrico Scarpa, Pauline Collon, Irina Panfilova, Guillaume Caumon

► To cite this version:

Enrico Scarpa, Pauline Collon, Irina Panfilova, Guillaume Caumon. Reproduction of channel stacking patterns in geomodeling: Metrics and impact of the modeling strategy on reservoir flow behavior. Marine and Petroleum Geology, 2025, 173, pp.107235. 10.1016/j.marpetgeo.2024.107235. hal-04873487

HAL Id: hal-04873487

https://hal.univ-lorraine.fr/hal-04873487v1

Submitted on 8 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Research paper

# Reproduction of channel stacking patterns in geomodeling: Metrics and impact of the modeling strategy on reservoir flow behavior

Enrico Scarpa<sup>a,*</sup>, Pauline Collon<sup>a</sup>, Irina Panfilova<sup>b</sup>, Guillaume Caumon<sup>a,c</sup>

## ARTICLE INFO

**Keywords:**

Channelized turbidite systems

Stacking patterns

Stochastic simulations

Hydraulic connectivity

Reservoir simulation

## ABSTRACT

Channelized turbidite systems are often grouped into complexes and exhibit various stacking patterns, which play a crucial role in controlling the connectivity between high-permeability and low-permeability sedimentary bodies. While some studies have analyzed the static connectivity of different stacking patterns, few have quantitatively evaluated the dynamic implications of these patterns on fluid flow circulation at large scale, and no method have started to compare the outcome of different modeling methods. This study addresses this gap by quantitatively investigating the impact of different classes of geostatistical modeling methods on both static and dynamic connectivity using multiple metrics. Using an original object-based simulation method, we have stochastically generated 3 different stacking pattern sets of 100 realizations each with similar facies proportions, that represent: (i) disorganized channels, (ii) independent channels conditioned to a vertical sand proportion map, and (iii) organized stacking that reproduce vertical and lateral migration of channels. Channel internal heterogeneities are neglected and homogeneous properties are attached to the three main facies: channels, inner levees and outer levees. To analyze the hydrodynamic responses, we set up a two-phase system of oil and water. Reservoir simulations are performed for all 300 realizations, providing production curves, saturation field and pressure field. Classical dynamic metrics (breakthrough time and recovery efficiency) are completed by original ones, designed to understand what governs the differences between the sets: at 30% of the Pore Volume Injected, saturation front shapes are extracted and we compute their surface area and sphericity to get objective comparison criteria. Results reveal significant differences in flow behavior across the three studied sets of stacking models: Disorganized stacking patterns, and, to a smaller extent, the conditioned disorganized stacking patterns, exhibit delayed water breakthrough times and more optimistic recovery than organized stacking patterns, which have, however, a higher static connectivity. This unexpected result seems to be related to gravity-driven fluid segregation in the considered reservoir production settings. Overall, these results quantitatively confirm that the ability of the channel simulation method to generate realistic stacking pattern is essential. In addition, the lack of a clear relationship between static and dynamic connectivity metrics suggests the need for further research to develop more effective forecasting techniques for reservoir behavior in turbiditic channel settings.

## 1. Introduction

Turbidite depositional systems are widely encountered along continental margins and comprise a large fraction of the shallow and deep marine sediments. In these systems, the reservoir rocks comprise meandering channels characterized by high longitudinal continuity. Channels are gathered into complexes and display stacked organizations arising from lateral migration and aggradation during deposition (Talling et al., 2015; Peakall et al., 2000; Labourdette et al., 2006). Despite decades of studies focusing on understanding the link

between depositional environment and meandering geometries, numerical reproduction of channelized systems while honoring subsurface data has proven challenging (Wang et al., 2018; Shmaryan and Deutsch, 1999; Li et al., 2023). This difficulty arises from the three-dimensional complexity and diversity of scales in the available data, that comprise seismic reflection images, and a limited amount of cores and well logs (Mayall et al., 2006; McHargue et al., 2011; Labourdette and Bez, 2010; Lemay et al., 2020).

<footnote><sup>a</sup>GeoResources, Université de Lorraine - CNRS, 54000, France</footnote>

<footnote><sup>b</sup>LEMTA, Université de Lorraine - CNRS, 54000, France</footnote>

<footnote><sup>c</sup>IUF, Institut Universitaire de France, 75000, France</footnote>

<footnote>* Corresponding author.
E-mail address: enrico.scarpa@univ-lorraine.fr (E. Scarpa).</footnote>
When high resolution seismic images are available, it is possible to directly interpret and draw the channel meanders (Labourdette et al., 2006; Covault et al., 2021; Morris et al., 2022). However, limited seismic frequency at depth calls for stochastically simulating channel facies locations to assess uncertainties below seismic resolution (Pyrcz and Deutsch, 2005; Larue and Hovadik, 2008). For this purpose, object-based methods, that represent a single channel with a parametric object, are often used as they guarantee the longitudinal continuity of channels, hence better reproduce reservoir connectivity (e.g., Deutsch and Tran, 2002; Lopez, 2003; Larue and Hovadik, 2008; Rongier et al., 2017b; Wang et al., 2018). A classical object-based approach is the Boolean model, which generates multiple independent channels or channel complexes, to simulate a complete depositional environment (e.g., Haldorsen and Damsleth, 1990; Wietzerbin and Mallet, 1994; Deutsch and Wang, 1996; Holden et al., 1998; Shmaryan and Deutsch, 1999; Viseur, 2004; Pyrcz and Deutsch, 2005; Bertoncello et al., 2013; Rongier et al., 2017b). Nevertheless, mutually independent channels do not accurately reproduce the stacking patterns resulting from the depositional processes (Labourdette, 2008). Physical-based methods (e.g., Teles et al., 1998; Cojan et al., 2005; Karssenberg and Törnqvist, 2001; Karssenberg and Bridge, 2008; Lopez, 2003; Sylvester et al., 2011) and rule-based methods (e.g., Wen, 2005; Labourdette, 2008; Pyrcz et al., 2015; Rongier et al., 2017a; Parquer et al., 2017, 2020) consider channel migration vectors, hence permit to build consistent channel stacking patterns. These elaborate tools require, however, larger computational times and resources, and honoring data is still not a fully resolved question. In the frame of an uncertainty analysis, in an operational case, this gain of geological realism is, so, often questioned at the light of its impact regarding the reservoir appraisal and development plans.

Few works quantitatively compare several geostatistical approaches for modeling channelized complexes in term of flow responses. Sensitivity studies based on Boolean models (Larue and Hovadik, 2008; Jager et al., 2009) show the flow impact of various reservoir parameters such as facies proportion or individual geometrical channel parameters. The importance of internal petrophysical heterogeneities is also highlighted on small scale models consisting of one to three channel stacking configurations (e.g., Labourdette et al., 2006; Alpak et al., 2013). More stacking arrangements for two successive channels are considered in a cylindrical 3D case by Meirovitz et al. (2021). All these studies suggest the impact of channel stacking arrangement on reservoir behavior, but it is unclear to what extent these results hold in other settings, as some averaging effects could occur at larger scale and when considering more variability in the configurations. In this paper, our first goal is to quantify this impact at field scale on a large sets of channel configurations obtained with different modeling approaches.

Rongier et al. (2017a) compare such large ensembles of geostatistical realizations reflecting different modeling strategies, but only with static metrics. Our study is similar, but quantifies the dynamic impact of the modeling strategies by performing flow simulations on each realization. Because of the large computational cost associated with this type of sensitivity study, an additional motivation for this work is to establish to what extent ranking methods based on static connectivity indicators may be used to select only a small number of representative models for flow simulation. Many studies have been made in this area (e.g., Hird and Dubrule, 1998; Hovadik and Larue, 2007; Scheidt and Caers, 2009), but they concern a single geostatistical methodology. Our second objective is, therefore, to assess to what extent static and dynamic connectivity metrics can be related, and whether these relations are affected by the type of geostatistical simulation method.

In summary, this paper objective is twofold. First, we investigate to what extent the modeling strategy used at a multi-kilometric scale for the channel stacking patterns affects reservoir behavior using both static and dynamic reservoir metrics, and with large sets of stochastic realizations. Second, we seek relationships between these two sets of metrics and try to assess whether these relationships depend on

the modeling method. For this, we extend the approach introduced by Rongier et al. (2016) to examine three sets of turbidite reservoir models in sparse data settings (Section 2). Each set considers one particular stacking pattern setting as proposed in Rongier et al. (2017a). The connectivity is quantified with various static and dynamic metrics, some of which were specially designed for this study (Section 3). The recovery process is performed for a reservoir appraisal scenario (a reservoir with no historical data available, Section 4). Lastly, we explore the connection between static and dynamic metrics while also offering potential insights into the suggested workflow (Section 5).

## 2. Data set generation

### 2.1. Geostatistical simulation of stacked channels

We consider different stochastic scenarios corresponding to a typical subsurface architecture similarly to Rongier et al. (2017a). The geological context is a channel complex with its internal elements inspired from many actual turbidite systems, such as those of West Africa described by Labourdette (2007) or Deptuck et al. (2007). In this context, intra channel elements are located in a deterministic large incision (i.e., channel complex) as could be identified from a low-resolution seismic image. We generate three sets of 100 stochastic realizations reproducing different stacking conditions. Differently from Rongier et al. (2017a), we adapt the number of channel elements in all stochastic realizations in order to have approximately the same range of channel facies proportion, and the same original Oil-In-Place.

In the first set, each realization includes 28 independent channel elements confined within the complex. Each channel is generated by a Lindemayer-system based approach (Rongier et al., 2017b), and placed randomly in the simulation domain. These realizations are termed *disorganized stacking* (DS, Fig. 1a), as channels are positioned independently one from another.

The second set, termed *conditioned disorganized stacking* (CDS, Fig. 1b-c), uses the same simulation method, and additionally employs a three-dimensional sand probability map as typically obtained by seismic inversion. This map controls the placement of the channels within the complex while maintaining their mutual independence (Fig. 1c). Because the amalgamation rate is higher than in the DS set, 33 channels are necessary to keep the same overall facies proportion.

The final set, *organized stacking* (OS, Fig. 1d), models two hypothetical phases of channel evolution, starting with low aggradation and extensive lateral migration, and followed by high vertical aggradation and reduced lateral migration. An initial channel is generated at the base of the domain with the same Lindemayer-based simulation as before. The more recent channels are then iteratively obtained by a forward SGS-based lateral migration algorithm combined with a global aggradation (Rongier et al., 2017a). Consistently with the probability cube used in the CDS scenario, the OS set uses 27 lateral migration steps and 12 vertical aggradation steps (Fig. 1c).

### 2.2. Simulation grid and petrophysical properties

For each geostatistical realization, the modeling method generates sets of channel centerlines. Volumetric channel envelopes are built by Non Uniform Rational Basis-Splines (NURBS) volumes (Ruiu et al., 2016), whose widths vary between 130 m and 330 m and thickness vary between 6 m and 10 m. NURBS-channel volumes are rasterized into a Cartesian grid covering a domain of approximately 11 km × 3.4 km × 300 m composed of 840,000 cells (240 × 70 × 50, Fig. 2). The rasterization process is applied to all 300 reservoir models. Whereas the Cartesian grid leads to a stair-step approximation of the channel geometry, it is chosen over a channel-conformable grid. This choice indeed supports the cell face orthogonality necessary for two-point flux-approximation (TPFA) discretization, which is crucial for accurate flow computation (Wu and Parashkevov, 2009).
<figure><img src="image_3.png" /><figcaption>**Fig. 1.** Example of geostatistical realizations for: disorganized stacking (DS, a), conditioned disorganized stacking (CDS, c), computed from a sand probability cube (b), and organized stacking (OS, d). The probability sand map used in the modeling method is the same for all CDS realizations.</figcaption></figure>

In the reservoir models, the rasterized channel elements are considered as coarse sandy sediments. The rest of the channel complex is assumed to correspond to the inner levees, composed of an alternation of interbedded fine sand and shale. Outside of the channel complex, the

lithofacies correspond to shale-rich outer levees. For simplicity, these three classes of lithofacies are modeled with constant petrophysical values for porosity $\phi$, horizontal permeability $k_h$ and vertical permeability $k_v$ (Table 1, and Fig. 2). The vertical permeability is ten times
<table><caption>**Table 1**
Lithofacies of reservoir model.</caption><thead><tr><th>Unit</th><th>φ</th><th>k<sub>h</sub> [mD]</th><th>k<sub>v</sub> [mD]</th></tr></thead><tbody><tr><td>Channel</td><td>0.35</td><td>1000</td><td>100</td></tr><tr><td>Inner Levees</td><td>0.25</td><td>10</td><td>1</td></tr><tr><td>Outer Levees</td><td>0.1</td><td>0</td><td>0</td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Example of channelized reservoir model. Turdibite channels are rasterized into a Cartesian grid of 6 m × 46 m × 48 m (840,000 cells). The simulation grid cells are populated by three lithofacies corresponding to the channels, the inner levees, and the outer levees. Wells are horizontal to ensure they cross channels in all realizations.</figcaption></figure>

lower than the horizontal permeability to reflect subscale layering. The choice of effective constant values per lithofacies is consistent with the approach used in numerous studies on turbidite deposits, including Labourdette (2007), Stewart et al. (2008), Alpak et al. (2013), Meirovitz et al. (2021).

## 2.3. Fluid flow model

We consider the hypothetical water flooding scenario of a confined reservoir initially fully saturated with oil. We assume the oil is undersaturated, so only two phases (water *w* and oil *o*) are present during simulation. The water is injected at three horizontal injector wells located at one side of the model. Three horizontal producer wells are located at the opposite side of the reservoir model. Each well crosses 70 cells of the grid along the *y*–*axis*. As the exact location of channels varies from one realization to another, such a setting optimizes the chances of perforating permeable channel facies. These wells are equidistant, with the same radius, and perforate the overall width of the reservoir (Fig. 2).

The injectors and the producers are set up with the same flow rate of 3200 m³/day. In the ECLIPSE simulator, when the flow rate cannot be maintained, an alternative condition has to be set: here we choose a minimum bottom-hole pressure of 31.7 MPa and 140 kPa, for the injectors and producers, respectively. These boundary conditions provide a quasi-constant reservoir pressure. The reference pressure is 24.15 MPa at the reservoir top, located at 2.4 km below the surface. Pressure calculations are based on contact positions and hydrostatic equilibrium. In our setting, the viscous forces are one of the main mechanisms for oil displacement and the capillary pressure is neglected. The fluid properties and relative permeability $k_{rγ}$ ($γ = o, w$) curves are given in Table 2 and Fig. 3. Note that the irreducible water saturation and the residual oil saturation are set to zero for simplicity. This ensures that there is no non-mobile oil, and that the recovery only reflects the non-displaced oil phase due to grid-scale reservoir heterogeneities. The production is simulated for a duration of 165 years, so that water breakthrough is reached on all reservoir models.

<table><caption>**Table 2**
Rock and fluid parameters.</caption><thead><tr><th>Reservoir parameters</th><th>Condition</th><th>Values</th><th>Unit</th></tr></thead><tbody><tr><td>Oil density</td><td>Stock-tank condition</td><td>800</td><td>kg/m³</td></tr><tr><td>Oil viscosity</td><td>Reference pressure</td><td>2.85</td><td>MPa s</td></tr><tr><td>Oil formation volume factor</td><td>–</td><td>1.02</td><td>Rm³/Sm³</td></tr><tr><td>Water density</td><td>Stock-tank condition</td><td>1000</td><td>kg/m³</td></tr><tr><td>Water viscosity</td><td>Reference pressure</td><td>0.5</td><td>MPa s</td></tr><tr><td>Reservoir rock compressibility</td><td>Reference pressure</td><td>4e<sup>-10</sup></td><td>Pa<sup>-1</sup></td></tr></tbody></table>

<figure><img src="image_5.png" /><figcaption>**Fig. 3.** Relative permeability curves ($k_{rγ}$-functions). These relative permeabilities correspond to optimal conditions for displacement without leaving any residual oil or water trapped in the porous structure.</figcaption></figure>

## 3. Reservoir metrics

Renard and Allard (2013) emphasize the significant role of connectivity on fluid flow within subsurface reservoirs. To assess and compare turbidite reservoir models described in Section 2.1, we employ static and dynamic metrics. Static metrics, which are computed on the channel facies in the simulation grid, consist of six distinct values that quantify the characteristics of each geostatistical realization (Table 3). As suggested by Jager et al. (2009), the dynamic metrics should contain information about: (i) the production behavior in time, (ii) the behavior over the entire simulation time, and (iii) the behavior over the entire model area. Therefore, four dynamic metrics, derived from the outcomes of flow simulations, are computed for each reservoir model (Table 3).

### 3.1. Static metrics

The static connectivity is depicted from a set of indicators based on the proportion and topology of the facies realizations. These indicators include: the channel facies proportion $φ_c$, the connection probability of channel facies cells $Γ_c$, the number of geobodies $ε_c$, the volume of the largest geobody $V_g$, the static sphericity of the largest geobody $ψ_s$, and the proportion of the channel facies perforated by the wells $δ_c$.

The channel facies proportion $φ_c$ is equal to the ratio between the number of cells of channel facies $n_c$ and the total number of reservoir cells $n_t$:

$$φ_c = \frac{n_c}{n_t} \qquad (1)$$

We compute the number of connected components or geobodies $N_c$, determined by face-connected cells, within the channel facies $c$. From a reservoir modeling perspective, we refer to $N_c$ as the number of geobodies $ε_c$. To compute the connectivity of channel facies $c$, we utilize the connection probability $Γ_c$, which was introduced by Renard
and Allard (2013). This probability measures the proportion of pairs of connected cells among all the pairs of cells within c:

$$ \Gamma_c = \frac{1}{n_c^2} \sum_{i=1}^{N_c} (n_c^i)^2, \qquad (2) $$

where $n_c^i$ the number of cells of the connected component i associated to channel facies c.

From percolation theory, the largest permeable region has implications for transport phenomena and fluid flow in porous media. In our settings, the largest geobody corresponds to the flow backbone and represents the permeable sand region through which water can percolate, influencing the overall flow behavior and connectivity. We compute the volume $V_g$ and surface area $B_g$ of the largest geobody to obtain the static sphericity $\psi_c$. It measures the deviation of geobody reservoir shape from a sphere and is commonly used to describe the compactness of objects (Wadell, 1935). Formally,

$$ \psi_c = \beta \frac{V_g^{2/3}}{B_g}, \qquad (3) $$

where $\beta = (36\pi)^{1/3} \approx 4.836$ is the minimum scaling factor value for a sphere. The static sphericity ($\psi_c$) is equal to 1 for a sphere and $\psi_c \in (0, 1)$ for any other shape. The surface area $B_g$ is estimated by counting the number of voxel faces that discretize the boundary of the object.

The proportion of the channel facies perforated by the wells $\delta_c$ is
$$ \delta_c = \frac{n_{cw}}{n_w} \times 100, \qquad (4) $$

where $n_{cw}$ is the total number of channel facies cells intersected by wells, and $n_w = 420$ is the total number of well cells.

## 3.2. Dynamic metrics

Dynamic metrics are used to characterize the two-phase flow behavior and the production data in each turbidite model. The analysis presented here incorporates four time-series analyzes: the water cut (WC), the oil recovery efficiency (RE), the voidage-replacement ratio (VRR) and the production performance profile.

The WC is defined at any time $t$ as the ratio of the volume of water $V_w$ at stock-tank conditions produced at a well at time $t$ against the volume of the total liquids produced $V_{tot}(t)$:

$$ WC(t) = \frac{V_w(t)}{V_{tot}(t)}. \qquad (5) $$

The oil recovery efficiency (RE) curve, defined at any time $t$ during the simulation, is equal to the total recovered oil (FOPT) at time $t$ divided by the original Oil-In-Place (OIP):

$$ RE(t) = \frac{\text{FOPT}(t)}{\text{OIP}} = \frac{\text{OIP}(t) - \text{OIP}}{\text{OIP}}. \qquad (6) $$

The FOPT is reported in unit volume at which oil is produced in all wells. The oil production is measured in stock tank barrel (stb) or cubic meters (sm³). The RE value is an adimensional number which characterizes the overall efficiency of oil exploitation. This value depends not only on reservoir properties but also reservoir production strategies such as the number of wells, the distance between the wells, and the flow injection rates. From the WC and RE analysis, two metric values are inferred: the water breakthrough time $\omega_c$ that quantifies the earliest travel time of injected water to reach the producers, and the value of oil recovery efficiency $\nu_\omega$ at the water breakthrough time $\omega_c$.

To quantify the effectiveness of fluid injection and displacement processes in the reservoir, we consider the voidage-replacement ratio (VRR) as indicator of waterflood performance at stock-tank conditions (Clark and Ludolph, 2003). The VRR is the ratio of the total volume of fluid injected $V_f$ to the volume of fluids produced $V_{tot}(t)$, and it is defined, for any time $t$, as

$$ VRR(t) = \frac{V_f(t)}{V_{tot}(t)} = \frac{FWIT(t)}{FOPT(t) + FWPT(t)}, \qquad (7) $$

<table><caption>**Table 3**
Reservoir metrics.</caption><thead><tr><th>Metric</th><th>Symbol</th><th>Description</th></tr></thead><tbody><tr><td rowspan="6">Static</td><td>φ<sub>c</sub></td><td>Channel facies proportion</td></tr><tr><td>Γ<sub>c</sub></td><td>Channel facies connection probability</td></tr><tr><td>c<sub>c</sub></td><td>Number of geobodies</td></tr><tr><td>V<sub>g</sub></td><td>Volume of the largest geobody</td></tr><tr><td>ψ<sub>s</sub></td><td>Static sphericity of the largest geobody</td></tr><tr><td>δ<sub>c</sub></td><td>Proportion of channel facies perforated by the wells</td></tr><tr><td rowspan="4">Dynamic</td><td>ω<sub>c</sub></td><td>Water breakthrough time</td></tr><tr><td>ν<sub>ω</sub></td><td>Oil recovery efficiency at ω<sub>c</sub></td></tr><tr><td>α<sub>s</sub></td><td>Surface area of saturation field at 30% PVI</td></tr><tr><td>ψ<sub>d</sub></td><td>Dynamic sphericity of saturation field at 30% PVI</td></tr></tbody></table>

where FWIT is the total volume of water injected, and the FWPT is the volume of water produced. Whereas a VRR of zero indicates solution-gas drive (Tang et al., 2006), recovery processes with a VRR of unity rely solely on viscous forces to displace oil (Clark and Ludolph, 2003). The latter case indicates pressure maintenance, suggesting that the reservoir pressure is maintained.

In conjunction with the VRR, the cumulative production curves of oil and water over the lifetime of the reservoir are considered to interpret the reservoir behavior. This time-series analysis indicates the production performance and effective fluid displacement.

The above metrics satisfy the first two conditions proposed by Jager et al. (2009). To fully capture the dynamics of oil displacement within the reservoir, two additional dynamic metric values are considered. These are based on the saturation front of injected water: the surface area of propagation field $\alpha_s$, and the dynamic sphericity of propagation field $\psi_d$.

The saturation field term $S_\gamma$ is expressed in each grid block as the volume of a fluid $V_\gamma$ divided by the pore volume $V_p$:

$$ S_\gamma(x) = \frac{V_\gamma(x)}{V_p(x)}, \quad \text{with } \gamma = o, w. \qquad (8) $$

The total saturation of fluid is $S_w + S_o = 1$. To summarize the saturation field, we consider the volume $V_{sat}$ where the water saturation is larger than $0.5: V_{sat} = \{x|S_w > 0.5\}$. The boundary $\{B_w = x|S_w(x) = 0.5\}$ of $V_{sat}$ approximates the front of injected water. The shape of water saturation front $S_w$ is characterized with the following metrics:

* the surface area

$$ \alpha_s = |B_w|, \qquad (9) $$

* the dynamic sphericity

$$ \psi_d = \beta \frac{V_{sat}^{2/3}}{|B_w|}, \qquad (10) $$

where $\beta$ is the same geometric correction term as in Eq. (3). The surface area $\alpha_s$ essentially describes the absolute area of the saturation front, whereas the dynamic sphericity $\psi_d$ measures the deviation of water-saturated reservoir shape from a sphere. The saturation field is determined for each reservoir model at a specific pore volume injected (PVI). The PVI term is the ratio between the volume of injected water $V_{inj}$ divided by the total pore volume $V_p$ of the considered reservoir. Instead of computing the saturation fields as a function of time, we calculate the saturation fields as a function of PVI. This way, we compare the saturation field of models under the same injected conditions to compensate for possible differences in net-to-gross between different reservoir models. The PVI value is chosen at 30%, for which no water breakthrough is yet observed on any model.

## 4. Results

The reservoir metrics introduced in Section 3 provide a quantitative characterization to compare the differences across three types of stacking patterns: disorganized (DS), conditioned disorganized (CDS),
<figure><img src="image_6.png" /><figcaption>Fig. 4. Violin and histogram plot representation of static metrics (a–f). DS = disorganized stacking, CDS = conditioned disorganized stacking and OS = organized stacking.</figcaption></figure>

and organized (OS). These metrics enable for a detailed evaluation of
how each stacking set affects reservoir connectivity and reservoir flow
dynamics.

4.1. Static metrics

Fig. 4 illustrates the static metrics for the three datasets using violin plots, with the exception of the number of geobodies, which is presented as discrete values. As mentioned in Section 2.1, all the realizations have approximately the same range of channel facies proportion and the median is equal to 0.123 in each set (Fig. 4a). However, the different number of channels between the sets, and the variable degree of overlap between channels for each setting affect the spread and the shape of the distributions. Nevertheless, all reservoir models have, on average, an equivalent Oil Initially In Place (OIIP).

The connection probability of channel facies illustrates the distri-
bution of channel volumes within the reservoir (Fig. 4b). Disorganized
stacking patterns show a more fragmented channel distribution with
median connectivity of 0.77, indicating the lowest connectivity for
channel facies cells. In contrast, conditioned disorganized stacking
patterns achieve a higher median connectivity of 0.96, suggesting that
nearly all channel facies cells are interconnected.

To assess channel fragmentation, the number of connected geobod-
ies across different stacking models is analyzed in a histogram plot
(Fig. 4c). As all channels go from one side of the domain to the other,
the fragmentation of channels in the disorganized and conditioned
disorganized reservoir models results in multiple preferential pathways
(Fig. 5). DS models typically feature two or three connected geobodies,
while CDS models show one or two, indicating increased connectivity
due to seismic attributes. This results in less variability in channel facies
proportion and connectivity for CDS scenario compared to DS scenario
(Fig. 4a–b). Despite the general trend towards higher connectivity, a
few conditioned disorganized realizations exhibit lower connectivity,
highlighting variability within the set. Organized stacking models, in
contrast, consistently show a single, fully connected channel geobody,
demonstrating the effectiveness of coherent migration algorithms in
enhancing connectivity (Rongier et al., 2017b).

The median volumes for the largest geobody across different stack-
ing sets are 1.20 Gm³ for DS, 1.37 Gm³ for CDS, and 1.39 Gm³ for OS,
indicating a trend towards larger cohesive geobodies in more organized
systems (Fig. 4d). This aligns with the observation that all realizations
have roughly the same volume of channel facies. Fig. 4d also reveals
that, in most cases, the largest geobody constitutes the majority of
<figure><img src="image_7.png" /><figcaption>**Fig. 5.** Example of connected component representation in channel region (vertical exaggeration ×10) for a disorganized (model 32), a conditioned disorganized (model 91) and an organized stacking realization (model 89). The largest connected component is represented in blue.</figcaption></figure>

the reservoir volume, although in some DS and CDS instances, smaller geobodies can hold about half of the sand facies volume.

The static sphericity metric, assessing the shape of the largest geobody, shows a distinct geometric difference between organized and disorganized models (Fig. 4e). Organized stacking models, with a median sphericity of 0.105, demonstrate a more compact shape due to fewer boundary voxel faces as compared to their disorganized counterparts.

The lower sphericity in the DS and CDS sets can be explained by the combination of two factors. First, many crossing channel trajectories yield a decrease in the channel boundary area. Second, the largest geobody often has a smaller vertical extent, hence a lower sphericity than in the OS models (Fig. 5).

The connectivity of channel facies to injector and producer wells distinguishes organized stacking models from the DS and CDS models, with organized models showing better perforation by wells (Fig. 4f). This is due to the OS channel cells being more frequently intersected by wells, particularly in zones of high lateral migration at the bottom of the channel complex.

In summary, the static connectivity metrics indicate the lowest connectivity for the disorganized set. The conditioning by seismic data increases the degree of channel amalgamation, hence the connectivity. The organized set distinctively features a single connected geobody encompassing all channel cells, a compact geobody shape, and enhanced well perforation of the channel facies due to significant lateral migration in the early stages of the deposition process.

## 4.2. Production data analysis

We ran flow simulations on the 300 reservoir models and provide the water-cut (WC), the oil recovery efficiency (RE), the voidage replacement ratio (VRR) and the production profile curves from the start $t_0 = 0$ years to $t_{end} = 165$ years (Fig. 6).

In the DS and CDS cases, many water-cut curves are characterized by several inflection points and some instances where the water saturation locally decreases (Fig. 6a–b). This is due to digitated water saturation front (Fig. 7) reaching the producer wells. Disorganized and conditioned disorganized scenarios exhibit a wider range of water
<figure><img src="image_8.png" /><figcaption>**Fig. 6.** For the three sets of stacking patterns, water cut curves (a-c) and violin plot representation of $\omega_c$ metric at the water breakthrough time (d). In the water cut curves, the black line highlights the representative model presented for each stacking group in Fig. 5.</figcaption></figure>

breakthrough times $\omega_c$, with median values occurring around 70 years and 65 years, respectively (Fig. 6d). The water breakthrough times $\omega_c$ for organized stacking realizations fall within a narrow range (58 years to 65 years) indicating compact unique water flow path with the median value equals to 63.5 years (Fig. 6c-d). These results reflect the impact of channel organization and channel element fragmentation on fluid dynamics (Fig. 6 and Fig. 7).

The maximum oil recovery efficiency ranges from 0.15 to 0.2 (Fig. 8a-c), reflecting the significant of reservoir heterogeneity and the considerable distance between injectors and producers. The three sets show a relatively restricted dispersion in RE curves. The disorganized stacking realizations exhibit the largest uncertainty among the three sets of realizations with a standard deviation of 0.05 at 165 years (Fig. 8). The recovery factors at water breakthrough times also highlight a difference between the organized and the other two stacking models (Fig. 8d). The median value is larger for the disorganized and conditioned disorganized models (0.1) than for organized stacking models (0.085), underscoring that increased structural organization in stacking patterns leads to reduced oil recovery efficiency.

In terms of uncertainty, the variability of the oil recovery is also lowest for the OS realizations. Indeed, flow paths may change laterally from one OS realization to the next, but these changes tend to be mitigated by the use of the horizontal wells because the channels follow the same stacking patterns.

In the voidage replacement ratio (VRR) analysis, an imbalance is noted between $t_0$ and $t = 30$ years for all scenarios, with VRR falling below 0.9, signaling a pressure decline during the early stages of production (Fig. 9a-b). Disorganized models show a more pronounced minimum VRR (0.55 over a median average duration of 3.47 years) compared to conditioned models (0.60 over 3.72 years), indicating less effective early oil displacement by injectors. Organized models experience an even lower VRR (0.49 in 4.59 years), reflecting a significant

pressure imbalance (Fig. 9c). After 50 years, all scenarios stabilize with a VRR near unity, signifying that viscous forces have taken over as the primary oil displacement mechanism. The higher VRR values in disorganized and conditioned disorganized models indicate more effective displacement, as a larger portion of the original fluid is being replaced by the injected fluid. This corresponds to improved sweep efficiency and increased recovery of hydrocarbons, as observed in the efficiency curves (Fig. 8). Consequently, one of the main factors contributing to lower VRR values in the organized models is gravity drainage.

By plotting the cumulative oil production against the cumulative water production (Fig. 9d-f), we observe larger volumes of produced oil without any water production for the DS and CDS cases (on average 230 $Mm^3$ and 220 $Mm^3$, respectively) than for the OS case (on average 180 $Mm^3$). In terms of overall behavior, we also note that the slope of the curves is smallest in the OS case, which confirms the lower sweeping efficiency on the injected water than in the DS and CDS cases (Fig. 9d-f).

### 4.3. Spatial dynamic metrics

This subsection considers the saturation field to gain more insights about the varying production behavior between the three considered modeling methods. The median swept area values at 30% of pore volume injected (PVI) are around 52 $Mm^2$, 49 $Mm^2$ and 33 $Mm^2$ for DS, CDS, and OS models, respectively (Fig. 10a). Spatially, the areas of the water-saturated front at 30% of pore volume injected tend to decrease when the amalgamation between channels increases. An almost opposite trend is observed for the sphericity of the saturation front at 0.3 PVI (Fig. 10b), which is consistent with (but slightly less marked than) the static channel facies sphericity (Section 4.1, Fig. 4e). The near symmetry between the saturation front area and sphericity
<figure><img src="image_9.png" /><figcaption>**Fig. 7.** Water saturation fronts for a disorganized, a conditioned disorganized and an organized stacking models. This visual representation shows only the fluid propagation in the channel regions (vertical exaggeration ×10). The water saturation fields are shown at the initial stage, for 30% of PVI and for 100% of PVI.</figcaption></figure>

can be explained by the fact that the volume of the saturation front is significantly controlled by the amount of injected water. Overall, Fig. 10 confirms that using geostatistical methods that constrain the channel realizations through a sand probability map or morphological rules tends to reduce the complexity of the saturation field and its variability from one realization to another. This is consistent with the conclusions of Larue and Hovadik (2008) which highlights the importance of connectivity and path lengths. Moreover, this suggests that seismically derived facies-proportion cubes and conceptual gradual migration rule-based models reduce the need for generating a large number of channel architectures to capture a representative range

of uncertainties, even at the considered large scale. Accounting for channel avulsion processes and, to a lesser extent, for variations in the horizontal and vertical channel migration rates in the rule-based simulation could, however, have a larger impact on field connectivity.

Fig. 11 illustrates the 0.5 isosurfaces of water saturation fronts for three typical realizations. Their horizontal extension ranges between 9 km and 10 km. In the disorganized and conditioned disorganized scenarios, multiple injected water phase paths are formed. Conversely, the organized scenario showcases a singular geobody, with oil displacement predominantly occurring at the lower region of reservoir due to gravitational forces. This pattern aligns with production data and dynamic connectivity metrics, corresponding to enhanced efficiency in the disorganized scenarios which tend to reduce gravity-driven fluid segregation. A practical approach for improving production performance in the organized scenario could involve closing perforations at the well's lowest point post-water breakthrough, a strategy potentially less effective in disorganized and conditioned disorganized scenarios due to more vertically dispersed water breakthroughs.

### 4.4. Relationship between static and dynamic connectivity metrics

The use of static metrics or fast dynamic proxies is an important avenue to reduce the computational burden of uncertainty quantification (Scheidt and Caers, 2009; Larue and Hovadik, 2006; Bardy et al., 2019; Hird and Dubrule, 1998). In this section, we investigate about the relationships between the various reservoir metrics used in this paper. For each stacking pattern class, realizations with a single connected component are considered to ensure a consistent comparison. The linear correlation coefficient *R* is computed to detect the correlation strength between a static metric and a dynamic metric inferred from bivariate analysis (Fig. 12). As a result, only the key plots are presented here. As the proportion of channel facies directly reflects the volume of this facies in all realizations characterized by a single connected component, cross plots featuring the channel facies proportion are redundant with those indicating the volume of the largest body, hence they are not presented. In all cross-plots, dots represent model values, delineated by both static and dynamic metrics.

The metric designated for gauging the perforation of the channel facies (δc, first column of Fig. 12) does not exhibit a significant correlation with any of the dynamic metrics. This observation justifies a *posteriori* the well placement strategy, by showing that it has no significant impact on the reservoir behavior for the considered scenarios.

The metric quantifying the volume of the largest geobody, $V_g$, demonstrates a notable correlation with the surface area of the saturation field at 30% of pore volume injected ($R = 0.88$ for organized models, $R = 0.64$ for disorganized models, and $R = 0.58$ for conditioned disorganized models, Fig. 12). The high-permeability zones (channel facies here) are indeed a key parameter for fluid circulation. In the OS scenarios, the stacking produces a unique channel facies connected zone (Fig. 4c). For CDS and DS models, there could be several channel geobodies (Fig. 5), but to ensure a consistent comparison, only the realizations with a single component were included in this correlation analysis. As such, in all these cases, $V_g$ is equivalent to the volume of the connected high-permeability zone, which explains its direct impact on the shape of the water front. For the DS and CDS models, this also explains the moderate correlation between $V_g$ and the water breakthrough time ($R = 0.47$ and $R = 0.39$, respectively) and with the oil recovery efficiency ($R = 0.60$ and $R = 0.55$, respectively). Its weak correlation for organized stacking (OS) models may appear less intuitive (Fig. 12: $R = 0.19$ for water breakthrough time $\omega_c$, and $R = 0.14$ for recovery efficiency RE at $\omega_c$). This can be attributed to the gravitational effects discussed earlier, which reduce the sweeping efficiency in OS model realizations (Fig. 7).

The right column of Fig. 12 displays the correlations between static sphericity and the four dynamic metrics. In this study, reduced values of
<figure><img src="image_10.png" /><figcaption>**Fig. 8.** For the three sets of stacking patterns, the oil recovery efficiency (a–c) and violin plot representation of $v_{\omega_c}$ metric at the water breakthrough time $\omega_c$ (d). In the oil recovery efficiency curves, the black line highlights a representative model for each group.</figcaption></figure>

static sphericity correspond to decreased compactness and to the presence of more meandering channel branches, which are characteristic of the DS and CDS models. This leads to a tortuous distribution of injected water fronts and to increased volumes of swept oil. In these scenarios, the static sphericity metric falls short of encapsulating the intricacies of reservoir geometry that govern front propagation, hence the absence of significant statistical correlation between static sphericity and recovery metrics in these two scenarios. Conversely, in organized models, static sphericity shows a negative correlation with oil recovery efficiency. Indeed, a more elongated or intricate permeable channel-shape indicates poorer oil displacement. Due to gravity forces, water invades the bottom of the reservoir, leading to faster displacement with a relatively small surface area. Finally, the interface between the advancing water phase and the remaining oil phase is mainly defined by the interface between channel facies and the inner levees. This is consistent with the observed correlation between static and dynamic sphericity in all models, with a correlation coefficient (R) between 0.55 and 0.60.

In summary, the shape and volume of the largest connected component impact the reservoir behavior. The amalgamation of turbidite channels significantly affects sweeping efficiency in black oil scenarios. The tested static metrics offer only a marginal preview of fluid-flow complexities prior to conducting dynamic simulations. The relationship between static and dynamic metrics is equivalent in disorganized and conditioned disorganized sets. It differs in organized set, because the channel amalgamation type, through gravity effects, impacts the reservoir sweeping process. It is also important to note that this analysis exclusively considers reservoirs constituted by a single connected component; however, in DS and CDS scenarios, other realizations featuring multiple connected components may result in fluid flows in all geobodies. Therefore, there is a need for considering static metrics capable of accommodating these multi-component scenarios as well.

# 5. Discussion and perspectives

## 5.1. Uncertainty quantification and impact of reservoir modeling strategy

Modeling the reservoir behavior of channelized turbidite deposits presents a significant conceptual and computational challenge. Indeed, one must choose among many existing modeling methods which integrate various degrees of conceptual geological knowledge to produce a geocellular reservoir representation at a given scale (Ringrose and Bentley, 2016). Moreover, uncertainty quantification requires, after choosing a modeling method deemed appropriate, to identify key input parameters, their variation ranges, and then to perform multiple simulations to forecast ranges of the dynamic reservoir response (e.g., Larue and Hovadik, 2008). As shown by several studies (e.g., Labourdette et al., 2006; Jager et al., 2009; Alpak et al., 2013; Issautier et al., 2016; Soni et al., 2020; Meirovitz et al., 2021; Alpak and Xue, 2022), the level of realism of the geological heterogeneity representation can have a significant impact on flow. The objective of this paper was to evaluate this impact by comparing the static and dynamic behavior of an object-based Boolean model and a rule-based model on a large-scale reservoir problem including a large number of channels.

In previous studies using the Boolean model, Larue and Hovadik (2008) and Jager et al. (2009) changed channel width, sinuosity and orientation between simulations. Their results showed the limited effect of these geometric parameters, and the important role of the net-to-gross values and of the static connectivity of the reservoir. However, the relative contribution of net-to-gross and connectivity to the flow uncertainty can be challenging to assess from a limited set of simulation samples. In this paper, we have chosen to focus the sensitivity analysis on the connectivity aspects, taking care to use comparable net-to-gross to ensure that static and dynamic metrics were mainly connected to the spatial arrangement of channels.
<figure><img src="image_11.png" /><figcaption>**Fig. 9.** Voidage replacement ratio (VRR) (a–c, left column) and production profile (d–f, right column) for the three sets of stacking patterns. The black line highlights a representative model for each group.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 10.** Violin plot representation of spatial dynamic metrics: (a) the surface area of saturation field at 30% PVI, and (b) the dynamic sphericity metrics. DS = disorganized stacking, CDS = conditioned disorganized stacking and OS = organized stacking.</figcaption></figure>

Another difference concerns the modeling method itself. In Larue and Hovadik (2008), eight realizations are generated by a Boolean approach to reproduce different stacking configurations, from dispersed

single storey channel deposits to more amalgamated channels. Each configuration is represented by one single realization of channel geometry for each stacking trend. Combined with a fixed vertical well
<figure><img src="image_13.png" /><figcaption>Fig. 11. Iso-surfaces for 0.5 water saturation at 30% of PVI for a disorganized, a conditioned disorganized and an organized stacking models (vertical exaggeration x10).</figcaption></figure>

positioning, this confers a crucial role to the channel perforation occurrence, as underlined by the authors. In our settings, we avoid potential bias induced by an arbitrary position of a single vertical well, so our simulation design uses three horizontal injector and producer wells. Whereas we used unconditional simulations, we do not observe significant correlation between dynamic metrics and the proportion of channel facies perforated by wells $\delta_c$ (Fig. 12). This corroborates that our methodology manages to isolate the effect of global connectivity on the results. As a consequence, our results can reliably be interpreted in term of stacking pattern and connectivity variability.

Another difference in this study is the use of state-of-the-art rule-based simulation that produces more realistic stacking patterns than Boolean simulation. For example, the cross-sections in Fig. 8C of Larue and Hovadik (2008) show that the patterns differ in the abundance of channel storeys from base to top, but do not reproduce the progressive evolution that can be generated with more recent migration-based algorithms like the one used in the OS set. The CDS set used in this paper is analogous to the setting of Larue and Hovadik (2008), as it uses a 3D probability cube which constrain the geometry and stacking of independently simulated channels. The DS set is closer to the stationary case used by Jager et al. (2009).

Previous studies using rule-based models integrate much more geological concepts. The modeling methods used by Labourdette et al. (2006) and Alpak et al. (2013) are based on careful outcrop studies, and consider channel complexes including fine scale heterogeneities

thanks to advanced shale drape modeling strategies. Alpak et al. (2013) additionally consider channel stories regrouped in channel story sets. As in the OS set, these works using rule-based simulation methods imply a migration process. As Larue and Hovadik (2008), these studies show the predominant influence of the net-to-gross value, and also of the number of meander belts, and of the shale drape coverage. More recently, Meirovitz et al. (2021) studied the influence of small scale heterogeneities, and they combined their effects with 12 stacking arrangements. The studied geometries correspond to fixed juxtapositions of two successive channels observed on outcrop and extrapolated cylindrically in 3D. The OS set used in our study is simplified as compared to these previous works, as it only considers one scale of genetically related channels, and neglects internal channel heterogeneity. However, it is also more complete in the sense that it considers 100 realizations of channel configurations in a large modeling domain. Even though the net-to-gross variability is kept small and all simulations are confined to a deterministic channel complex, the observed dispersion of static and dynamic metrics highlights the importance of considering a large set of channel configurations in multi-parameter sensitivity studies.

The assumption of homogeneous facies is similar to that of previous studies based on the Boolean model. With the chosen 1/10 permeability anisotropy ratio used in this work, we see that this assumption entails significant gravity-driven fluid segregation. If this assumption is correct, this suggests that the “simplest” Boolean modeling approaches,
<figure><img src="image_14.png" /><figcaption>**Fig. 12.** Bivariate analysis with their correlation coefficients *R* between static and dynamic metrics in all sets (refer to Table 2 for the description of the symbols). Metric values are provided in the supplementary material.</figcaption></figure>

which are currently more widely used, tend to be overly optimistic. This could explain the positive forecasting bias observed in disorganized stacking, as well as in many turbidite reservoir studies (Hempton et al., 2005; Gainski et al., 2010; Meddaugh et al., 2011, 2017).

Nonetheless, the behavior of the OS set should not be transposed directly to areas in which significant basal shale drapes are known to exist, and the present study should then be extended to consider the effect of internal heterogeneities and channel avulsion. Even though this should be confirmed quantitatively, one may qualitatively expect that such features, by increasing the vertical compartmentalization, could make the flow outcomes for the OS case closer to those of the DS and CDS cases, which display less gravity-driven fluid segregation. This would imply that the classical Boolean model with homogeneous facies might, in some cases, compensate vertical compartmentalization (due to internal low permeability facies) by the random independent placement of channels in the reservoir. This can explain why it is possible in some cases to history match reservoir models even though the modeling concept is inappropriate, illustrating the ill-posed nature of the history matching problem.

In principle, increasing the degree of realism of the numerical representation can only help to better manage risks, especially for

spatial problems such as infill well placement, but it often faces computational costs constraints and limitations in practice. We think that more analysis methodologies as in this paper are needed to address this type of questions and better understand the required level of geomodel complexity to appropriately quantify uncertainties. Also, practical sensitivity studies should in principle integrate more static reservoir parameters such as net-to-gross, and also analyze the effect of dynamic parameter uncertainties (e.g., capillary pressure, relative permeability, other reservoir drive mechanisms, reactive transport), but this calls for effective strategies to harness the computational cost (Scheidt et al., 2018).

To summarize, the results presented in this paper are distinct from previous studies because they consider several modeling approaches and a large modeling domain to represent various stacking patterns. We took precautions to isolate the effect of connectivity associated with these methods from other factors known to affect the reservoir behavior in turbiditic channelized reservoirs. This strategy was motivated by computational reasons, as integrating more parameters would call for an impractical number of simulations to ensure statistical representativeness. For the selected settings, we see that different modeling methods
induce ensembles of models with different connectivity and flow behavior distributions. More sensitivity studies, however, are needed to further scrutinize the effect of algorithmic choices on the modeling outcomes as compared to the parameter uncertainty (e.g., for a wider range of net-to-gross values).

## 5.2. Metrics and simulation settings

In an uncertainty quantification study, objective analyzes and comparisons require to use adapted metrics. In this study, we followed the recommendations of Jager et al. (2009) and used four time-series: the water cut (WC), the oil recovery efficiency (RE), the voidage replacement ratio (VRR) and the production performance profile. From these curves, we extracted the water breakthrough time $\omega_c$, and the value of oil recovery efficiency $v_\omega$ at $\omega_c$. All these curves and metrics are classical tools regularly used in comparison studies (e.g., Labourdette et al., 2006; Alpak et al., 2013; Alpak and Barton, 2014; Meirovitz et al., 2021).

Nonetheless, they do not allow to understand the full dynamics of oil displacement within the reservoir, especially the spatial evolution of fluid phases. To do so, we introduced simple metrics to describe the shapes of the observed saturation fronts at a specific time, chosen as 30% of Pore Volume Injected (PVI) in this study. These new metrics allow us to objectify the visual observations made on three-dimensional saturation fields, and highlight the specific reservoir behavior of organized stacking (OS) models as compared to the DS and CDS models produced by the Boolean method (Fig. 10).

The static metrics used here had two goals: (i) to separate connectivity and shape from other factors (see discussion in Section 5.1 on the NTG and the channel facies perforation by the wells $\delta_c$), and (ii) to potentially find proxies of reservoir dynamic behavior.

The latter objective is essential to reduce the computational cost of practical sensitivity studies (Scheidt et al., 2018). For this, we defined four static metrics, which were expected to be potentially discriminant between modeling approaches and flow behaviors. The number of connected geobodies being a discrete value, it was not directly used in a simple Pearson correlation analysis even though it helps to understand the processes. The largest geobody volume $V_g$ appears to be the most effective metric of our test. This is not surprising as it directly refers to the connected high-permeability volume, which has been proven from percolation theory to have a crucial role. As we analyzed the realizations consisting of only one channel geobody for consistency, and as this geobody is connected to the well doublet, $V_g$ becomes in our case an equivalent measurement of the static reservoir connectivity as defined by Larue and Hovadik (2006, 2008). However, this metric is clearly insufficient to characterize the late dynamic behavior of the field, which is affected by smaller connected bodies. The two other metrics $\psi_s$ and $\Gamma_c$ show only weak correlations with dynamic metrics. The channel facies connection probability $\Gamma_c$ is probably inadequate in our case, because the focus on the largest connected component makes it difficult to compare between the organized set and the disorganized sets. If the static sphericity $\psi_s$ is moderately correlated to the shape of the saturation front, the reservoir performance seems to result from several interlaced factors that cannot be approached by such a simple geometric measurement. When net-to-gross is a fixed “input” parameter and well location is optimized to fit several channel configurations, our study suggests that it becomes difficult to forecast the dynamic reservoir behavior from simple static metrics. Note, however, that the small number of connected reservoir components in all realizations may explain this lack of correlation. By increasing the number of components, experiments with lower net-to-gross values and rule-based models considering avulsion processes could possibly highlight better relationships between static and dynamic metrics.

Nonetheless, we observe distinct statistical relationships between static and dynamic metrics for the three considered modeling approaches (Fig. 12): introducing more constraints through probability

maps or morphological rules using a migration algorithm, results in less variability (both static and dynamic) than using the common assumption of independent channels in a stationary domain. In the organized set, this reduction in variability can be partially attributed to the use of fixed migration and aggradation parameters for all realizations. Varying these rates between realizations would likely increase the diversity of the organized set. It is probable that perturbing these rates from one realization to the next one would increase the diversity of the OS set, but probably not to the same level as in the DS case, as migration-based approaches are bound to produce a higher degree of connectivity than classical Boolean models. As discussed in Section 5.1, these observations on metric relationships should be taken with caution, as they could be affected by the representation of internal heterogeneities in the channels. This aspect needs to be further investigated in future studies.

## 6. Conclusions

In sparse data settings, channel envelopes cannot be directly inferred from field observations, and various geostatistical methods are developed and used instead to generate their geometry. The results presented in Section 4 show that, for comparable net-to-gross values and no internal heterogeneities, different channel simulation methods yield different flow responses, different ranges of uncertainty, and different relationships between static and dynamic metrics.

In organized models, where the channel migration process is simulated, early water breakthrough and lower recovery efficiency suggests inefficient oil displacement. This unexpected behavior, despite simpler shapes and better connectivity of the channels, can be attributed to gravity-driven fluid segregation during reservoir production. The new metrics characterizing the shape of the saturation front provide a way to quantify this effect. The disorganized stacking patterns as obtained by classical software packages that implement the Boolean Model and, to a smaller extent, conditioned disorganized stacking patterns, exhibit delayed water breakthrough times and more optimistic recovery efficiencies. In terms of uncertainty assessment, the variability among the different stochastic channel architectures tends to decrease as seismic-based proportions and process-based gradual migration concepts are considered in the modeling method. However, in complement to previous work (Labourdette and Bez, 2010; Alpak et al., 2013) further studies should be made at global scale to assess under what circumstances more geological realism (e.g., internal facies heterogeneities, integration of channel avulsion in the rule-based simulation) affects reservoir forecasts in terms of average behavior and uncertainty range. Checking whether our results hold for various net-to-gross scenarios and considering internal petrophysical heterogeneities, appears as an obvious but computationally challenging way forward.

Concerning static and dynamic metrics, no clear correlation was found independently of the modeling method, but relationships appear to be distinct for the three considered modeling methods. This mitigates the generalization of previous results which often fixed the channel geometry modeling strategy, and underscores the challenges in forecasting dynamic behavior based on static connectivity metrics independently of the net-to-gross values.

Overall, our results highlight the role of modeling methods both in average and range of reservoir behavior, even for a narrow range of net-to-gross values. This illustrates the impact of the integration of data (seismic-based proportion) and of geological concepts (channel migration) in reservoir modeling methods.

### CRediT authorship contribution statement

**Enrico Scarpa:** Writing – review & editing, Writing – original draft, Software, Data curation, Conceptualization. **Pauline Collon:** Writing – review & editing, Validation, Supervision, Conceptualization. **Irina Panfilova:** Writing – review & editing, Validation, Software. **Guillaume Caumon:** Writing – review & editing, Validation, Supervision, Conceptualization.
Declaration of competing interest

The authors declare that they have no known competing finan-
cial interests or personal relationships that could have appeared to
influence the work reported in this paper.

Acknowledgments

This work was performed in the frame of the RING project at Université de Lorraine. We would like to thank for their support the industrial and academic sponsors of the RING-Gocad Consortium managed by ASGA (www.ring-team.org/Consortium). We also acknowledge AspenTech for the SKUA-Gocad software and SLB for the ECLIPSE software. The simulations are performed with the ECLIPSE100 software.

Appendix A. Supplementary data

Supplementary material related to this article can be found online
at https://doi.org/10.1016/j.marpetgeo.2024.107235.

Data availability

Data will be made available on request.

References

Alpak, F.O., Barton, M.D., 2014. Dynamic impact and flow-based upscaling of the estuarine point-bar stratigraphic architecture. J. Pet. Sci. Eng. 120, 18-38.

Alpak, F.O., Barton, M.D., Naruk, S.J., 2013. The impact of fine-scale turbidite channel architecture on deep-water reservoir performance. AAPG Bull. 97, 251-284.

Alpak, F.O., Xue, G., 2022. Effects of fine-scale turbidite lobe stratigraphic architecture on dynamic reservoir performance. Mar. Pet. Geol. 139, 105540.

Bardy, G., Biver, P., Caumon, G., Renard, P., 2019. Oil production uncertainty assessment by predicting reservoir production curves and confidence intervals from arbitrary proxy responses. J. Pet. Sci. Eng. 176, 116-125.

Bertoncello, A., Sun, T., Li, H., Mariethoz, G., Caers, J., 2013. Conditioning surface-based geological models to well and thickness data. Math. Geosci. 45, 873-893.

Clark, R., Ludolph, B., 2003. Voidage replacement ratio calculations in retrograde condensate to volatile oil reservoirs undergoing EOR processes. In: SPE Annual Technical Conference and Exhibition. SPE, SPE-84359.

Cojan, I., Fouché, O., Lopéz, S., Rivoirard, J., 2005. Process-based reservoir modelling in the example of meandering channel. In: Geostatistics Banff 2004. Springer, pp. 611-619.

Covault, Jacob A., Sylvester, Zoltán, Ceyhan, Can, Dunlap, Dallas B., 2021. Giant meandering channel evolution, Campos deep-water salt basin, Brazil. In: Geosphere 17. Geological Society of America, pp. 1869-1889.

Deptuck, M.E., Sylvester, Z., Pirmez, C., O'Byrne, C., 2007. Migration-aggradation history and 3-D seismic geomorphology of submarine channels in the pleistocene Benin-major Canyon, Western Niger Delta slope. Mar. Pet. Geol. 24, 406-433.

Deutsch, C., Tran, T., 2002. Fluvsim: a program for object-based stochastic modeling of fluvial depositional systems. Comput. Geosci. 28, 525-535.

Deutsch, C.V., Wang, L., 1996. Hierarchical object-based stochastic modeling of fluvial reservoirs. Math. Geol. 28, 857-880.

Gainski, M., Macgregor, A., Freeman, P., Nieuwland, H., 2010. Turbidite Reservoir Compartmentalization and Well Targeting with 4d Seismic and Production Data: Schiehallion Field, Uk. vol. 347, Geological Society, London, Special Publications, pp. 89-102.

Haldorsen, H.H., Damsleth, E., 1990. Stochastic modeling (includes associated papers 21255 and 21299). J. Pet. Technol. 42, 404-412.

Hempton, M., Marshall, J., Sadler, S., Hogg, N., Charles, R., Harvey, C., 2005. Turbidite reservoirs of the sele formation, central north sea: geological challenges for improving production. In: Geological Society, London, Petroleum Geology Conference Series. The Geological Society of London, pp. 449-459.

Hird, K., Dubrule, O., 1998. Quantification of reservoir connectivity for reservoir description applications. SPE Reserv. Eval. Eng. 1, 12-17.

Holden, L., Hauge, R., Skare, O., Skorstad, A., 1998. Modeling of fluvial reservoirs with object models. Math. Geol. 30, 473-496.

Hovadik, J.M., Larue, D.K., 2007. Static characterizations of reservoirs: refining the concepts of connectivity and continuity. Petrol. Geosci. 13, 195-211.

Issautier, B., Viseur, S., Audigane, P., Chiaberge, C., Le Nindre, Y.M., 2016. A new approach for evaluating the impact of fluvial type heterogeneity in co₂ storage reservoir modeling. C. R. Geosci. 348, 531-539.

Jager, G.de., Van Doren, J.F., Jansen, J.D., Luthi, S.M., 2009. An evaluation of relevant geological parameters for predicting the flow behaviour of channelized reservoirs. Petrol. Geosci. 15, 345-354.

Karssenberg, D.E., Bridge, J.S., 2008. A three-dimensional numerical model of sediment transport, erosion and deposition within a network of channel belts, floodplain and hill slope : extrinsic and intrinsic controls on floodplain dynamics and alluvial architecture. Sedimentology 55, 1717-1745.

Karssenberg, D., Törnqvist, J.S., 2001. Conditioning a process-based model of sedimentary architecture to well data. J. Sediment. Res. 71, 868-879.

Labourdette, R., 2007. Integrated three-dimensional modeling approach of stacked turbidite channels. AAPG Bull. 91, 1603-1618.

Labourdette, R., 2008. 'LOSCS' lateral offset stacked channel simulations: Towards geometrical modelling of turbidite elementary channels. Basin Res. 20, 431-444.

Labourdette, R., Bez, M., 2010. Element migration in turbidite systems: Random or systematic depositional processes. AAPG Bull. 94, 345-368.

Labourdette, R., Poncet, J., Seguin, J., Temple, F., Hegre, J., Irving, A., 2006. Three-dimensional modelling of stacked turbidite channels in west africa: impact on dynamic reservoir simulations. Petrol. Geosci. 12, 335-345.

Larue, D.K., Hovadik, J., 2006. Connectivity of channelized reservoirs: a modelling approach. Petrol. Geosci. 12, 291-308.

Larue, D., Hovadik, J., 2008. Why is reservoir architecture an insignificant uncertainty in many appraisal and development studies of clastic channelized reservoirs? J. Petrol. Geol. 31, 337-366.

Lemay, M., Grimaud, J.L., Cojan, I., Rivoirard, J., Ors, F., 2020. Geomorphic variability of submarine channelized systems along continental margins: Comparison with fluvial meandering channels. Mar. Pet. Geol. 115, 104295.

Li, W., Li, S., Yu, J., Shi, J., 2023. A well conditioning method of object-based models. Interpretation 11, SA115-SA125.

Lopez, S., 2003. Modélisation de réservoirs chenalisés méandriformes: approche génétique et stochastique. (Ph.D. thesis). Ecole des Mines de Paris.

Mayall, M., Jones, E., Casey, M., 2006. Turbidite channel reservoirs—key elements in facies prediction and effective development. Mar. Pet. Geol. 23, 821-841.

McHargue, T., Pyrcz, M.J., Sullivan, M.D., Clark, J., Fildani, A., Romans, B., Covault, J., Levy, M., Posamentier, H., Drinkwater, N., 2011. Architecture of turbidite channel systems on the continental slope: patterns and predictions. Mar. Pet. Geol. 28, 728-743.

Meddaugh, W.S., Champenoy, N., Osterloh, W.T., Tang, H., 2011. Reservoir forecast optimism-impact of geostatistics, reservoir modeling, heterogeneity, and uncertainty. In: SPE Annual Technical Conference and Exhibition. SPE.

Meddaugh, W., Meddaugh, W., McCray, B., 2017. Quantitative assessment of the impact of sparse data and decision bias on reservoir recovery forecasts. In: SPE Annual Technical Conference and Exhibition. OnePetro.

Meirovitz, C.D., Straight, L., Hubbard, S.M., Romans, B.W., 2021. The influence of inter-and intra-channel architecture on deep-water turbidite reservoir performance. Petrol. Geosci. 27, petgeo2020-005.

Morris, Paul, Sylvester, Zoltán, Jacob, A., Mohrig, David, 2022. Channel trajectories control deep-water stratigraphic architecture. In: The Depositional Record. Wiley Online Library.

Parquer, M., Collon, P., Caumon, G., 2017. Reconstruction of channelized systems through a conditioned reverse migration method. Math. Geosci. 49, 965-994.

Parquer, M., Yan, N., Colombera, L., Mountney, N., Collon, P., Caumon, G., 2020. Combined inverse and forward numerical modelling for reconstruction of channel evolution and facies distributions in fluvial meander-belt deposits. Mar. Pet. Geol. 104409.

Peakall, J., McCaffrey, B., Kneller, B., 2000. A process model for the evolution, morphology, and architecture of sinuous submarine channels. J. Sediment. Res. 70, 434-448.

Pyrcz, M.J., Deutsch, C.V., 2005. Conditioning event-based fluvial models. In: Geostatistics Banff 2004. Springer, pp. 135-144.

Pyrcz, M.J., Sech, R.P., Covault, J.A., Willis, B.J., Sylvester, Z., Sun, T., 2015. Stratigraphic rule-based reservoir modeling. Bull. Canad. Petrol. Geol. 63, 287-303.

Renard, P., Allard, D., 2013. Connectivity metrics for subsurface flow and transport. Adv. Water Resour. 51, 168-196.

Ringrose, P., Bentley, M., 2016. Reservoir Model Design. Springer.

Rongier, G., Collon, P., Renard, P., 2017a. A geostatistical approach to the simulation of stacked channels. Mar. Pet. Geol. 82, 318-335.

Rongier, G., Collon, P., Renard, P., 2017b. Stochastic simulation of channelized sedimentary bodies using a constrained L-system. Comput. Geosci. 105, 158-168.

Rongier, G., Collon, P., Renard, P., Straubhaar, J., Sausse, J., 2016. Comparing connected structures in ensemble of random fields. Adv. Water Resour. 96, 145-169.

Ruiu, J., Caumon, G., Viseur, S., 2016. Modeling channel forms and related sedimentary objects using a boundary representation based on Non-Uniform Rational B-Splines. Math. Geosci. 48, 259-284.

Scheidt, C., Caers, J., 2009. Representing spatial uncertainty using distances and kernels. Math. Geosci. 41, 397-419.

Scheidt, C., Li, L., Caers, J., 2018. Quantifying Uncertainty in Subsurface Systems. Vol. 236, John Wiley & Sons.

Shmaryan, L., Deutsch, C., 1999. Object-based modeling of fluvial/deepwater reservoirs with fast data conditioning: methodology and case studies. In: SPE Annual Technical Conference and Exhibition. SPE, SPE-56821.
Soni, K., Manzocchi, T., Haughton, P., Carneiro, M., 2020. Hierarchical characterization and modelling of deep-water slope channel reservoirs. In: SPE Norway Subsurface Conference. SPE, D021S005R005.

Stewart, J., Dunn, P.A., Lyttle, C., Campion, K., Oyerinde, A., Fischer, B., 2008. Improving performance prediction in deep-water reservoirs: Learning from outcrop analogues, conceptual models and flow simulation. In: International Petroleum Technology Conference. IPTC-12892.

Sylvester, Z., Pirmez, C., Cantelli, A., 2011. A model of submarine channel-levee evolution based on channel trajectories: Implications for stratigraphic architecture. Mar. Pet. Geol. 28, 716-727.

Talling, P.J., Allin, J., Armitage, D.A., Arnott, R.W., Cartigny, M.J., Clare, M.A., Felletti, F., Covault, J.A., Girardclos, S., Hansen, E., et al., 2015. Key future directions for research on turbidity currents and their deposits. J. Sediment. Res. 85, 153-169.

Tang, G.Q., Sahni, A., Gadelle, F., Kumar, M., Kovscek, A.R., 2006. Heavy-oil solution gas drive in consolidated and unconsolidated rock. SPE J. 11, 259-268.

Teles, V., de Marsily, G., Perrier, E., 1998. Sur une nouvelle approche de modélisation de la mise en place des sédiments dans une plaine alluviale pour en représenter l'hétérogénéité. Comptes Rendus l'Académie des Sci. - Ser. IIA - Earth Planet. Sci. 327, 597-606.

Viseur, S., 2004. Turbidite reservoir characterization: object-based stochastic simulation meandering channels. Bulletin de la Société Géologique de France 175, 11-20.

Wadell, H., 1935. Volume, shape, and roundness of quartz particles. J. Geol. 43, 250-280.

Wang, Y.C., Pyrcz, M.J., Catuneanu, O., Boisvert, J.B., 2018. Conditioning 3d object-based models to dense well data. Comput. Geosci. 115, 1-11.

Wen, R., 2005. 3D geologic modelling of channellized reservoirs: applications in seismic attribute facies classification. First Break 23.

Wietzerbin, L., Mallet, J., 1994. Parameterization of complex 3d heterogeneities: a new cad approach. SPE Comput. Appl. 6, 11-18.

Wu, X.H., Parashkevov, R., 2009. Effect of grid deviation on flow solutions. SPE J. 14, 67-77.