A screening assessment of the impact of sedimentological heterogeneity on CO2 migration and stratigraphic-baffling potential: Johansen and Cook formations, Northern Lights project, offshore Norway

We use a method combining experimental design, sketch-based reservoir modelling, and single-phase flow di- agnostics to rapidly screen the impact of sedimentological heterogeneities that constitute baffles and barriers to CO 2 migration in the Johansen and Cook formations at the Northern Lights CO 2 storage site. The types and spatial organisation of sedimentological heterogeneities in the wave-dominated deltaic sandstones of the Johansen-Cook storage unit are constrained using core data from the 31/5-7 (Eos) well, previous interpretations of seismic data and regional well-log correlations, and outcrop and subsurface analogues. Delta planform geometry, clinoform dip, and facies-association interfingering extent along clinoforms control: (1) the distribution and connectivity of high-permeability medial and proximal delta-front sandstones, (2) effective horizontal and vertical permeability characteristics of the storage unit, and (3) pore volumes injected at breakthrough time (which approximates the efficiency of stratigraphic baffling). In addition, the lateral continuity of carbonate-cemented concretionary layers along transgressive surfaces impacts effective vertical permeability, and bio- turbation intensity impacts effective horizontal and vertical permeability. The combined effects of these and other heterogeneities are also influential. Our results suggest that the baffling effect on CO 2 migration and retention of sedimentological heterogeneity is an important precursor for later capillary, dissolution and mineral trapping.


Introduction
Geological storage of CO 2 has the potential to significantly curb global carbon emissions, in order to meet the targets of the Paris Climate Agreement (Benson et al., 2005). Such storage will require large-scale injection of CO 2 into subsurface saline aquifers and depleted hydrocarbon reservoirs (e.g. Ringrose and Meckel, 2019;Ringrose et al., 2021). In saline aquifer sandstones, CO 2 acts as a low-buoyancy, non-wetting phase that forms a plume after injection (Benson et al., 2005). Initially, short-term (<100 years) CO 2 storage occurs predominantly within structural and stratigraphic traps. Structural and stratigraphic trapping potential, and CO 2 storage volume, are controlled by the porosity and permeability characteristics of the aquifer, the spill points of the trap, and the capillary entry pressure of the seal. High permeability and high effective pore volume allow dissipation of pressure away from the injection well(s) while maintaining a high CO 2 injection rate (Benson et al., 2005;Ajayi et al., 2019). High effective porosity in saline aquifer sandstones also tends to increase the capacity for capillary, dissolution and mineral trapping over medium-term (100-10,000 years) and long-term (>10,000 years) timescales (e.g. Gibson-Poole et al., 2009;Ajayi et al., 2019). However, geological heterogeneities that reduce permeability in the storage unit can be important in two ways. Firstly, they disperse CO 2 and thus decrease the rate at which the CO 2 plume reaches the limits of the storage site, and secondly, they create small-scale stratigraphic baffling and trapping configurations that increase storage efficiency (Flett et al., 2007;Gibson-Poole et al., 2009).
The Northern Lights project involves full-scale CO 2 storage in the saline aquifer sandstones of the Johansen and Cook formations, offshore western Norway (Bergmo et al., 2009;Riis, 2018;Furre et al., 2019) (Fig. 1). These formations were deposited along the eastern margin of the northern North Sea Basin in the early Jurassic (Fig. 1A, B), during a period of thermal subsidence after Permo-Triassic rifting (Husmo et Husmo et al., 2003). Areas with poor data constraints are shown in semi-transparent shading. The study area (Fig. 1C) is located. (B) Sequence stratigraphy and lithostratigraphy of the upper Statfjord Group and lower-to-middle Dunlin Group, offshore Norway, based on the regional maxiumum flooding surfaces of Partington et al. (1993) and megasequences of Steel (1993) (after Husmo et al., 2003). (C) Simplified outline of Exploitation Licence EL001 and Troll Field, highlighting major faults and fault blocks that post-date deposition of the Johansen and Cook formations. The locations of the 31/5-7 (Eos) well and 2D seismic line HRTRE-00212 ( Fig. 2A) are shown.
2003). Subsequent late-Jurassic-to-early-Cretaceous rifting resulted in the generation of north-south-trending, westward-dipping extensional faults that bound eastward-tilted fault blocks (Fig. 1C), thus defining the structural trap of the Troll Field and the structural configuration of the Johansen and Cook formations in the Northern Lights storage site (Fig. 1D). The Johansen and Cook formations occur at 2000-3000 m depth here (Fig. 1D). After seismic mapping of the Johansen and Cook formations over the Northern Lights storage site, confirmation well 31/5-7 (Eos) was drilled and tested between December 2019 and March 2020 to appraise the primary storage unit and its primary seal (www. equinor.com/en/news/20201019-sharing-data-northern-lights.html; Meneguolo et al., 2022). CO 2 injection into the primary storage unit of the Johansen and Cook formations is planned at a rate of 1.5 Mt per year for 25 years, from 2024 (Riis, 2018). The Johansen and Cook formations in the Northern Lights storage site have previously been characterised as shallow-marine sandstones, based on 2D and 3D seismic data, regional well correlations and limited core data from the surrounding area (Sundal et al., 2013(Sundal et al., , 2016Meneguolo et al., 2022) and building on older regional work (Marjanac, 1995;Marjanac & Steel, 1997;Charnock et al., 2001). Given sparse data coverage of the storage site, there remains uncertainty in detailed aspects of the depositional model of the Johansen and Cook formations and the impact of sedimentological heterogeneity on CO 2 migration and storage. Published reservoir modelling and flow simulation studies have generally focussed on CO 2 storage capacity and plume migration in the full storage site, via consideration of simple porosity-depth trends, fault transmissibilities, and CO 2 relative permeability-saturation relationships (Class et al., 2009;Eigestad et al., 2009;Wei and Saaf, 2009). Aspects of sedimentological heterogeneity have been incorporated in various published reservoir modelling studies, including: (1) subdivision of the Johansen Formation aquifer by laterally extensive mudstone layers above marine flooding surfaces (Wei and Saaf, 2009;Sundal et al., 2015), (2) porosity and permeability distributions derived from seismic inversion and attribute analysis (Sundal et al., 2015), (3) permeability anisotropy values (k v /k h ratio) that represent the effects of laterally extensive carbonate-cemented concretions and pervasive micaceous laminae (Sundal et al., 2015), and (4) variogram-based permeability distributions conditioned to interpreted depositional trends, representing variations in sandstone extent and connectivity, and permeability anisotropy values (k v /k h ratio), representing the effects of heterogeneities within sandstones (Meneguolo et al., 2022). However, as outlined in Sections 2 and 3, additional sedimentological heterogeneities are documented in the Johansen and Cook formations and their depositional analogues; the effects of these heterogeneities on CO 2 migration and retention have yet to be fully evaluated.
The aims of this paper are threefold: (1) to evaluate the depositional model of the Johansen and Cook formations CO 2 storage unit using highresolution core photographs from the 31/5-7 (Eos) well; (2) to develop and assess scenarios for sedimentological heterogeneity in the Johansen and Cook formations that reflect uncertainty in the depositional model; and (3) to identify the key sedimentological heterogeneities that control single-phase flow in the Johansen and Cook formations, as a proxy for CO 2 migration and stratigraphic-baffling potential. We address these aims by designing and interpreting numerical experiments that quantify the impact of heterogeneity on flow using three-dimensional (3D) reservoir models. The models are prototypes constructed using a sketchbased approach implemented in Open Source research code (Rapid reservoir modelling, RRM), which allowed scenarios for multiple types of sedimentological heterogeneity to be modelled quickly in a geologically intuitive manner (Costa Jacquemyn et al., 2021a). Scenarios are based on subsurface and outcrop data from potential depositional analogues, in addition to data from the Johansen and Cook formations. The influence of heterogeneity on CO 2 migration and retention is assessed using single-phase flow diagnostics, which are computationally cheap and allow key flow properties and behaviour to be assessed using a single pressure solution (Zhang et al., 2017(Zhang et al., , 2020Jacquemyn et al., 2021b). The resulting methodology is fast, efficient and allows a large number of geological scenarios to be explored prior to more detailed flow simulation. In addition to the CO 2 storage application presented in this paper, the methodology is generally applicable to assessing the impact of heterogeneity on flow in the subsurface (e.g. in groundwater, hydrocarbon and geothermal reservoirs).

Stratigraphic and sedimentological framework from previous work
We use the stratigraphic framework of Sundal et al. (2016) for the Northern Lights storage site and neighbouring areas to define large-scale sedimentological heterogeneity in our study. The framework comprises base-and top-Johansen Formation and top-Cook Formation surfaces that have been mapped over seven 3D seismic cubes and contiguous 2D seismic lines (e.g. Fig. 2A), combined with long-distance correlation between 41 wells with wireline-log data (Sundal et al., 2016). Wireline-log correlations indicate that the Johansen Formation  Fig. 1C) showing the southward (down-dip) pinchout of the Johansen-Cook CO 2 storage unit (shaded yellow and orange) and overlying Sognefjord and Krossfjord formations (shaded blue) (after Sundal et al., 2016). (B, C) Two alternative depositional models for the Johansen Formation: (B) westward progradation of a southward-deflected wave-dominated delta fronted by a spit that lies seaward of a sheltered lagoon (after Sundal et al., 2015Sundal et al., , 2016; and (C) westward progradation of an arcuate wave-dominated delta with subaqueous clinoforms fed by north-to-south sand transport, by analogy with the Sognefjord Formation in the Troll Field (Patruno et al., 2015). The approximate size and position is shown of the models described and analysed in this paper. 4 comprises a single regressive-to-transgressive wedge (Steel, 1993;Marjanac and Steel, 1997;Sundal et al., 2016), bounded at its base and top by the regional J12 and J14 maximum flooding surfaces (Partington et al., 1993;Steel, 1993). Within this wedge, individual wells commonly contain two or three vertically stacked, upward-sandier successions that are interpreted as upward-shallowing parasequences (Sundal et al., 2016). Although vertical trends in parasequence stacking allow the lower, regressive and upper, transgressive parts of the wedge to be correlated regionally, individual parasequences cannot be confidently correlated over distances greater than c. 8 km between some wells (Sundal et al., 2016). Locally, the Johansen Formation regressive-to-transgressive wedge contains a single seismically resolved clinoform set, in which clinoforms dip to the west and to the south ( Fig. 2A) (Sundal et al., 2016). However, clinoforms are not resolved in most of the seismic data, which instead contain reflections that appear to be sub-parallel to the base-and top-Johansen Formation surfaces (Sundal et al., 2016). Only 18 m of core from the Johansen Formation in well 31/2-3 was available to Sundal et al. (2016), which limited the detail of their sedimentological interpretation. Sedimentological interpretation of the limited core data suggested a coarse-grained, high-energy delta front setting (Sundal et al., 2016). More recent interpretation of core and wireline-log data from the well 31/5-7 support the distinction of the lower, regressive and upper, transgressive parts of the Johansen Formation (Meneguolo et al., 2022).
The overlying Cook Formation is bounded at its base and top by the regional J14 and J18 maximum flooding surfaces (Partington et al., 1993;Steel, 1993). Regional well correlation panels have been used to interpret four regressive-to-transgressive tongues in the Cook Formation (Steel, 1993;Marjanac and Steel, 1997), although the formation is thinner than the Johansen Formation over the Troll Field and Northern Lights storage site (Sundal et al., 2016), where it locally contains a single seismically resolved, southward-dipping clinoform set ( Fig. 2A). Regional wireline-log data show three motifs: upward-sandier successions, sharp-based "blocky" sandstones, and sharp-based upward-shalier successions (Marjanac and Steel, 1997). Cores contain a range of facies, including wave-ripple cross-laminated heteroliths and low-angle cross-stratified sandstones assigned to wave-dominated shorefaces (Dreyer and Wiig, 1995;Marjanac and Steel, 1997;Churchill et al., 2017), and cross-bedded sandstones with mud drapes and current-ripple cross-laminated heteroliths assigned to tidal channels, estuarine deposits and tidal bars (Dreyer and Wiig, 1995;Marjanac and Steel, 1997;Gupta and Johnson, 2001;Churchill et al., 2017), including in the 31/5-7 well (Meneguolo et al., 2022).

Analysis of data from 31/5-7 (Eos) well
The 31/5-7 (Eos) well was drilled after Sundal et al. (2016) completed their study, and provides continuous core and wireline logs through the Johansen and Cook formations. Meneguolo et al. (2022) provide the first sedimentological description of cores and wireline-log data from the 31/5-7 well. Data from the well are publically available (www.equinor.com/en/news/20201019-sharing-data-northern-lights. html). We used high-resolution photographs of 138 m of core from the well to interpret depositional facies associations, facies successions, stratigraphic surfaces and associated small-and medium-scale sedimentological heterogeneities. Our facies analysis of core photographs is based on observations of lithology, grain size, sedimentary structures, body fossils, and trace-fossil assemblages, and is independent of the analysis conducted by Meneguolo et al. (2022). The Bioturbation Index (BI) scheme of Taylor and Goldring (1993) is used to describe bioturbation intensity, which ranges from 0 (no bioturbation) to 6 (completely bioturbated). Porosity and permeability of the depositional facies were characterised using publically available conventional core analysis data (www.equinor.com/en/news/20201019-sharing-data-nor thern-lights.html) from core plugs taken at 0.25 m spacing in the cored interval.

Subsurface and outcrop analogue data
The sedimentological characteristics and depositional models of the Johansen and Cook formations in the Northern Lights storage site are similar in many aspects to those of the middle-to-upper Jurassic Krossfjord, Fensfjord and Sognefjord formations in the nearby Troll Field, which contains high-resolution 3D seismic data and core, wireline-log, biostratigraphic and production data from numerous wells (Fig. 1C, D). Based on published descriptions of these formations (e.g. Stewart et al., 1995;Holgate et al., 2013Holgate et al., , 2015Patruno et al., 2015;Sundal et al., 2016), similarities include: (1) the areal extent and thickness of seismically imaged clinoform sets; (2) the occurrence of coarsening-upward, sandstone-rich successions in wireline logs, the thickest of which correspond to seismically imaged clinoform sets; (3) the grain size range, proportion of sandstone and sedimentary structures observed in core; (4) the abundance of carbonate-cemented concretions; and (5) geological context, characterised by shallow-marine deposition on a narrow shelf adjacent to the eastern margin of the North Sea rift basin. Depositional models of the Johansen, Cook, Krossfjord, Fensfjord and Sognefjord formations have developed through time, as the quantity and quality of available data have increased and as concepts for interpretation have evolved. For example, early depositional models of the Johansen Formation portray a simple, linear-to-arcuate shoreline (Fig. 1A, after Husmo et al., 2003) while later depositional models show a morphologically complex shoreline with a wave-dominated delta and contiguous spit (Fig. 2B, after Sundal et al., 2015. Additional constraints on sedimentological heterogeneity at length scales larger than core and wireline-log sampling but smaller than seismic resolution are taken from outcrop analogues. In particular, the lateral extent and continuity of carbonate-cemented concretions is uncertain in the Johansen and Cook formations; these may form concretionary layers of high lateral continuity (e.g. in the Bridport Sand Formation; Hampson et al., 2015;Sundal et al., 2015) or much less extensive, discoidal concretions (e.g. in the Frewens Sandstone Formation; Dutton et al., 2000).

Sedimentological heterogeneity in the Johansen and Cook formations
As the first part of our analysis, we characterise sedimentological heterogeneity in the Johansen and Cook formations in the Northern Lights storage site, from larger length scales imaged in seismic data to smaller length scales sampled in core. Our characterisation includes original interpretations of facies associations and facies successions in the 31/5-7 (Eos) well, and also draws on subsurface and outcrop analogues.

Seismically imaged stratigraphic architecture
At the largest length scales considered for the Johansen and Cook formations (1s to 10s km laterally, 10s to 100s m vertically), sedimentological heterogeneity comprises seismically imaged stratigraphic architecture, geomorphology and porosity trends (Sundal et al., 2015(Sundal et al., , 2016. Seismic reflection geometries and configurations indicate westward and southward progradation of clinoforms in both formations ( Fig. 2A) (Sundal et al., 2016). Clinoform dips are calculated using a vertical distance of 150 ms two-way time (TWT) (e.g. Fig. 2A) to equate to 180-200 m of combined thickness for the Amunsden, Johansen, Burton and Cook formations (e.g. Fig. 3A). Steep clinoform geometries (up to 16 ) are observed locally (Sundal et al., 2016), although gentler dips (1-3 ) are more common (e.g. clinoform sets highlighted in Fig. 2A). Where clinoforms are not seismically resolved, reflections appear to be shingled or sub-parallel to the base-and top-Johansen Formation surfaces (Sundal et al., 2016). The sedimentologically similar Sognefjord Formation occurs at shallower depths in the Troll Field, and is thus imaged at higher resolution; here, similar sub-horizontal reflector packages correspond to coarsening-upward successions that thicken and pass laterally into shingled clinoform sets (Patruno et al., 2015). By analogy, shingled or sub-horizontal reflector packages in the Johansen and Cook formations may represent one or more vertically stacked clinoform sets that are too thin to resolve fully and have uncertain orientations. Seismic amplitude and acoustic impedance volumes calibrated to wireline-log data indicate that the lower Johansen Formation is associated with high porosity beneath the Troll Field, and that high porosities extend south of the eastern part of the Troll Field (Sundal et al., 2016). Porosity decreases markedly to the west across a gently arcuate, northnorthwest-southsoutheast-oriented boundary that coincides approximately with the regional palaeoseaward pinchout of Johansen Formation sandstones (Fig. 1A). In the upper Johansen Formation, high porosity in the Troll Field is contiguous with an elongate body of high porosity that is oriented northnorthwest-southsoutheast and protrudes south of the western part of the Troll Field (Sundal et al., 2016). These porosity distributions are interpreted to record westward progradation of an arcuate wave-dominated delta in the lower Johansen Formation, followed by progradation and aggradation in the upper Johansen Formation of a southward-deflected wave-dominated delta fronted by a spit that lay seaward of a sheltered lagoon (Fig. 2B) (Sundal et al., 2015(Sundal et al., , 2016. Seismically derived porosity distributions in the Cook Formation are not published. Similar seismic-stratigraphic pinchout geometries, clinoform orientations and porosity distributions have been documented in the overlying Krossfjord, Fensfjord and Sognefjord formations in the Troll Field (Dreyer et al., 2005;Patruno et al., 2015), although these units are thicker because they contain a greater number of vertically stacked clinoform sets than the Johansen and Cook formations (e.g. Fig. 2A). Depositional models developed for the data-rich Krossfjord, Fensfjord and Sognefjord formations are therefore a useful comparison of those developed for the data-poor Johansen and Cook formations. The Fensfjord, Krossfjord and Sognefjord formations are interpreted as the deposits of westward-prograding wave-dominated deltas with contiguous southward-prograding spits (Dreyer et al., 2005;Holgate et al., 2015;cf. Fig. 2B) or westward-prograding wave-dominated deltas with subaqueous clinoforms fed by north-to-south transport of sand (Patruno et al., 2015;cf. Fig. 2C). Core data indicate that relatively mud-rich, tidally influenced deposits lie eastward (palaeolandward) of the interpreted spit or southward-deflected subaqueous clinoforms, because wave energy was lower during early delta progradation (Stewart et al., 1995;Holgate et al., 2013Holgate et al., , 2015Patruno et al., 2015). The Krossfjord, Fensfjord and Sognefjord deltas each underwent multiple episodes of regression and transgression, and contain multiple vertically-stacked parasequences (Stewart et al., 1995;Holgate et al., 2013Holgate et al., , 2015Patruno et al., 2015). Clinoforms occur in each parasequence, but are most clearly seismically imaged where a parasequence extends and thickens beyond the downdip pinchout of underlying parasequences (Holgate et al., 2015;Patruno et al., 2015). Clinoform dips vary in each parasequence, indicating that clinoforms steepened as the delta evolved from relatively mud-rich and tidally influenced during early progradation to mud-poor and wave-dominated during late progradation (Patruno et al., 2015).
In summary, depositional models of the Johansen and Cook formations and analogous Krossfjord, Fensfjord and Sognefjord formations contain two contrasting components (Fig. 2B, C). The first component is the westward (palaeoseaward) progradation of a linear-to-arcuate deltaic shoreline that became progressively more wave-dominated and less tidally influenced through time. The second component is the southward (shoreline-parallel) transport of sand in the wave-dominated delta.

Facies associations and facies successions in cores and wireline logs
Wireline-log data from individual wells indicate that the Johansen and Cook formations contain several vertically stacked, upward-sandier successions (Marjanac and Steel, 1997;Sundal et al., 2016). Using core data from the 31/5-7 (Eos) well, we identify and interpret five sedimentological facies associations within the formations, document their vertical stacking into facies successions, and outline their wireline-log characteristics (Figs. 3,4). Porosity and permeability characteristics of the facies associations are summarised in Table 1. Note that FA2 was omitted from our models because it is present in only a small proportion of the core (Fig. 3), and its porosity and permeability are measured in a small, and potentially unrepresentative, number of samples (Table 1).
The first facies association (FA1) comprises siltstones that contain thin (1-5 cm) beds of parallel-laminated and asymmetrical-and symmetrical-ripple cross-laminated fine-grained sandstones (Figs. 3, 4A, B). Bioturbation varies from low (BI of 2) to intense (BI of 5) by a diverse trace fossil assemblage (Planolites, Palaeophycus, Teichichnus, Rhizocorallium, Thalassinoides, mud-filled Arenicolites, Phycosiphon, Chondrites). Thick (>3 m) intervals of FA1 occur in the lower part of upward-sandier wireline-log successions, although thin (<3 m) intervals occur in their middle and upper parts. Intervals of FA1 are intercalated with and gradationally overlain by the second and third facies associations (FA2, FA3) (Fig. 3). FA1 is characterised by moderate-to-high gamma ray values (70-140 API), and low density (RHOB) and high neutron porosity (NPHI) curves that are widely separated. FA1 is interpreted as offshore transition zone deposits, characterised by deposition of mud from suspension and episodic influxes and reworking of sand by currents and waves (cf. "A" and "D1" facies of Dreyer and Wiig, 1995, "distal transition zone" deposits of Churchill et al., 2017). Trace fossil assemblages constitute the archetypal Cruziana ichnofacies, developed under fully marine conditions (MacEachern and Bann, 2008), with variations in bioturbation intensity recording bed-scale variations in sedimentation rate (Gingras et al., 2011). Most intervals of this facies association are interpreted as tidal flat deposits by Meneguolo et al. (2022).
The second facies association (FA2) comprises intercalated siltstones, thin (1-5 cm) beds of current-and wave-ripple cross-laminated fine-grained sandstones, and thick (5-20 cm) structureless and graded beds of poorly sorted fine-grained silty sandstones (Figs. 3, 4C). Bioturbation intensity is low to intense (BI of 2-5) in the thinly bedded intervals, with the same diverse trace fossil assemblage as FA1, but absent (BI of 0) in the thick sandstone beds. Intervals of FA2 are intercalated with FA1 and FA3 (Fig. 3). Intervals of FA2 are thin (<0.3 m), and thus not well resolved in wireline-log data, but they have moderate gamma ray values (70-90 API). FA2 is interpreted as offshore transition zone deposits, as FA1, but containing gravity-flow event beds. In the context of a delta (Fig. 2B, C), gravity-flow event beds most likely record sediment-laden hyperpycnal flows developed during periods of high fluvial discharge (Mulder and Syvitski, 1995), potentially enhanced by storm waves (Macquaker et al., 2010). FA2 therefore constitutes distal delta front and prodelta deposits. Most intervals of this facies association are interpreted to occur in tidal flat deposits by Meneguolo et al. (2022).
Facies association 3 (FA3) consists of highly to intensely bioturbated (BI of 4-5), micaceous sandstones. The facies association contains a diverse trace fossil assemblage (Planolites, Palaeophycus, Teichichnus, Thalassinoides, Rosselia, Asterosoma, mud-filled Arenicolites, Phycosiphon) (Figs. 3, 4D-F). Where preserved, primary sedimentary structures comprise asymmetrical-ripple cross-lamination with mudstone drapes along foresets and cross-set boundaries (Figs. 3, 4F). FA3 gradationally overlies FA1 and FA2, and is gradationally overlain by FA4 (Fig. 3). FA3 is characterised by moderate gamma ray values (60-80 API), and RHOB and NPHI curves that are closely spaced and/or overlap. FA3 is interpreted as medial delta front deposits, characterised by episodic influxes of sand by currents with intervening periods of extensive biogenic reworking (cf. "B" facies of Dreyer and Wiig 1995). The occurrence of mudstone-draped current-generated ripples implies that repeated fluctuations in flow velocity, potentially due to tides, played a role in episodic sand transport. Direct evidence of wave action, such as symmetrical ripples and hummocky cross-stratification, is absent, but the abundance of sand is consistent with offshore transport of suspended mud above wave base. Trace fossil assemblages constitute the proximal Facies association 4 (FA4) comprises cross-bedded, medium-to coarse-grained sandstone (Figs. 3, 4G, H). Cross-sets are c. 10 cm thick and contain sparse mudstone drapes along toesets and occasionally foresets (Fig. 4H). Bioturbation intensity is sparse to low (BI of 1-2) by a trace fossil assemblage of low diversity (Skolithos, Planolites, mottling). Shell fragments and other bioclasts occur locally along erosional bed bases (Fig. 4G). FA4 overlies FA3, either gradationally or across a sharp bed boundary, and is intercalated with FA5 ( Fig. 3). FA4 has uniformly low-to-moderate gamma ray values (40-60 API) and cross over of RHOB and NPHI curves, consistent with relatively clean sandstone. Calcite cemented concretions are relatively abundant in FA4, marked by high RHOB and low NPHI values. FA4 is interpreted as proximal delta front deposits, potentially including distributary channels, characterised by migration of sand dunes in response to unidirectional currents (Sundal et al., 2016; cf. "E" and "F" facies of Dreyer and Wiig, 1995, "upper shoreface" deposits of Churchill et al., 2017). The coarse, moderately sorted texture of sandstones implies fluvial sediment supply. Repeated variations in flow velocity, potentially indicating tidal influence, are recorded by mudstone drapes along cross-set foresets and toesets, while the scarcity of mudstone is consistent with offshore transport of suspended mud above wave base. The trace fossil assemblage constitutes the Skolithos ichnofacies (MacEachern and Bann, 2008), while calcite concretions are interpreted to result from local redistribution of bioclastic shell material (Sundal et al., 2016). Intervals of facies association 4 are interpreted as subaqueous distributary channel, mouth bar and compound dune deposits by Meneguolo et al. (2022).
The final facies association (FA5) consists of cross-bedded, mediumgrained to granular sandstone (Figs. 3, 4I, J). Cross sets are several tens of centimetres up to one metre in thickness, exhibit clear upward steepening of laminae from toesets to foresets. Toesets are commonly rhythmically draped by mudstone, while mudstone drapes occur more rarely along foresets (Fig. 4I). Cross-set bases are lined by coarse, bioclastic lags, and are variably calcite cemented (Fig. 4J). Bioturbation is absent to sparse (BI of 0-1) by a restricted trace fossil assemblage (Planolites, Palaeophycus, mud-filled Arenicolites). FA5 is intercalated with FA4 ( Fig. 3). FA5 exhibits a low gamma ray values (40-50 API) and cross over of RHOB and NPHI curves, consistent with clean sandstone. FA5 is interpreted as proximal delta front deposits, as FA4, but with more pronounced tidal influence recorded by more abundant, rhythmic mudstone drapes along cross-set foresets and toesets. Intervals of facies association 5 are interpreted as subaqueous distributary channel, mouth bar, tidal inlet, flood tidal delta, tidal channel, upper subtidal and compound dune deposits by Meneguolo et al. (2022).
The five facies associations, FA1-5, are consistent with both seismicscale depositional models of the Johansen Formation ( Fig. 2B, C). Wireline-log data from the 31/5-7 (Eos) well indicate that the Amunsden, Johansen, Burton and Cook formations together comprise three large-scale successions (51-77 m thick) that each contain an overall coarsening-upward trend from mudstone into overlying sandstone, capped by a thinner fining-upward trend (Fig. 3A). Each of these successions is described below. (1) The first succession comprises the Amunsden Formation and lower Johansen Formation (2832-2753 m in The third succession is interpreted to record progradation and subsequent aggradation of the "Cook delta" (cf. stacked tidal sand ridges interpreted by Meneguolo et al. (2022).

Outcrop analogues of core-scale heterogeneities
Many of the heterogeneities sampled by the 31/5-7 (Eos) well are not represented adequately in cores and core plugs, but require larger rock volumes for characterisation. These representative rock volumes can be provided by outcrop analogue data, although the selection of one or more appropriate analogues may not be straightforward. Interfingering between facies associations is observed to occur over distances of tens to several hundreds of metres in various shallow-marine and deltaic sandstone reservoir analogues (e.g. Kjønsvik et al., 1994, Howell et al., 2008a. Facies-association interfingering is observed to occur along clinoforms in numerous shallow-marine and deltaic sandstones at outcrop (e.g. Willis et al., 1999;Howell et al., 2008b;Sech et al., 2009), and is inferred to conform to clinoforms in the subsurface Krossfjord, Fensfjord and Sognefjord formations (Dreyer et al., 2005;Holgate et al., 2013Holgate et al., , 2015Patruno et al., 2015). Different outcrop analogues demonstrate that carbonate cement can occur as isolated concretions (e.g. Frewens Sandstone; Dutton et al., 2000) or as more continuous concretionary layers (e.g. Bridport Sand Formation; Hampson et al., 2015) in shallow-marine and deltaic sandstones. In general, carbonate-cemented concretions and concretionary layers developed directly below transgressive surfaces and flooding surfaces are more laterally extensive than those occurring within the intervening facies successions (e.g. Taylor et al., 1995;Sixsmith et al., 2008;Morad et al., 2010;Martinius, 2017). The abundance, lateral extent and geometrical configuration of mudstone beds and drapes are highly variable in heterolithic tidal sandstones, but control the effective permeability and permeability anisotropy in these sandstones. Their effects have been characterised in several outcrop analogues (Jackson et al., 2005;Massart et al., 2016) and related forward stratigraphic modelling studies (Ringrose et al., 2005). The effects of bioturbation fabric on effective permeability properties have been similarly studied using outcrop analogues and samples (Gingras et al., 2012).

Hierarchical arrangement of heterogeneity
The multi-scale heterogeneity described above is synthesised into an interpreted hierarchy of sedimentological heterogeneity in the Johansen and Cook formations (Fig. 5), which provides a framework to identify, organise and model these heterogeneities. The hierarchy of heterogeneity is modified from similar schemes for other shallow-marine strata (Kjønsvik et al., 1994;Sech et al., 2009;Graham et al., 2015). Heterogeneity resolved in seismic data ( Fig. 5A-D), sampled in wireline logs ( Fig. 5B-E) and captured in core photographs (Fig. 5E) is considered in the hierarchy. Heterogeneity at the smaller scale of pores and grains is important in assessing textural and mineralogical controls on CO 2 trapping in the Johansen and Cook formations (Sundal and Hellevang, 2019), but lies beyond the scope of this study. Note that our hierarchical scheme does not accommodate the channelised and landward-tapering, wedge-shaped elements (e.g. distributary channels, tidal channels, tidal inlets, flood tidal deltas) interpreted by Meneguolo et al. (2022).

Modelling methodology
We use a screening approach to assess the influence of sedimentological heterogeneity on CO 2 migration and storage in the Johansen-Cook CO 2 storage unit. Our methodology has three key elements, which in combination require only a small number of models (via experimental design) that can be constructed quickly in a geologically intuitive way (via sketch-based modelling) and analysed in a computationally cheap manner (via single-phase flow diagnostics). These three elements and other aspects of the modelling methods are described below. Our approach is scenario-based and deterministic (Bentley and Smith, 2008), and it is appropriate for screening the most influential sedimentological heterogeneities prior to more detailed follow-up studies, including stochastic modelling of selected scenarios.

Design of modelling experiment
Eight sedimentological heterogeneities were selected for investigation in this study (Table 2). Experimental design techniques (Box et al., 1978;Damsleth et al., 1992;Kjønsvik et al., 1994;White and Royer, 2003) and analysis of variance (Box et al., 1978) were used to efficiently explore the resulting parameter space. A two-level fractional-factorial design (2 IV 8-3 ) was used, in which each factor (heterogeneity) can take either a low or high setting, with settings chosen to reflect the range of uncertainty in the Johansen-Cook CO 2 storage unit (Table 2), as outlined below (Section 4.2). The experimental design allows us to efficiently quantify the effect of varying each factor from its low setting to its high setting. Its resolution IV design ensures that the main effects, due to each of the eight studied heterogeneities, are not confounded with interactions between two heterogeneities (Box et al., 1978). The experimental design requires 32 models to be constructed for the screening study (Table 3). This experimental design has limitations, however; additional models are required to assess systematically the effects of interactions between two or more heterogeneities (Box et al., 1978), and the heterogeneities are assumed to vary independently of each other. Note also that only heterogeneities captured in our hierarchical scheme (Fig. 5) are considered in the modelling experiments.

Modelled heterogeneities and settings
The sedimentological heterogeneities under investigation are listed below in order of decreasing length scale (Fig. 5), and summarised in Table 2. Values of low and high settings for the sedimentological heterogeneities are chosen to encompass uncertainty in their interpretation given sparse data constraints in the Johansen-Cook CO 2 storage unit. Combinations of parameter settings in each model are shown in Table 3.
At the largest length scale, the "lower Johansen, upper Johansen and Cook deltas" have been interpreted to exhibit either a simple arcuate geometry (Figs. 2C, 5A) or a more complex deflected, elongate geometry (Figs. 2B, 5B). These contrasting configurations are taken as the two settings for delta planform geometry (Table 2). Clinoforms resolved in seismic data (e.g. Fig. 2A) vary in dip; 1and 3were used as representative low and high settings, respectively (Table 2; cf. Section 3.1). Steeper clinoform dips are noted locally (up to 16 ; Section 3.1), thus our choice of 3as a high setting may underestimate the effects of clinoform dip.
Carbonate-cemented concretions are common in the Johansen and Cook formations, along transgressive surfaces and maximum flooding surfaces and also in between such surfaces (Fig. 3). Carbonate-cemented concretions along transgressive surfaces and maximum flooding surfaces are interpreted to be laterally extensive sheets that contain holes ( Fig. 5D; cf. Gibbons et al., 1993). The lateral continuity of concretionary sheets between holes is taken as c. 100 m and c. 1000 m for the low and high settings, respectively (Table 2), and the concretionary sheets are assumed to be non-porous and impermeable (i.e. they are assigned porosity and permeability values of 0% and 0 mD in our models).
Carbonate-cemented concretions in between transgressive surfaces are interpreted to have limited lateral extents and variable vertical spacings (Fig. 5E). Mean values of 5.1 m, measured in the Frewens Sandstone outcrop analogue (Dutton et al., 2000), and 32.7 m, measured in the Bridport Sand Formation outcrop analogue , are used as the low and high settings for the concretion lateral extent (Table 2). Based on variability in vertical spacing noted in the core (Fig. 3) (1985), and assuming that the vertical and horizontal permeability of the non-cemented sandstone matrix are the same: where F s is the fraction of cemented sandstone, s is the inverse of the mean vertical spacing of concretions, and l is the mean length of concretions. Assuming a mean concretion thickness of 30 cm, values of F s are 0.3 and 0.05 for mean vertical concretion spacings of 1 m and 6 m, respectively. The resulting estimates of k v /k h ratio that account for carbonate-cemented concretions in between transgressive surfaces are: (1) 0.1 for the low settings of mean lateral extent and vertical spacing of concretions; (2) 0.6 for the low setting of mean lateral extent and high setting of vertical spacing; (3) 0.005 for the high setting of mean lateral extent and low setting of vertical spacing; and (4) 0.1 for the high settings of mean lateral extent and vertical spacing of concretions. We also consider the effects of variable mudstone drape continuity and extent in heterolithic cross-bedded sandstones (FA4, FA5) by estimating values of k v /k h ratio that are assigned to grid cells (Table 2; Section 4.3). The mean proportions of sandstone in FA4 and FA5 are estimated using core photos at 1.00 and 0.97, respectively representing low and high mudstone drape continuity and extent, and mudstone drapes generally occur in cross-bed toesets (Fig. 4G-J). Corresponding outcrop-based reservoir models of similar heterolithic, trough crossbedded sandstones are used to estimate values of k v /k h ratio of 1.0 and 0.1 for low and high settings of mudstone drape continuity and extent, respectively (Fig. 6, after Massart et al., 2016). The low setting corresponds to poorly amalgamated mud patches >10 cm in length draped along the cross-bed toesets (Massart et al., 2016), representing FA5, whereas the high setting corresponds to no mud drapes (Massart et al., 2016), representing FA4. It should be noted that a few cross-beds Table 2 Summary of investigated sedimentological heterogeneities (factors) and their low and high settings in the screening study. The impact of these eight factors on simulated fluid flow is assessed by observing the percentage change in average response when each factor is varied from its low setting to its high setting.

Table 3
Design and parameter settings (l = low, h = high; Table 2) of the 32 models constructed for the screening study. in FA4 contain lower proportions of sandstone and exhibit mudstone drapes along both foresets and toesets (e.g. in parts of Fig. 4I); k v /k h ratio in these cross-beds may be an order of magnitude lower than the estimated value used here for the lower setting (Fig. 6), but these cross-beds are not considered representative of the facies associations. At the smallest length scale, bioturbation acts to destroy physical sedimentary structures (e.g. cross-bedding, lamination) and may either homogenise the sediment or introduce a new permeability fabric defined by the connectivity of distinct burrow networks (Gingras et al., 2012). In FA1, FA2 and FA3 (and to a lesser extent FA4 and FA5), bioturbation in the Johansen and Cook formations has served to homogenise the sediment. In such cases, the geometric mean of core-plug permeability measurements provides an appropriate value of effective permeability for the facies association (e.g. La Croix et al., 2017) at grid-cell scale (Section 4.3). We use the arithmetic and geometric means of core-plug permeability measurements (Table 1), respectively, as the low and high settings for bioturbation intensity in each facies association (Table 2). Thus, for the low setting, the arithmetic mean is used to represent the effects of less intense bioturbation on horizontal permeability (k h ), and the k v /k h ratio representing variable mudstone drape continuity and extent is retained in FA4 and FA5. For the high setting, the geometric mean is used to represent the effects of more intense bioturbation on both k h and k v .

Sketch-based construction of reservoir models
Reservoir models are constructed using an intuitive, sketch-based approach that allows geological concepts and scenarios to be rapidly captured by non-experts (Costa Jacquemyn et al., 2021a). Our approach adapts Sketch-Based Interface and Modelling (SBIM) methods developed for non-geological CAD and CFD applications (e.g. Olsen et al., 2009). Geological architectures and heterogeneities (e.g. faults, sequence stratigraphic surfaces, facies boundaries, diagenetic boundaries) are represented by surfaces (cf. Denver and Phillips, 1990;Hamilton and Jones, 1992), which 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 widely used by geoscientists to conceptualise and represent geological interpretations (e.g. in maps, cross-sections and block diagrams; Fig. 5), and can be generated and manipulated readily using SBIM methods. Interactions between SBIM-generated surfaces are controlled by operators that necessitate geological viability in the resulting models. Our sketch-based modelling approach is implemented in Open Source research code (Rapid reservoir modelling, RRM), and is described in detail in Costa  and Jacquemyn et al. (2021a). 32 models were constructed, following the experimental design (Table 3). Our approach allows surfaces to be sketched in any order (Jacquemyn et al., 2021a), although models were sketched in order of the hierarchy of heterogeneity (Fig. 5) in this case; heterogeneities at large length scales were sketched first, followed by heterogeneities at progressively smaller length scales. Four heterogeneities, which occur at relatively large length scales, were sketched explicitly in the models: (1) planform geometry; (2) clinoform dip; (3) lateral continuity of carbonate-cemented concretions along transgressive surfaces and maximum flooding surfaces; and (4) dip extent of facies-association interfingering down clinoforms (Table 2). Where combinations of the same settings for these heterogeneities recur in two or models, the surfaces that represent the heterogeneities were re-used. Re-use of the surfaces aids rapid model construction, and also retains consistency in geometrical architectures and geological-domain volumes between models. Four further heterogeneities, which occur at relatively small length scales, were represented implicitly by assigning different values of effective permeability (k h , k v /k h ) to geological domains: (5) mean vertical spacing of carbonate-cemented concretions in between transgressive surfaces; (6) mean lateral extent of carbonate-cemented concretions in between transgressive surfaces; (7) mudstone drape continuity and extent in heterolithic cross-bedded sandstones (FA4, FA5); and (8) bioturbation intensity (Table 2). Because these four heterogeneities are not explicitly represented by sketched surfaces, their setting does not affect geometrical architectures and geological-domain volumes.
Each model has dimensions of 2500 m (west-east) x 3000 m (northsouth) x 200 m (thickness) (Fig. 7). The areal extent of the models is thus smaller than that of the Johansen-Cook storage unit in the Northern Lights storage site (Fig. 1C, D), and the models are intended to investigate only a representative part of the storage site. Structural elements of the Northern Lights storage site, including faults and tectonic dip, are not incorporated in the models so that the effects of sedimentological heterogeneity are not obscured; additional models would be required to  Fig. 7). Data are derived from small models (3 m x 3 m x 1 m) constructed using outcrop data, and the proportion of foreset-toeset surfaces that are draped by mudstone patches ranges from 100% (at a sandstone proportion of 0.74) to 0% (for a sandstone proportion of 1.00).
investigate the interaction between sedimentological and structural heterogeneity. Models are conditioned to sedimentological data (i.e. maximum flooding surfaces, transgressive surfaces, facies successions, and distributions of carbonate-cemented concretions; Fig. 3) and petrophysical data (i.e. porosity, permeability; Table 1) from the 31/5-7 (Eos) well, which is included in the model volume. Fig. 7 shows an example of sketched surfaces and surface-bounded domains taken from model number 32 (Table 3). Models are generated without reference to an underlying grid (cf. Jacquemyn et al., 2019), although a grid is created to visualise them (Fig. 7A) or to perform numerical calculations (Zhang et al., 2018). Grid-cell size and grid resolution are discussed in Section 4.4.

Flow diagnostics
Sketched models are quantitatively analysed to determine key reservoir properties of the Johansen-Cook CO 2 storage unit. Volumetric calculations can be obtained directly from the sketched model after porosity values have been assigned to different facies associations (Table 1). Our implementation of SBIM for reservoir modelling is integrated with computationally cheap flow diagnostics, which allow key flow properties and behaviours to be assessed rapidly using a reducedphysics, single-phase pressure solution (Shahvali et al., 2012;Møyner et al., 2014). The key outputs of flow diagnostics are 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) wells. These outputs highlight flow paths through connected, highly permeable facies associations. Additional parameters such as sweep efficiency, proxies for storage efficiency and the Lorenz coefficient can then be calculated (Møyner et al., 2014). The effects of fluid compressibility, transient flow (e.g. gravity segregation), multiphase fluid interaction (e.g. relative permeability, dissolution) and fluid-rock interactions (e.g. mineralisation) are not included in flow diagnostics. Our use of flow diagnostics therefore allows rapid, preliminary assessment of the impact of different geological concepts and scenarios on flow properties and behaviours, as a precursor to more detailed but time-consuming full-physics, multiphase simulations (Zhang et al., 2017(Zhang et al., , 2020Jacquemyn et al., 2021b). In the present study, we use flow diagnostics to assess the effects of sedimentological heterogeneity and related stratigraphic baffling and trapping on CO 2 migration and storage over the length scales of heterogeneities represented in our models (Figs. 6, 7). Capillary, dissolution and mineral trapping mechanisms are not simulated by flow diagnostics, and are therefore not accounted for in our results.
Volumetric and flow-diagnostic calculations require a grid to be generated for our models, while flow-diagnostic calculations require specification of boundary conditions and the number, location, perforation interval and bottom-hole pressure of injection and offtake wells. The six faces of each model are set as no-flow boundaries, and a single injection well and single offtake well are used (Fig. 8). Wells are placed in the centre of either the west and east faces (coloured red in Fig. 8A, B) or the north and south faces (coloured blue in Fig. 8A, B) of the models, and they are perforated over the entire model thickness. The pressure differential between injection and offtake wells is set at 100 bar. These well placements, perforations and bottom-hole pressure constraints are not indicative of the plan for development and operation for the Northern Lights storage site, but are instead chosen to investigate CO 2 plume migration, dissipation and potential trapping over the length scale of the entire model volume.
An orthogonal grid is used, to ensure numerical stability.   (Table 3). Cross-sections intersect the 31/5-7 (Eos) well in the model volume. Maximum flooding surfaces, transgressive surfaces, facies-association boundaries, and the boundaries of carbonate-cemented concretions along transgressive surfaces and maximum flooding surfaces are represented by sketched surfaces (A, B). Facies associations and carbonate-cemented concretions along transgressive surfaces and maximum flooding surfaces are represented by geological domains bounded by sketched surfaces (C, D). Model grids are shown in 3D perspective views for visualisation purposes only. See Fig. 3 for key to facies association colours. that this grid resolution is sufficient to calculate flow diagnostics with reasonable accurately (e.g. for south-to-north breakthrough time for tracer between injection and offtake wells; Fig. 8C).
We use four metrics to compare the volumetric and flow-diagnostic calculations for different models. (1) Total pore volume describes the maximum potential for fluid storage. Not all of this pore volume will be available to store injected CO 2 , because some of it is unconnected (e.g. Payton et al., 2021) or holds residual water (Krevor et al., 2012). (2) Effective permeability is computed for the whole model volume in three major directions (x, y, z) using flow-based upscaling with no-flow boundaries.
(3) Pore volume injected at breakthrough time provides a measure of how much injected fluid is retained in the model volume as a result of stratigraphic baffling and trapping. Since flow diagnostics are calculated for tracer flow, values of pore volume injected at breakthrough time do not account for relative permeability or residual water saturation; calculated values are thus likely to be systematic overestimates, but are appropriate as indicative values for screening purposes. (4) The Lorenz coefficient describes the degree of heterogeneity under dynamic conditions within the storage unit, and varies between 0 for a homogeneous unit and 1 for a completely heterogeneous unit (Schmalz and Rahme, 1950). Pore volume injected at breakthrough time and the Lorenz coefficient are calculated for both west-to-east (coloured red in Fig. 8A, B) and south-to-north displacements (coloured blue in Fig. 8A, B).

Stratigraphic architectures
Representation of stratigraphic architecture is assessed by visual inspection of the models (e.g. Fig. 9). Stratigraphic architecture reflects a combination of four heterogeneities that were sketched explicitly in the models: (1) planform geometry; (2) clinoform dip; (3) lateral continuity of carbonate-cemented concretions along transgressive surfaces and maximum flooding surfaces; and (4) dip extent of facies-association interfingering down clinoforms (Table 2). Visual inspection shows that the sketched models (e.g. Fig. 9) match closely with the geological concepts and data on which they are based (Figs. 2-5A-D). Thus, deterministic scenarios are represented faithfully in sketch-based models.

Total pore volume
The mean value of total pore volume in the suite of 32 models is 268 × 10 6 m 3 . There is little variation between models, with a minimum value of 259 × 10 6 m 3 in models number 5 and 21 (Table 3) and a maximum value of 277 × 10 6 m 3 in models number 4 and 20 (Table 3). Changes in total pore volume attributed to changes of the setting for each heterogeneity (Table 2) are small (<2.5%; Fig. 10). Clinoform dip and planform geometry are the heterogeneities with the largest impact on total pore volume (>1%; Fig. 10), because they change the areal extent and thickness of high-porosity, sandstone-rich facies associations (FA3, FA4, FA5). For example, the planform geometry of southwarddeflected, elongate deltas (e.g. Fig. 9B, C) incorporates a larger areal extent of FA3 than that of westward-prograding, arcuate deltas (e.g. Fig. 9A). Models containing steep clinoforms have greater crosssectional areas of high-porosity facies associations FA3-FA5 than models containing gently dipping clinoforms (e.g. compare models with steep clinoform dips in Fig. 9B and gentle clinoform dips in Fig. 9A, C), which increases total pore volume (Fig. 10). Only clinoform dip generates changes in total pore volume that are larger than those for the confounded interactions of two heterogeneities (Fig. 10).

Effective permeability
Effective permeability at the scale of the model volume is highly anisotropic (Fig. 11). Mean values of effective permeability in the suite of models are 15 mD, 220 mD and 0.076 mD in the north-south (k x ), west-east (k y ) and vertical (k z ) orientations, respectively. Planform geometry has the greatest impact on k x and k y (Fig. 11A). Changing the planform geometry from a westward-prograding, arcuate delta (e.g. Fig. 9A) to a southward-deflected, elongate delta (e.g. Fig. 9B) introduces an additional flow path through connected, highpermeability facies associations FA3-FA5 in the western part of the model, which increases k x by 137% on average (Fig. 11A), and also a region of low-permeability facies associations FA1-FA2 in the eastern part of the model, which decreases k y by 70% on average (Fig. 11A). The corresponding increase in k z is 46% on average (Fig. 11B), which is attributed to the increased areal extent of high-permeability facies associations FA3-FA5 in a southward-deflected, elongate planform geometry (e.g. Fig. 8B, C).  (Table 3, Fig. 7) showing well placements for south-to-north tracer flow (blue) and west-to-east tracer flow (red). The model is shown for coarse (A) and fine (B) grid resolutions. (C) Plot of simulated south-to-north breakthrough time in model 32 for different grid resolutions. The selected grid resolution (labelled "c" in red) is the coarsest for which breakthrough time is considered to be reasonably accurate.
Using the geometric mean for the permeability of each facies association, as a proxy for increasing bioturbation intensity, decreases effective permeability in all three orientations, by 70% (k x ), 48% (k y ) and 44% (k z ) (Fig. 11). The horizontal permeability of each facies association (k h ) is decreased by changing from the arithmetic to geometric mean (Table 1), which in turn results in decreased k x and k y (Fig. 11A). The vertical permeability (k h ) of heterolithic cross-bedded sandstones (FA4, FA5) is also decreased by using the geometric mean in models with the low setting of mudstone drape continuity and extent in these sandstones (Table 2, Fig. 11B).
The dip extent of facies-association interfingering has a significant effect on k x , which is increased by the change from small to large downdip interfingering distances by 90% on average (Fig. 11A). The corresponding increase in k y is smaller (9% on average; Fig. 11A). Increasing the extent of facies-association interfingering increases the area of connected, high-permeability facies associations FA3-FA5 that are in contact between opposite faces of the model, thus widening flow paths through connected, highly permeable facies associations FA3-FA5 and increasing effective horizontal permeability (Fig. 11A). This effect is most pronounced in the north-south orientation, hence the change in k x is greater. Increasing the extent of facies-association interfingering decreases k z by 31% on average (Fig. 11B), because fingers of lowpermeability facies associations FA1-FA2 become more laterally extensive as well as fingers of high-permeability facies associations FA3-FA5 (e.g. compare models with large facies-association interfingering extent in Fig. 9A and small facies-association interfingering extent in Fig. 9B, C).
The change from gentle to steep clinoform dip increases k y and k z by 51% and 28% on average, respectively, but decreases k x by 15% on average (Fig. 11). Models containing steep clinoforms have a greater cross-sectional areas of highly permeable facies associations FA3-FA5 at their western face, particularly in models with a southward-deflected, elongate planform geometry (e.g. compare models with steep clinoform dips in Fig. 9B and gentle clinoform dips in Fig. 9A, C), which increases k y and k z (Fig. 11). k z is decreased by increasing the lateral continuity of carbonatecemented concretions along transgressive surfaces and maximum flooding surfaces (by 57% on average), by increasing the mean lateral extent of carbonate-cemented concretions in between transgressive surfaces (by 23% on average) and by increasing the mean vertical spacing of carbonate-cemented concretions in between transgressive surfaces (by 13% on average) (Fig. 11B). These results are expected given the method by which we calculate k v /k h ratio (Eq. (1); Begg and King, 1985). These heterogeneities have the three smallest impacts on horizontal permeability (i.e. mean of k x and k y ; Fig. 11A), although still appear to be influential (e.g. increasing the mean vertical spacing of carbonate-cemented concretions in between transgressive surfaces results in a 34% decrease in k x on average; Fig. 11A). There is no apparent reason why the lateral continuity, lateral extent and vertical spacing of carbonate-cemented concretions should influence effective horizontal permeability (e.g. Eq. (1); Begg and King, 1985). Consequently, we tentatively attribute these relatively small impacts to higher order effects between three or more heterogeneities that describe carbonatecemented concretion distributions confounded with the interfingering, areal extents and pinch-out geometries of facies associations.
Decreasing k v /k h ratio as a proxy for increased mudstone drape continuity and extent in heterolithic cross-bedded sandstones (FA4, FA5) has a similarly small influence on k x (average decrease of 37%) and k y (average increase of 16%) (Fig. 11A), which is also attributed to higher order effects between three or more confounded heterogeneities. The corresponding influence on k z is 0% (Fig. 11B).
Our interpretation of these results implies that higher order effects between three or more confounded heterogeneities are significant for k x and k y . This inference is supported by the results presented in Fig. 11A, since only planform geometry, bioturbation intensity and faciesassociation interfingering extent generate changes in k x and k y that are larger than those for the confounded interactions of two heterogeneities (Fig. 11A). The effects of multiple confounded heterogeneities on k z are less significant (Fig. 11B).

Pore volume injected at breakthrough time
Mean values of pore volume injected at breakthrough time are 0.35, with a range of 0.20 (model number 14; Table 3) to 0.59 (model number 17; Table 3), for south-to-north displacements (coloured blue in Fig. 9) and 0.52, with a range of 0.18 (model number 10; Table 3) to 0.76 (model number 30; Table 3), for west-to-east displacements (coloured red in Fig. 9) in the suite of models. There is therefore a four-fold variation in the volume of injected fluid retained in the model volume as a result of stratigraphic baffling and trapping at simulated breakthrough time.  (Table 2): (A) model number 9, characterised by an arcuate planform geometry, gentle clinoform dip, and large facies-association interfingering extent (Table 3); (B) model number 4, characterised by an elongate planform geometry, steep clinoform dip, and small facies-association interfingering extent (Table 3); and (C) model number 2, characterised by an elongate planform geometry, gentle clinoform dip, and small facies-association interfingering extent (Table 3). See Fig. 3 for key to facies association colours.
Clinoform dip and the dip extent of facies-association interfingering down clinoforms are the two heterogeneities with the largest impact on pore volume injected at breakthrough time (Fig. 12). Increasing faciesassociation interfingering extent decreases the pore volume injected at breakthrough time for both south-to-north and west-to-east displacements (by 34% and 21% on average, respectively; Fig. 12), because it increases the area and lateral connectivity of high-permeability facies associations FA3-FA5 that are in contact between injection and offtake wells, thus increasing k x and k y (Fig. 11A) and decreasing breakthrough time. Increasing clinoform dip increases the cross-sectional areas of sandstone-rich facies associations FA3-FA5 (e.g. compare models with steep clinoform dips in Fig. 9B and gentle clinoform dips in Fig. 9A, C) and thus total pore volume (Fig. 10), while also increasing k y and k z (Fig. 11). Consequently, pore volume injected at breakthrough time is decreased for west-to-east displacements (by 26% on average; Fig. 12) but increased for south-to-north displacements (by 13% on average; Fig. 12). Changing the planform geometry from a westward-prograding, arcuate delta (e.g. Fig. 9A) to a southward-deflected, elongate delta (e.g. Fig. 9B, C) also increases total pore volume (Fig. 10), increases k x and k z (Fig. 11) and decreases k y (Fig. 11), which in combination result in increased pore volume injected at breakthrough time for both south-tonorth and west-to-east displacements (by 17% and 7% on average, respectively; Fig. 12). Fig. 13 uses two contrasting models to illustrate further that simulated tracer dispersal and volumetric sweep are strongly influenced by heterogeneities that control the distribution and connectivity of high-permeability facies associations FA3-FA5. Poorly connected volumes of high-permeability facies associations FA3-FA5 are not contacted by the injected tracer (e.g. in the eastern part of model number 4; Figs. 9B, 13G-L), while connected volumes of these facies associations serve to channel injected tracer towards the offtake well (e. g. south-to-north displacement along the southward-deflected, elongate planform geometry of model number 4; Figs. 9B, 13G-L).
Increasing the mean lateral extent of carbonate-cemented concretions in between transgressive surfaces decreases the pore volume injected at breakthrough time for both south-to-north and west-to-east displacements (by 17% and 12% on average, respectively; Fig. 12). The increased lateral extent of the concretions decreases k y and k z but increases k x (Fig. 11), thus promoting lateral flow from injection to offtake wells, particularly for south-to-north displacements (Fig. 12). For similar reasons, increasing the lateral continuity of carbonatecemented concretions along transgressive surfaces and maximum flooding surfaces also decreases the pore volume injected at breakthrough time for south-to-north displacements (by 13% on average; Fig. 12), although it increases slightly for west-to-east displacements (by 7% on average; Fig. 12).
Changes in mudstone drape continuity and extent in heterolithic cross-bedded sandstones, mean vertical spacing of carbonate-cemented concretions in between transgressive surfaces, and bioturbation intensity all have a relatively minor influence on the pore volume injected at breakthrough time (Fig. 12). The effects of these heterogeneities individually are smaller than those for the confounded interactions of two heterogeneities (Fig. 12).

Lorenz coefficient
Values of the Lorenz coefficient are higher for south-to-north displacements (mean of 0.71, and range of 0.59-0.78) than for west-to-east displacements (mean of 0.58, and range of 0.48-0.72) in the suite of models. The models are therefore more heterogeneous in the northsouth orientation than the west-east orientation.
The dip extent of facies-association interfingering down clinoforms has the greatest impact on the Lorenz coefficient, but planform geometry and clinoform dip are both also influential (Fig. 14). Increasing faciesassociation interfingering extent increases the Lorenz coefficient in the  Fig. 10.. Tornado chart showing the average percentage changes in total pore volume that result from varying each factor from its low setting to its high setting (Table 2). If the bar lies to the right then the change is positive. For example, modelling steep clinoform dips (high setting) increases total pore volume by c. 2% compared with modelling gentle clinoform dips (low setting). The largest response of confounded 2-factor interactions is shown for comparison with the main effects due to individual factors. Changes in total pore volume are small (<2.5%) for all factors. Effective horizontal permeability in model volume: kx (North-South) ky (West-East) (planform geometry + dip extent of facies interfingering along clinoforms) + (clinoform dip + mud-drape continuity and extent in heterolithic cross-bedded sandstones) mean value for modelled scenarios = 15 mD mean value for modelled scenarios = 220 mD A B Fig. 11.. Tornado charts showing the average percentage changes in (A) effective horizontal permeability in north-south (k x ) and west-east (k y ) orientations, and (B) effective vertical permeability (k z ) that result from varying each factor from its low setting to its high setting (Table 2). If the bar lies to the right then the change is positive. For each tornado chart, the largest response of confounded 2-factor interactions is shown for comparison with the main effects due to individual factors. Three and five individual factors have main effects greater than the largest response of confounded 2-factor interactions for effective horizontal permeability and effective vertical permeability, respectively. west-east orientation (by 14% on average; Fig. 14), but decreases it in the north-south orientation (by 8% on average; Fig. 14). Changing the planform geometry from a westward-prograding, arcuate delta (e.g. Fig. 9A) to a southward-deflected, elongate delta (e.g. Fig. 9B) increases the Lorenz coefficient in the north-south orientation (by 7% on average; Fig. 14) and decreases it in the west-east orientation (by 5% on average; Fig. 14). Increasing clinoform dip increases the Lorenz coefficient in both north-south and west-east orientations (by 3% and 8% on average, respectively; Fig. 14). Facies-association interfingering extent, planform geometry and clinoform dip all control facies-association distribution and connectivity (e.g. Figs. 9, 13), thus changing these parameters affects values of mean k x and k y in successive layers of the storage unit (which in this study are represented by successive layers of the orthogonal model grid). The Lorenz coefficient reflects the cumulative frequency distributions of layer-specific mean k x and layer-specific mean k y , and is thus modified accordingly. Facies-association interfingering extent, planform geometry and clinoform dip also exert significant influence on effective horizontal permeability at the scale of the entire model volume (Fig. 11A), although the rank order of these heterogeneities for effective horizontal permeability differs to that for the Lorenz coefficient (Fig. 14).
Changing from the arithmetic mean (k h ) and k v /k h modifiers to the geometric mean for the permeability of each facies association, as a proxy for increasing bioturbation intensity, increases the Lorenz coefficient in both north-south and west-east orientations (by 5% and 10% on average, respectively; Fig. 14). The effects of other individual heterogeneities, which describe the distribution of carbonate-cemented concretions in sandstones (FA3-FA5) and the continuity and extent of mudstone drapes in heterolithic cross-bedded sandstones (FA4, FA5), are smaller than those for the confounded interactions of two heterogeneities (Fig. 14). These results are similar to those for effective horizontal permeability (k x , k y ) at the scale of the entire model volume (Fig. 11A). Other heterogeneities, which impact k z (e.g. lateral continuity of carbonate-cemented concretionary layers along transgressive surfaces; Fig. 11B), exert little influence on the cumulative frequency distributions of layer-specific mean k x and layer-specific mean k y that are used to calculate the Lorenz coefficient.

Which heterogeneities control CO 2 migration and stratigraphicbaffling potential?
Our results indicate that for the heterogeneities and settings chosen ( Table 2), most of the heterogeneities investigated have a significant impact on effective permeability (k x , k y , k z ), pore volumes injected at breakthrough time, and/or the Lorenz coefficient (Figs. 11,12,14). Heterogeneities that control the distribution and connectivity of highpermeability facies associations FA3-FA5 (i.e. planform geometry, clinoform dip, facies-association interfingering extent along clinoforms) impact all of these flow-diagnostic measures. Heterogeneities that control the permeability characteristics of stratigraphic surfaces (i.e. lateral continuity of carbonate-cemented concretionary layers along transgressive surfaces) strongly influence k z (Fig. 11B). Of the heterogeneities that control the internal permeability characteristics of highpermeability facies associations FA3-FA5, only bioturbation intensity significantly impacts k x , k y , k z and the Lorenz coefficient (Figs. 11,14), with other heterogeneities (i.e. vertical spacing and lateral extent of carbonate-cemented concretions, continuity and extent of mudstone drapes) appearing to play a minor role. The interactions of two factors are also influential for all of the studied flow-diagnostic measures (Figs. 11,12,14).
These results are consistent with previous work evaluating the role of heterogeneities on sweep and recovery from shallow-marine sandstones  (Table 2). If the bar lies to the right then the change is positive. The largest response of confounded 2-factor interactions is shown for comparison with the main effects due to individual factors. Three individual factors have main effects greater than the largest response of confounded 2-factor interactions.
that act as hydrocarbon reservoirs. It is widely appreciated that the choice of depositional model (Fig. 2B, C) exerts a fundamental control on the geometry, extent and distribution of high-permeability sandstones in shallow-marine settings (e.g. Howell et al., 2008a;Willis et al., 2021), while clinoform geometry and the distribution of associated heterogeneities may have significant impacts on flow tortuosity (e.g. White et al., 2004;Howell et al., 2008b;Jackson et al., 2009;Graham et al., 2015). Laterally extensive cemented horizons below transgressive surfaces have also been noted to act as barriers to vertical flow (e.g. Morad et al., 2010).
The impact of bioturbation in homogenising the texture and permeability characteristics of sandstone facies associations has been recognized previously (e.g. Gingras et al., 2012). The small apparent impact of other heterogeneities that control the internal permeability characteristics of sandstone facies associations (i.e. vertical spacing and lateral extent of carbonate-cemented concretions, continuity and extent of mudstone drapes) appears surprising, but may reflect the limited range of values assigned to parameters that describe these heterogeneities based on cores from the 31/5-7 (Eos) well. For example, mudstone drapes may be more pervasive along cross-set toesets and foresets in facies associations FA4 and FA5 than those sampled in core, and this would then increase their impact on flow tortuosity (cf. Jackson et al., 2005;Ringrose et al., 2005;Massart et al., 2016).
Assuming that the studied heterogeneities are independent of each other, then their most favourable combination to slow migration and retain injected fluid in the modelled rock volume is: (1) a southwarddeflected, elongate delta planform (e.g. Fig. 2B); (2) low continuity of carbonate-cemented concretions along transgressive surfaces and maximum flooding surfaces; (3) low values of dip extent of faciesassociation interfingering; (4) high values of mean vertical spacing of carbonate-cemented concretions in between transgressive surfaces; and (5) low continuity and extent of mudstone drapes in heterolithic crossbedded sandstones (FA4, FA5) (Fig. 12). These heterogeneity settings result in greater homogeneity in the areal extent, connectivity and  Fig. 3 for key to facies association colours.
internal permeability structure of sandstones (FA3-FA5), and thus allow relatively isotropic plume migration and dispersal. In contrast, rapid and anisotropic plume migration is favoured by: (1) a westward-prograding, arcuate delta planform (e.g. Fig. 2C); (2) high continuity of carbonatecemented concretions along transgressive surfaces and maximum flooding surfaces; (3) high values of dip extent of facies-association interfingering; (4) low values of mean vertical spacing of carbonatecemented concretions in between transgressive surfaces; and (5) high continuity and extent of mudstone drapes in heterolithic cross-bedded sandstones (FA4, FA5) (Fig. 12). Whether or not clinoform dip, mean lateral extent of carbonate-cemented concretions in between transgressive surfaces, and bioturbation intensity act to slow migration and increase retention of injected fluid depends on the orientation of fluid movement (Fig. 12).
In summary, our screening results imply that it is important to consider the effects of multiple sedimentological heterogeneities on the stratigraphic trapping potential of CO 2 in the Johansen-Cook storage unit. The effects of heterogeneities acting in combination are also important (Figs. 11, 12, 14), but have not been explored here. The baffling effect of sedimentological heterogeneity on CO 2 migration and retention is likely to be an important precursor for later capillary, dissolution and mineral trapping. Our results are consistent with the recently published modelling work of Meneguolo et al. (2022), which includes drainage and imbibition relative permeability curves for CO 2 and thus captures the effects of capillary trapping, for the Johansen and Cook formations over the entire Northern Lights CO 2 storage site. Their findings indicate that rapid areal migration of the CO 2 plume is favoured by high lateral sandstone connectivity and that vertical plume migration is aided by high k v /k h ratios across stratigraphic surfaces and within the units they bound (Meneguolo et al., 2022).

Advantages and limitations of screening methodology
The methodology presented here is fast and efficient in quickly screening the impact of sedimentological heterogeneity on flow properties. We attribute this efficiency to a combination of three elements: (1) experimental design; (2) a sketch-based modelling approach; and (3) single-phase flow diagnostics. All three elements combine to make the methodology more efficient than the sum of its parts, but we emphasise that each element, and thus the combined methodology, has its limitations.
While experimental design allows economical exploration of a wide parameter space, the range of heterogeneities and their values is defined by the geological scenarios considered. If a viable geological scenario is not considered, then the resulting parameter space may be too narrow. Examples of additional scenarios for the Johansen-Cook CO 2 storage unit include: (1) alternative depositional models to those shown in Fig. 2B, C (e.g. containing the channelised and landward-tapering, wedge-shaped depositional elements interpreted by Meneguolo et al., 2022); (2) alternative definitions of facies associations and thus parameterisations of their porosity and permeability characteristics; (3) inclusion of steeper clinoform dips; and (4) inclusion of more abundant and/or laterally continuous mudstone drapes (Fig. 6) (Section 4.2). In addition, experimental design does not allow the effects of interactions between factors (heterogeneities) to be investigated fully (Kjønsvik et al., 1994;White and Royer, 2003). In the fractional factorial design applied in this study, the effects of stochastic variation in each heterogeneity are also not considered. Nonetheless, the design is sufficient to identify the most significant of the studied heterogeneities that influence a particular response.
Our sketch-based modelling approach allows rapid construction of prototype models of sedimentological heterogeneity (Costa Jacquemyn et al., 2021a). However, the resulting models are deterministic representations of a particular interpretation or scenario, (planform geometry + dip extent of facies interfingering along clinoforms) + (clinoform dip + mud-drape continuity and extent in heterolithic cross-bedded sandstones) mean value for modelled scenarios = 0.71 mean value for modelled scenarios = 0.58 Fig. 14.. Tornado chart showing the average percentage changes in Lorenz coefficient that result from varying each factor from its low setting to its high setting (Table 2). If the bar lies to the right then the change is positive. The largest response of confounded 2-factor interactions is shown for comparison with the main effects due to individual factors. Four individual factors have main effects greater than the largest response of confounded 2-factor interactions. and do not allow stochastic realisations to be generated of this scenario. As a result, probability distributions cannot be generated for the model responses (i.e. the metrics used to compare the models). This approach is well suited to a screening study that explicitly addresses the wide range of uncertainty encapsulated in the modelled scenarios (Bentley and Smith, 2008), but is not appropriate for characterising the uncertainty inherent to a particular scenario. Sketch-based models are also limited by the user's ability to sketch their interpreted geological scenario, and by the degree to which this scenario is reasonable given the supporting dataalthough the latter is true for any reservoir model, irrespective of the method used to generate it. Sketch-based models can be constructed by users who have domain expertise in sedimentology, but are not specialists in reservoir modelling. In principle, this is an advantage in generating models of geologically reasonable scenarios.
Flow diagnostics are computationally cheap (Shahvali et al., 2012;Møyner et al., 2014), but are calculated using a single fluid phase and a steady-state pressure field for a given combination of fluid injection and offtake (production) wells. As a result, the effects of relative permeability between injected CO 2 and in-place water (and potentially hydrocarbons) are not considered in flow diagnostics. Consequently, we do not consider capillary trapping, which may result in underestimation of the absolute effects of heterogeneities as baffles and barriers. Supercritical CO 2 is compressible above its critical pressure (74 bar), although only to a relatively small degree under the conditions required for geological storage (Hassanzadeh et al., 2008), and is less dense than water, giving rise to gravity segregation where gravity forces dominate over viscous forces (Ide et al., 2007). Neither fluid compressibility nor transient flow due to gravity segregation are considered in flow diagnostics, and structural dip also needs to be included in models that consider the effects of gravity segregation. Consequently, we regard our results as indicative of flow behaviour rather than accurate predictions. More detailed modelling of combined structural-stratigraphic trapping and capillary trapping of CO 2 require full-physics, multiphase simulations. Further, dissolution and mineral trapping mechanisms that operate over medium-term (100-10,000 years) and long-term (>10,000 years) timescales also need to be modelled.
The limitations described above do not undermine the results of this screening study, but they do emphasise that our work serves as a precursor to more detailed but time-consuming modelling studies. The screening study serves to illustrate that multiple types of sedimentological heterogeneity are important in controlling CO 2 migration and storage by stratigraphic baffling and trapping, and that the effects of these heterogeneities need to be incorporated into more detailed future modelling work.

Conclusions
The saline aquifers of the Johansen and Cppook formations constitute the primary CO 2 storage unit in the Northern Lights project, but are sparsely sampled in the storage site. Sedimentological analysis of core photos from the recent 31/5-7 (Eos) confirmation well support the previously interpreted depositional model of a wave-dominated delta front that underwent three episodes of progradation, aggradation and abandonment ("lower Johansen, upper Johansen and Cook deltas"). However, uncertainty in the depositional model of the Johansen-Cook storage unit remains, particularly in regions not sampled by wells and below seismic resolution. Uncertain sedimentological heterogeneities include: delta planform geometry; clinoform dip; the interfingering extent of facies associations down clinoforms; the lateral continuity, lateral extent and vertical spacing of carbonate-cemented concretions; mudstone-drape continuity and extent in cross-bedded, proximal deltafront sandstones; and bioturbation intensity. A method combining experimental design, sketch-based reservoir modelling, and single-phase flow diagnostics was used to rapidly and efficiently screen the impact of these sedimentological heterogeneities on simulated CO 2 migration and retention by stratigraphic baffling and trapping in the absence of capillary effects and gravity segregation. Data from outcrop and subsurface depositional analogues were used to constrain values of the investigated heterogeneities, in addition to data from the Johansen-Cook storage unit. Models are 2.5 km in south-north extent, 3.0 km in west-east extent and 200 m thick, and lack faults and tectonic dip in order to isolate the effects of sedimentological heterogeneity. Flow was simulated in south-to-north and west-to-east orientations between a single injection and a single offtake well.
Our results indicate that heterogeneities which control the distribution and connectivity of high-permeability medial and proximal deltafront sandstones (i.e. delta planform geometry, clinoform dip, and facies-association interfingering extent along clinoforms) significantly impact effective horizontal and vertical permeability, the Lorenz coefficient, and pore volumes injected at breakthrough time. In addition, the lateral continuity of carbonate-cemented concretionary layers along transgressive surfaces strongly influences effective vertical permeability, and bioturbation intensity significantly impacts effective horizontal and vertical permeability and the Lorenz coefficient. Heterogeneities acting in combination are also influential. Thus, the baffling effect on CO 2 migration and retention of sedimentological heterogeneity is significant, and provides a precursor template for capillary, dissolution and mineral trapping.
The method used in this screening assessment has the potential to be used in many other evaluations of subsurface geological heterogeneity and its impact on flow. Experimental design allows efficient exploration of a wide parameter space, the sketch-based modelling approach enables rapid construction of deterministic models of interpreted scenarios, and single-phase flow diagnostics provide computationally cheap approximations of full-physics, multiphase simulations that are reasonable for many subsurface-flow conditions. In combination, these three elements result in a method that is more efficient than the sum of its parts.

Declaration of Competing Interest
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 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 funding this work and granting permission to publish this paper. Geiger acknowledges partial funding for his Chair from Energi Simulation. We also thank Jafar Alshakri for practical discussion of the sketch-based models presented herein. WAJ constructed and analysed the models as part of the Petroleum Geoscience MSc course at Imperial College London.