Impact of sedimentological heterogeneity on subsurface storage in net-transgressive, shallow-marine sandstones: Cliff House Sandstone, New Mexico, USA Open Access

The Cliff House Sandstone records approximately 170 km of shoreline retreat over a period of c. 4 Myr (Fig. 1; Molenaar 1983). The unit consists of a retrogradationally stacked succession of marginal-marine and shallow-marine deposits, recording the development of regressive strandplain and wave-dominated deltaic shorelines that alternated with transgressive barrier islands (Sears et al. 1941; Palmer and Scott 1984; Donselaar 1989; Olsen et al. 1999; Jordan et al. 2016). It represents deposition along a NW–SE-trending stretch of shoreline that fringed part of a broad embayment (‘Utah Bight’) along the western margin of the Cretaceous Western Interior Seaway (Fig. 1; Kauffman and Caldwell 1993; Van Cappelle et al. 2018). Shoreline orientation and patterns of shoreline retreat reflect a combination of tectonic subsidence, sea-level fluctuations, and sediment supply from the Sevier fold-and-thrust belt to the west and Mogollon Highlands to the SW (Fig. 1; e.g. Krystinik and DeJarnett 1995; Lawton et al. 2003; Painter and Carrapa 2013; Szwarc et al. 2015; Van Cappelle et al. 2018; Li and Aschoff 2022). The combination of these controls is reflected in the locally variable thickness (up to 240 m) and stratigraphic-architectural style of the Cliff House Sandstone (e.g. Fig. 1) and of net-transgressive sandstone successions more generally. Thick successions comprising vertically stacked, regressive–transgressive sandstone tongues (parasequences) separated locally by mudstones are developed where the net-transgressive shoreline trajectory (i.e. overall angle of shoreline retreat) was relatively steep (e.g. Jordan et al. 2016; cf. Cattaneo and Steel 2003; Sixsmith et al. 2008). In contrast, thin successions of erosionally amalgamated sandstones in which mudstones are poorly preserved characterize relatively shallow net-transgressive shoreline trajectories (cf. Cattaneo and Steel 2003; Sixsmith et al. 2008).

The upper part of the Cliff House Sandstone, which is Middle Campanian in age, is exposed in deeply incised canyons in Chaco Culture National Historical Park, NW New Mexico, USA, resulting in 3D outcrop control (Scott et al. 1984; Donselaar 1989; Jordan et al. 2016) (Fig. 2). Here, the Cliff House Sandstone contains 11 parasequences of shoreface-shelf sandstones (wave-dominated facies association in Table 1) that are aggradationally-to-retrogradationally stacked into four parasequence sets (Jordan et al. 2016). Each shoreface-shelf sandstone parasequence is laterally extensive along the (NW–SE) depositional strike of the palaeoshoreline, interfingers with and grades into offshore mudstones down depositional dip (towards NE), and pinches out into mudstone-dominated coastal plain deposits and/or channelized tidal sandbodies up depositional dip (towards SW) (continental and tide-dominated facies associations in Table 1; Fig. 2; Jordan et al. 2016; cf. Sixsmith et al. 2008). Internally, shoreface-shelf sandstone parasequences contain a shallowing-upward facies succession that comprises, from base to top: offshore mudstones, distal lower shoreface sandstone and mudstone heteroliths, proximal lower shoreface sandstones, and top-truncated upper shoreface and foreshore sandstones (Table 1; Fig. 2; Jordan et al. 2016). Each parasequence is bounded at its base and top by laterally extensive, wave ravinement surfaces that are interpreted to have been cut by wave erosion during transgression (Fig. 2; Jordan et al. 2016; cf. Swift 1968; Sixsmith et al. 2008). Channelized tidal sandbodies (facies T1 in Table 1) erode into, and define the updip termination of, most shoreface-shelf sandstone tongues. In the latter examples, the basal erosion surfaces of the channelized tidal sandbodies are interpreted as tidal ravinement surfaces (Fig. 2; Jordan et al. 2016; cf. Swift 1968; Sixsmith et al. 2008). Tidally influenced fluvial channel-fill sandbodies (facies T2 in Table 1) are located below and palaeolandward of the lower three shoreface-shelf sandstone tongues (Fig. 2; Jordan et al. 2016).

We use three correlation panels along major, WSW–ENE-oriented canyons (Fig. 2), oblique to depositional dip, to construct a suite of reservoir models (Fig. 3). The three correlation panels are part of the dataset documented in Jordan et al. (2016), and share a stratigraphic framework that is also based on a WNW–ESE-oriented correlation panel, oblique to depositional dip. The WSW–ENE-oriented correlation panels are treated as three parallel planes (Fig. 2), to facilitate their use in sketch-based model construction. There is little uncertainty in the distribution and internal architecture of shoreface-shelf sandstone tongues between the three correlation panels, because the sandstone tongues are laterally extensive along depositional strike (wave-dominated facies association in Table 1). In contrast, there is uncertainty in the geometry and distribution of channelized tidal and tidally influenced fluvial sandbodies (facies T1 and T2 in Table 1) between the three panels.

In order to assess the influence of sedimentological heterogeneity on subsurface fluid migration and storage in the Cliff House Sandstone example, we selected six heterogeneities for investigation (Table 2). Three of these heterogeneities reflect the geometry, distribution and spatial organization of facies units (Table 1), while three further heterogeneities reflect the values of porosity and permeability assigned to selected facies. The six heterogeneities are selected to assess the impact of uncertainties in a typical subsurface dataset and/or the outcrop dataset, as outlined below.

1. The plan-view distribution and geometry of tidally influenced fluvial channel-fill sandbodies (facies T2) is uncertain in regions between the three correlation panels along major, WSW–ENE-oriented canyons (Fig. 2). We develop three different interpreted scenarios of the orientation, width and plan-view distribution of these channelized sandbodies, each of which is consistent with palaeocurrent data in figure 6C of Jordan et al. (2016) (Table 1). Palaeocurrents exhibit a wide spread of orientations, with mean palaeoflow orientations towards the SW (N237°) and north (N006°) in lower and upper units of tidally influenced fluvial channel-fill sandbodies (figure 6C of Jordan et al. 2016). Channelized sandbodies with similar vertical positions relative to a widespread datum (wave ravinment surface wRS 600 in Fig. 2) are correlated between the three panels. Scenario 1 contains three channelized sandbodies that represent channel belts: the oldest is oriented east–west (Fig. 4a), overlain by a second north–south-oriented channelized sandbody and a third SW–NE-oriented channelized sandbody (Fig. 4b). The youngest sandbody does not extend to the northwestern side of the study area, where it is eroded by overlying tidal channels. Scenario 2 also contains three channelized sandbodies representing channel belts. The oldest sandbody is oriented SW–NE, the second sandbody is sinuous with a mean NW–SE orientation (Fig. 4c), while the youngest sandbody is similar to that in Scenario 1 (Fig. 4d). Scenario 3 contains four channel-belt sandbodies; the first, second and third sandbodies are oriented NW–SE and stacked with lateral offsets towards the SW and NE (Fig. 4e, f), while the youngest sandbody is similar to that in Scenario 1 (Fig. 4f). Channel-belt sandbodies in each scenario are 5–7 km wide with width:thickness ratios of 5–30, which are typical of deltaic distributaries (Gibling 2006).

2. Two interpretations are developed for the thickness and vertical connectivity of tidal (facies T1) and tidally influenced fluvial (facies T2) channel-fill sandbodies in between the three cliff-face correlation panels (Table 1). The channel-belt sandbodies in the three scenarios described above are either stacked vertically to form a larger multistorey body (Fig. 4), or are separated vertically by intervals of mudstone-dominated coastal plain deposits. In combination with their plan-view distribution, the thickness and vertical connectivity of channelized sandbodies are expected to influence the proportion of tidal and tidally influenced fluvial channel-fill sandstone facies T1 and T2, and, by corollary, the proportion of non-channelized continental facies Cp1 and T3. The connectivity of channelized sandbodies has also been shown to control fluid flow patterns and recovery in hydrocarbon reservoirs (Jones et al. 1995; Larue and Hovadik 2006; Villamizar et al. 2015).

3. Mouth bar deposits are treated either as distinct, localized facies bodies (facies D1) developed at the down-dip terminations of tidally influenced fluvial channel-fill sandstones (facies T2) or are grouped together with proximal lower shoreface facies belts (facies W3) (Table 1). These two contrasting treatments are considered to reflect uncertainty in the penetration and subsequent identification of localized mouth bar deposits by wells in the subsurface.

4. Coastal plain (facies Cp1) and tidal flat (facies T3) deposits are relatively poorly exposed in the Cliff House Sandstone outcrops, and the two facies are grouped together in this study. Both facies are dominated by mudstones, but contain subordinate sandstones and siltstones (Table 1). In hydrocarbon reservoir studies, the two facies have been considered either as non-reservoir deposits (ø = 10%, k_h = k_v = 0.001 mD) (e.g. Livera 1989) or as poor-quality reservoir deposits (ø = 20%, k_h = 20 mD, k_v = 0.02 mD) (e.g. Abbots and Van Kuijk 1997; Colombera and Mountney 2021) (Table 2).

5. Tidally influenced fluvial channel-fill sandbodies (facies T2) are comparable in their grain-size and textural characteristics as tidal channel-fill sandbodies (facies T1) (Table 1). We consider two interpretations, one in which the two facies are assigned the same reservoir quality (ø = 27%, k_h = 2000 mD, k_v = 200 mD) and one in which tidally influenced fluvial channel-fill sandbodies (facies T2) are assigned higher reservoir properties (ø = 28%, k_h = k_v = 4000 mD) (Table 2). Such differences in reservoir quality between channelized sandbodies have been demonstrated to influence recovery in hydrocarbon reservoirs (e.g. Jones et al. 1995).

6. Distal lower shoreface (facies W4) and offshore (facies W5) deposits are also grouped together in this study. The former comprise interbedded fine-grained sandstones, siltstones and mudstones, while the latter is dominated by mudstones (Table 1). Both facies contain sandstone beds with steep-sided erosional scours or gutter casts at their base (Jordan et al. 2016), which can enhance sandstone connectivity and effective permeability (Onyenanu et al. 2019). The two facies are therefore considered either as non-reservoir deposits (ø = 10%, k_h = k_v = 0.001 mD) or as poor-quality reservoir deposits (ø = 20%, k_h = 20 mD, k_v = 0.02 mD) (Table 2).

We do not vary the porosity and permeability of proximal lower shoreface (facies W3), upper shoreface (facies W2), foreshore (facies W1) and tidal channel-fill (facies T1) deposits (Tables 1 and 2). Therefore, the impact of variations in these properties on reservoir character and behaviour is not investigated in our study. Additional heterogeneities at small (centimetre-to-decametre) scales, which are important in generating differences in relative permeability and capillary entry pressure during multiphase flow (e.g. Saadatpoor et al. 2009; Gershenzon et al. 2016; Jackson and Krevor 2020), can aid trapping of CO₂ and H₂, for example. Consideration of these heterogeneities lies beyond the scope of this work, but is the focus of ongoing work.

Twelve models have been constructed of the Cliff House Sandstone (Fig. 3), representing all possible combinations of the interpreted scenarios for heterogeneities at large length scales that impact the geometry, distribution and spatial organization of facies bodies (labelled 1–12 in Table 3). Each of the 12 models has eight sets of porosity and permeability values assigned to it, representing all possible combinations of the interpreted scenarios for heterogeneities at relatively small length scales that are represented implicitly by the values of porosity and permeability assigned to facies Cp1/T3, T2 and W4/W5 (labelled A–H in Table 4). Thus, a suite of 96 models and porosity-permeability sets (labelled 1A–H through to 12A–H) were used (Tables 3 and 4). This suite constitutes a full factorial experimental design (Box et al. 1978).

The 12 models were constructed using a sketch-based approach that allows geological concepts and scenarios to be rapidly and intuitively captured by non-experts, and which is implemented in the Open Source research code, RRM (Costa Sousa et al. 2020; Jacquemyn et al. 2021). Geological architectures and heterogeneities (e.g. sequence stratigraphic surfaces, facies boundaries) are represented by surfaces (cf. Denver and Phillips 1990; Hamilton and Jones 1992) that define and bound geological domains (cf. Pyrcz et al. 2005; Caumon et al. 2009; Sech et al. 2009; Ruiu et al. 2016; Jacquemyn et al. 2019). Surfaces and surface-bounded geological domains are generated and manipulated using Sketch-Based Interface and Modelling (SBIM) methods that were developed for non-geological CAD and CFD applications. Operators that necessitate geological viability control interactions between SBIM-generated surfaces.

Our sketch-based modelling approach allows surfaces to be sketched in any order (Jacquemyn et al. 2021), although we follow the approach taken by Jackson et al. (2022) and Alshakri et al. (2023) in constructing comparable suites of sketch-based models. To ensure consistent and rapid model construction, surfaces that are shared across all 12 models were sketched initially and then reused. Surfaces that define heterogeneities that impact the geometry and distribution of facies units in different models (Table 3) are sketched later. Laterally extensive stratigraphic surfaces of relatively simple 3D geometry, including wave ravinement and tidal ravinement surfaces, were sketched in and linearly extruded between each cliff-face panel (Costa Sousa et al. 2020; Jacquemyn et al. 2021). More geometrically complex and less extensive surfaces, including the erosional bases of channelized, tidally influenced fluvial sandbodies (facies T2; Fig. 4), are sketched in one cliff-face panel and their cross-sectional geometries then extruded along a plan-view trajectory (Costa Sousa et al. 2020; Jacquemyn et al. 2021).

Each model has dimensions of 10 770 m (NE–SW) × 5130 m (NW–SE) × 102 m (thickness), and is visualized with a vertical exaggeration of ×35 (Fig. 3). Since our focus is on investigating the effects of sedimentological heterogeneity, structural elements such as faults and tectonic dip are not considered in the models. Models are generated without reference to an underlying grid. Once the model has been generated, a grid is then created to visualize the models (e.g. Fig. 3; Jacquemyn et al. 2021) or to perform numerical calculations (Petrovskyy et al. 2023). In the latter context, grid-cell size and grid resolution are discussed below.

Visual inspection of the sketch-based models confirms consistency between the modelled stratigraphic architectures (Figs 3, 4) and the underlying geological data and interpretations (Fig. 2). The modelled architectures thus represent the combinations and settings of heterogeneities outlined in Table 3.

Calculations of pore volume are directly obtained from the sketched models after assigning porosity and permeability values to each facies (Table 1). Our implementation of sketch-based reservoir modelling is integrated with computationally efficient flow diagnostics, which allow key flow properties and behaviours to be assessed rapidly using controlled numerical experiments that rely on a reduced-physics, single-phase pressure solution (Shahvali et al. 2012; Møyner et al. 2014; Rasmussen and Lie 2014; Lie et al. 2015). The fundamental governing equations in flow diagnostics are based on Darcy's law and assume incompressible single-phase flow, absence of gravity and mass conservation (Møyner et al. 2014; Rasmussen and Lie 2014; Lie et al. 2015); Petrovskyy et al. (2023) provide detailed documentation of how flow diagnostics are implemented in the RRM code. Flow diagnostic outputs include the time of flight and stationary distribution of tracer fluid obtained from a steady-state pressure field for a given combination of fluid injection and offtake (production) boundary conditions (Petrovskyy et al. 2023). Flow paths through connected, highly permeable facies are highlighted by the time-of-flight output.

An orthogonal grid is used for flow diagnostic calculations, to ensure numerical stability (Petrovskyy et al. 2023). Grid cells measure 107.7 m (NE–SW) × 51.3 m (NW–SE) × 1.02 m (thickness), and there are 1 000 000 grid cells in each model. This grid resolution is sufficient to capture the geometry and connectivity of sketched facies units, and, based on sensitivity tests, to calculate flow diagnostics with reasonable accuracy (i.e. <5% difference in flow diagnostic results compared to higher-resolution models with up to 10 000 000 grid cells); a finer grid brings no significant increase in accuracy of the flow diagnostic results in the models considered here. The time of flight and displacement by injected fluid are calculated between opposing faces of each model, with the other four model faces set as no-flow boundaries. Injection and offtake (production) of tracer fluid take place over an entire model face, rather than in wells, which implies that these model boundaries are open to flow. Displacement by injected fluid is calculated in four directions: (1) up depositional dip, from NE to SW (x_max to x_min); (2) down depositional dip, from SW to NE (x_min to x_max); (3) along depositional strike, from NW to SE (y_max to y_min); and (4) along depositional strike, from SE to NW (y_min to y_max). Pressure conditions at the model-face boundaries are determined from the distance between opposing model faces, with pressure differentials of 1 MPa between SW and NE faces (along the x-axis) and 0.48 MPa between SE and NW faces (along the y-axis), respectively, in order to maintain a consistent pressure gradient for all four evaluated flow directions. These pressure gradients are chosen principally to investigate migration and storage of an incompressible tracer fluid over the length scale of the model volume, consistent with the use of flow diagnostics, but are comparable to the pressure gradients expected in regions of a CO₂ storage reservoir that are not near to injection wellbores (cf. De Simone and Krevor 2021; Bump and Hovorka 2024). As outlined below, we use a metric to compare fluid migration and storage in different models that normalizes the tracer-fluid injection rate (relative to the pore volume in the model and to the time required for tracer fluid to flow between the open boundaries at opposite faces of the model), such that this comparative metric is independent of the absolute value of pressure differential. For each flow direction, a simulated observation well is placed in the centre of the model face at which injection occurs; this well shares the same pressure conditions as the rest of the model-face boundary at which it is located.

Future work will focus on incorporating the effects of smaller-scale heterogeneities than those considered here, via the use of flow diagnostics to derive effective properties from models that characterize representative elementary volumes of individual facies (cf. Hossain et al. 2025). Effective properties for multiphase flow derived from models of the characteristic small-scale heterogeneities in each facies will also be the focus of future work, in order to generate results that are specific to CO₂ and H₂ storage.

Three metrics are used to compare the volumetric and flow-diagnostic calculations for different models (cf. Jackson et al. 2022; Alshakri et al. 2023). (1) Total pore volume describes the maximum potential for subsurface fluid storage. We do not use a low-porosity or net-to-gross cut-off(s) to exclude the low permeability units from the calculations of total pore volume. (2) Effective permeability is computed, using flow-based upscaling with no-flow boundaries, over the entire model volume in three orthogonal directions (x, y, z). (3) Pore volume injected (PVI) at breakthrough time provides a proxy for the normalized volume of injected fluid stored in a model due to stratigraphic baffling and trapping. PVI at breakthrough time is computed in each of the four directions described above (i.e. from NE to SW, from SW to NE, from NW to SE and from SE to NW), for tracer fluid between open boundaries at two opposite faces of the model. Snapshots of contacted reservoir volumes at different PVI, which are obtained from the tracer distributions and the time-of-flight field, are also used to visualize flow paths in relation to facies distributions in different models.

The mean value of total pore volume across the suite of 12 models and 8 sets of porosity and permeability values is 1.27 × 10⁹ m³ (corresponding to porosity of 22.5%). There is little variation in total pore volume across the suite (Fig. 5a), with a minimum value of 1.18 × 10⁹ m³ (porosity of 21.0%) for models 4C and 10C (Tables 3 and 4) and a maximum value of 1.35 × 10⁹ m³ (porosity of 23.9%) for models 5G and 11G (Tables 3 and 4). Most of this variation in porosity (8% around mean value; Fig. 5a) is accounted for by the porosity and permeability values assigned to coastal plain and tidal flat deposits (facies Cp1 + T3) (Table 4), although some variation (2% around mean value; Fig. 5a) is attributed to the thickness of tidal and tidally influenced fluvial channel-fill sandbodies (facies T1, T2) (Table 4), plan-view geometry and distribution of tidally influenced fluvial channel-fill sandbodies (facies T2) (Table 4), and porosity and permeability values assigned to distal lower shoreface and offshore deposits (facies W4 + W5) (Table 4). Porosity values for proximal lower shoreface (facies W3), upper shoreface and foreshore (facies W1 + W2) and tidal channel-fill deposits (facies T1) are not varied in our study.

The mean value of total pore volume in the 12 models shows a strong negative correlation with the proportion of coastal plain and tidal flat deposits (facies Cp1 + T3) (R² = 1.00), a strong positive correlation with the proportion of tidally influenced fluvial channel-fill sandbodies (facies T2) (R² = 0.97), and a moderate positive correlation with the proportion of tidal channel-fill sandbodies (facies T1) (R² = 0.62) (Fig. 5b). The proportions of other facies either do not vary in the 12 models (facies W1 + W2, W4 + W5) or have only a minor effect on the mean value of total pore volume and are not highlighted in Figure 5b.

The mean effective permeability values across the suite of 12 models and 8 sets of porosity and permeability values are 270 mD (NE–SW orientation, k_x; Fig. 6a), 531 mD (NW–SE orientation, k_y; Fig. 6d) and 0.04 mD (vertical orientation, k_z; Fig. 6g). Values of k_x vary from 246 mD in model 2C to 295 mD in model 5G (Fig. 6c). Values of k_y tend to be larger than those of k_x, and vary from 400 mD in model 4C to 767 mD in model 5G (Fig. 6f). Values of k_z are much smaller than those of k_x and k_y, and vary from 0.01 mD in models 4C and 5C to 0.08 mD in model 5G (Fig. 6i).

The plan-view geometry and distribution of tidally influenced fluvial channel-fill sandbodies (facies T2) (Table 3; Fig. 4), and the porosity and permeability values assigned to these sandbodies (Table 4) have the greatest influence on values of k_x (6% around mean value; Fig. 6a). Note that the permeability of proximal lower shoreface (facies W3), upper shoreface and foreshore (facies W1 + W2) and tidal channel-fill deposits (facies T1) is not varied. The interactions of the plan-view geometry and distribution of tidally influenced fluvial channel-fill sandbodies (facies T2) with other heterogeneities are also important (Fig. 6a). Mean values of k_x increase progressively from models with scenario 1 of the plan-view distribution and geometry of tidally influenced fluvial channel-fill sandbodies (facies T1) (Fig. 4a, b; models 1, 2, 7, 8 in Table 3) to corresponding models with scenario 2 (Fig. 4c, d; models 3, 4, 9, 10 in Table 3) and scenario 3 (Fig. 4e, f; models 5, 6, 11, 12 in Table 3), reflecting an increase in the up-dip (SW) to down-dip (NE) connectivity of the channelized sandbodies (facies T1). There are only weak correlations between mean values of k_x and the proportions of coastal plain and tidal flat deposits (facies Cp1 + T3) (R² = 0.30), tidally influenced fluvial channel-fill sandbodies (facies T2) (R² = 0.27) and tidal channel-fill sandbodies (facies T1) (R² = 0.27) (Fig. 6b). Thus sandbody connectivity along depositional dip in the suite of models is not controlled predominantly by the proportion of channelized sandbodies, but by their geometry and spatial organization.

The same two heterogeneities, together with the thickness of tidal and tidally influenced fluvial channel-fill sandbodies (facies T1, T2) (Table 4), also control values of k_y (22–23% around mean value; Fig. 6d). Values of k_y increase progressively from models with scenario 2 of the plan-view distribution and geometry of channelized sandbodies (Fig. 4c, d; models 3, 4, 9, 10 in Table 3) to corresponding models with scenario 1 (Fig. 4a, b; models 1, 2, 7, 8 in Table 3) and scenario 3 (Fig. 4e, f; models 5, 6, 11, 12 in Table 3). However, mean values of k_y exhibit moderate-to-strong correlations with the proportions of tidally influenced fluvial channel-fill sandbodies (facies T2) (R² = 0.99), coastal plain and tidal flat deposits (facies Cp1 + T3) (R² = 0.98) and tidal channel-fill sandbodies (facies T1) (R² = 0.50) (Fig. 6e). By implication, sandbody connectivity along depositional strike is controlled predominantly by the proportion of channelized sandbodies.

Values of k_z depend almost entirely on the value of permeability assigned to coastal plain and tidal flat deposits (facies Cp1 + T3) (>150% around mean value; Fig. 6g), and have only weak correlations with facies proportions (R² = 0.03–0.39; Fig. 6h). Coastal plain and tidal flat deposits (facies Cp1 + T3) form a continuous layer in the lower part of each model (Figs 3, 4 and 6i), which determines k_z. Less laterally continuous lenses and wedges of these deposits, and of distal lower shoreface and offshore deposits (facies W4 + W5) that are also assigned reservoir or non-reservoir properties (Table 4), are present in the upper part of each model (Figs 3, 6i), but exert little influence on k_z.

The mean values of PVI at breakthrough vary between the four flow directions within the suite of models. Values are larger for flow along depositional dip, along the longest dimension of the model (Figs 7, 8). For flow up depositional dip (x_max to x_min), mean values of PVI at breakthrough range from 0.52 (model 11) to 0.74 (model 10), with a mean of 0.64 (Fig. 7a). For flow down depositional dip (x_min to x_max), mean values of PVI at breakthrough range from 0.46 (model 9) to 0.67 (model 7), with a mean of 0.54 (Fig. 8a). Values are smaller for flow along depositional strike, along the shortest dimension of the model (Figs 9, 10), but there are small differences depending on whether flow is from NW to SE (y_max to y_min; Fig. 9) or from SE to NW (y_min to y_max; Fig. 10). Mean values of PVI at breakthrough range from 0.24 (model 10) to 0.31 (model 5), with a mean of 0.27, for flow from NW to SE (y_max to y_min; Fig. 9a), and from 0.26 (model 9) to 0.32 (model 5), with a mean of 0.29, for flow from SE to NW (y_min to y_max; Fig. 10a).

For flow up and down depositional dip, most variation in PVI at breakthrough is accounted for by the plan-view geometry and distribution of tidally influenced fluvial channel-fill sandbodies (facies T2) (Table 4), the porosity and permeability values assigned to these sandbodies (facies T2) (Table 4), the thickness of tidal and tidally influenced fluvial channel-fill sandbodies (facies T1, T2) (Table 4), and combinations of these three parameters (Figs 7a, 8a). The connectivity of high-permeability tidal and tidally influenced fluvial channel-fill sandbodies (facies T1, T2) controls flow paths, and thus PVI at breakthrough (Figs 7c–h and 8c–h). Models with low connectivity and/or low permeability of channelized sandbodies exhibit high values of PVI at breakthrough. There are, at best, only weak correlations between mean values of PVI at breakthrough and the proportions of coastal plain and tidal flat deposits (facies Cp1 + T3) (R² = 0.29–0.30), tidally influenced fluvial channel-fill sandbodies (facies T2) (R² = 0.31–0.33), and tidal channel-fill sandbodies (facies T1) (R² = 0.05–0.16) (Figs 7b, 8b). These results are consistent with those for effective permeability along depositional dip, k_x (Fig. 6a–c). Values of PVI at breakthrough are greater for flow up depositional dip (x_max to x_min) (Fig. 7a) than for flow down depositional dip (x_min to x_max) (Fig. 8a), because high-permeability channelized sandbodies (facies T1, T2) are abundant on the up-dip (SW) face of the model but rarely extend to the down-dip (NE) face of the model (perspective views in Fig. 3a–l). As a result, tracer fluid injected along the down-dip (NE) model face generates a piston-like displacement through lower shoreface sandstones (facies W3) in each parasequence, which are all connected to this face, before being focused through the channelized sandbodies (facies T1, T2) (Fig. 7c–h). In contrast, tracer fluid injected along the up-dip (SW) model face moves rapidly through the channelized sandbodies (facies T1, T2) before sweeping heterogeneously through lower shoreface sandstones (facies W3) in only those parasequences that occur at the down-dip terminations of the channelized sandbodies (Fig. 8c–h).

Along depositional strike, as along depositional dip, most variation in PVI at breakthrough is accounted for by the plan-view geometry and distribution of tidally influenced fluvial channel-fill sandbodies (facies T2) (Table 4), the porosity and permeability values assigned to these sandbodies (facies T2) (Table 4), the thickness of tidal and tidally influenced fluvial channel-fill sandbodies (facies T1, T2) (Table 4), and combinations of these three parameters (Figs 9a, 10a). Along depositional strike, there is less variation in the connectivity of high-permeability tidal and tidally influenced fluvial channel-fill sandbodies (facies T1, T2), such that this parameter is uniformly high and exerts less influence on PVI at breakthrough (Figs 9c–f and 10c–h). There is a corresponding increase in the influence of the proportion of these channelized sandbodies (facies T1, T2) for flow from NW to SE (y_max to y_min; Fig. 9). In this direction, there are strong correlations between mean values of PVI at breakthrough and the proportions of coastal plain and tidal flat deposits (facies Cp1 + T3) (R² = 0.67) and tidally influenced fluvial channel-fill sandbodies (facies T2) (R² = 0.74) (Fig. 9b); this correlation reflects the pronounced variations in the proportions of these high-permeability facies along the NW face of the models, along which tracer is injected (right-hand perspective views in Figs 2a–l and 9c–f). In contrast, the proportion of tidal and tidally influenced fluvial channel-fill sandbodies (facies T1, T2) exerts little influence on flow from SE to NW (y_min to y_max; Fig. 10). In this direction, there are no correlations between mean values of PVI at breakthrough and the proportions of coastal plain and tidal flat deposits (facies Cp1 + T3) (R² = 0.04), tidally influenced fluvial channel-fill sandbodies (facies T2) (R² = 0.06), and tidal channel-fill sandbodies (facies T1) (R² = 0.00) (Fig. 10b); this lack of correlation reflects the very limited variations in the proportions of these high-permeability facies along the SE face of the models, along which tracer is injected (left-hand perspective views in Figs 2a–l and 10c–h). These results are consistent with those for effective permeability along depositional strike, k_y (Fig. 6d–f), but demonstrate the variation that can result from the exact location of injector and producer faces in our analysis.

Only some of the heterogeneities and settings considered in our models (Tables 3 and 4) have a significant impact on porosity, effective permeability, and/or PVI at breakthrough. There is only a modest range of porosity (21.0–23.9%) and effective permeability (k_x = 246–295 mD; k_y = 400–767 mD; k_z = 0.01–0.08 mD) in the suite of 12 models (Table 3) and 8 sets of porosity and permeability values (Table 4) investigated for the Cliff House Sandstone example. However, this modest variability results in a comparatively large range of associated PVI at breakthrough values (0.35–0.83 along depositional dip; 0.17–0.42 along depositional strike) and associated flow patterns (Figs 7–10). Heterogeneities that have the largest impact on porosity and effective permeability are the plan-view distribution, thickness, and porosity and permeability characteristics of tidal and tidally influenced channel-fill sandbodies (facies T1, T2), and the porosity and permeability characteristics of coastal plain and tidal flat deposits (facies Cp1 + T3) (Figs 5, 6a, d and g).

PVI at breakthrough and associated flow patterns are controlled by the subdivision of shoreface sandstones into laterally extensive parasequences bounded by offshore mudstones in the down-dip, shallow-marine parts of the unit, and by the spatial distributions, connectivity and permeability of channelized sandbodies in the up-dip, marginal-marine parts of the unit (e.g. Figs 7c–h, 8c–h, 9c–h and 10c–h). There is little variability in the thickness, continuity, internal lithological character and fluid-flow patterns of shoreface parasequences along depositional strike (Figs 3, 9c–h and 10c–h), with the transition to localized mouth bar deposits having negligible influence on PVI at breakthrough (Figs 9a, 10a). Where channelized sandbodies are well connected along depositional strike, there is a positive correlation between PVI at breakthrough and the proportion of these sandbodies (Fig. 9b), although this relationship breaks down where channelized sandbodies are poorly connected (Fig. 10b). Variability in flow patterns is more pronounced along depositional dip, reflecting the localized distribution of channelized sandbodies and their connections with shoreface sandstones at the up-dip pinchouts of shoreface parasequences (Figs 3, 7c–h and 8c–h). Flow down depositional dip results in less uniform displacement in shoreface sandstones, and higher values of PVI at breakthrough, because it is more dependent on the localized connections between up-dip channelized sandbodies and down-dip shoreface sandstones (e.g. compare Fig. 7c–h with Fig. 8c–h). Such heterogeneities have been identified previously as important controls on flow patterns and drainage in generic frameworks for hydrocarbon reservoirs (e.g. ‘jigsaw puzzle’ reservoirs of Weber and Van Geuns 1990), but the detailed stratigraphic architecture and related connectivity between channelized sandbodies and shoreface sheet sandstones are specific to net-transgressive, shallow-marine sandstone reservoirs.

As illustrated in Figures 7–10 and discussed above, our results demonstrate that injected tracer fluid is not dispersed uniformly through the Cliff House Sandstone, even though the unit has a high net-to-gross ratio (80–100% for the suite of models illustrated in Fig. 3 and porosity and permeability values listed in Tables 1 and 4). Injected tracer fluid is dispersed in locally variable and anisotropic patterns, particularly through the channelized sandbodies in the up-dip, marginal-marine parts of the unit. The dispersal patterns of specific injected fluids (e.g. supercritical CO₂, warm or cold water) at specific storage sites will be further complicated by the interaction of structural heterogeneities (e.g. tectonic dip, fault juxtaposition and sealing), regional groundwater flow and well placement with the investigated sedimentological heterogeneities (cf. Skorstad et al. 2008; Jackson et al. 2009), together with the pore-scale interaction of injected and displaced fluid phases with each other and with grain surfaces (Krevor et al. 2015; Hashemi et al. 2021; Qu et al. 2024) and variations in capillary entry pressure associated with centimetre- to decametre-scale heterogeneity (Saadatpoor et al. 2009; Gershenzon et al. 2016; Jackson and Krevor 2020). Differences in the thermal conductivity of sandstones, mudstones and other rock types will also affect thermal energy storage to a modest extent (e.g. Clauser 2021; Baird et al. 2024). Our results thus provide generic insights into the effects of large (metre-to-kilometre) scale sedimentological heterogeneities on subsurface fluid flow and storage, but they need to be extended to include additional heterogeneities, fluid properties and rock-fluid interactions for specific storage sites and projects.

Storage efficiency could be reduced by differential sweep, since less of the available pore space is contacted by the injected fluid (cf. hydrocarbon recovery; Weber and Van Geuns 1990; Tyler and Finley 1991), or increased as a result of stratigraphic trapping at sandbody pinchouts (Flett et al. 2007; Gibson-Poole et al. 2009), and the areal footprint of the injected-fluid plume may differ from pre-injection predictions that do not account for sedimentological heterogeneity. Given the small thickness of individual shoreface parasequences and channelized sandbodies (<15 m; Figs 2, 3) and the sandstone-dominated character of the unit, it is likely that lithological contrasts will not be seismically imaged in corresponding subsurface strata. Differential fluid saturations and pressures that arise during fluid injection may also lie below the resolution of time-lapse seismic data. Wireline-log and core data from wells could enable the identification of shoreface sandstones in laterally extensive parasequences bounded by offshore mudstones and, where penetrated, of channelized sandbodies in mudstone-dominated marginal-marine strata. However, the spatial distribution and connectivity of channelized sandbodies with each other and with shoreface sandstones are poorly constrained by well data, particularly where wells are widely spaced (e.g. Weber and Van Geuns 1990; Tyler and Finley 1991). The resulting uncertainty in sandbody connectivity and reservoir behaviour can be evaluated by reservoir modelling studies (e.g. Jones et al. 1995; Larue and Hovadik 2006; Villamizar et al. 2015).

The Cliff House Sandstone is used as a well-documented outcrop example to investigate the impact of sedimentological heterogeneities on subsurface fluid migration and storage potential in net-transgressive, shallow-marine sandstone reservoirs. Previously published cliff-face correlation panels from the Cliff House Sandstone in Chaco Culture National Historical Park, New Mexico, USA are used to construct a suite of sketch-based reservoir models, which are then compared using flow-diagnostic metrics.

The Cliff House Sandstone outcrop example comprises mudstone-dominated coastal plain, lagoonal and tidal flat deposits that contain channelized tidal and tidally influenced fluvial sandbodies, which are overlain by and pass down depositional dip into 11 aggradationally-to-retrogradationally stacked parasequences composed of wave-dominated shoreface sandstones bounded by offshore mudstones. Six sedimentological heterogeneities are investigated: (1) the plan-view distribution and geometry of channelized tidally influenced fluvial sandbodies; (2) the thickness and vertical connectivity of channelized tidal and tidally influenced fluvial sandbodies; (3) the inclusion of mouth bar sandstones distinct from surrounding proximal lower shoreface sandstones; (4) porosity and permeability of mudstone-dominated coastal plain, lagoonal and tidal flat deposits; (5) porosity and permeability of channelized tidally influenced fluvial sandbodies; and (6) porosity and permeability of distal lower shoreface heteroliths and offshore mudstones.

Total pore volume exhibits little variation between the models. Effective horizontal permeability and PVI at breakthrough, a proxy for flow patterns and stratigraphic trapping potential, are controlled predominantly by: (1) the plan-view distribution and geometry of channelized tidally influenced fluvial sandbodies; (2) the porosity and permeability of channelized tidally influenced fluvial sandbodies; and (3) the interactions of these and other heterogeneities. In detail, flow patterns visualized using time-of-flight flow diagnostics are governed by the stratigraphic compartmentalization of shoreface sandstones by parasequence-bounding offshore mudstones; the spatial distribution, connectivity and permeability of channelized sandbodies; and the local positions of connections between channelized sandbodies and shoreface sandstones. Flow up depositional dip results in more uniform displacement in shoreface sandstones in down-dip locations. Flow down depositional dip occurs through channelized sandbodies in up-dip locations, which are variably connected with down-dip shoreface sandstones that are poorly swept as a result. Flow patterns along depositional strike reflect local variations in the connectivity of channelized sandbodies. Such sandbody connectivity is likely to be poorly constrained in subsurface seismic and well data, and characterization of the resulting uncertainty in subsurface fluid flow and storage potential requires reservoir modelling.

We are grateful for the constructively critical reviews of Gareth Williams and an anonymous reviewer, and the editorial handling of Philip Ringrose. We thank Oliver Jordan, Sanjeev Gupta and Howard Johnson for discussion of the sedimentology and stratigraphy of the Cliff House Sandstone. We also thank Julio Silva, Sicilia Judice, Fazilatur Rahman, Mario Costa Sousa and members of Phase 1 of the Rapid Reservoir Modelling Consortium (Equinor, ExxonMobil Upstream Research Company, Petrobras, Shell, and IBM Research Brazil/IBM Centre for Advanced Studies (CAS) Alberta, Canada) and Phase 2 of the Rapid Reservoir Modelling Consortium (Equinor, ExxonMobil Upstream Research Company, Petrobras, Petronas and Shell) for discussion of sketch-based modelling. Geiger acknowledges partial funding for his Chair from Energi Simulation, and Jacquemyn and Jackson acknowledge support from EPSRC via the ATESHAC project.

XW: conceptualization (supporting), investigation (lead), methodology (supporting), writing – original draft (lead); GJH: conceptualization (lead), data curation (lead), methodology (lead), supervision (lead), writing – review & editing (equal); CJ: conceptualization (supporting), methodology (supporting), supervision (supporting), writing – review & editing (supporting); SH: conceptualization (supporting), methodology (supporting), supervision (supporting), writing – review & editing (supporting); MDJ: conceptualization (supporting), methodology (supporting), supervision (supporting), writing – review & editing (supporting); DP: conceptualization (supporting), methodology (supporting), supervision (supporting), writing – review & editing (supporting); SG: methodology (supporting), supervision (supporting), writing – review & editing (supporting).

This research received no specific grant from any funding agency in the public, commercial, or not-for-profit sectors.

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

The geological data used to construct the models is published in Jordan et al. (2016). The Rapid Reservoir Modelling prototype (executable and source code) used to construct the models is available at: https://bitbucket.org/rapidreservoirmodelling/rrm. The 12 models used in this study are available at: https://figshare.com/articles/dataset/RRM_models_of_Cliff_House_Sandstone_Chaco_Culture_National_Historical_Park_USA_/25669491?file=45833250.