Abstract

Geological modelling is widely used to predict resource potential in subsurface reservoirs. However, modelling is often slow, requires use of mathematical methods that are unfamiliar to many geoscientists, and is implemented in expert software. We demonstrate here an alternative approach using sketch-based interface and modelling, which allows rapid creation of complex three-dimensional (3D) models from 2D sketches. Sketches, either on vertical cross-sections or in map-view, are converted to 3D surfaces that outline geological interpretations. We propose a suite of geological operators that handle interactions between the surfaces to form a geologically realistic 3D model. These operators deliver the flexibility to sketch a geological model in any order and provide an intuitive framework for geoscientists to rapidly create 3D models. Two case studies are presented, demonstrating scenarios in which different approaches to model sketching are used depending on the geological setting and available data. These case studies show the strengths of sketching with geological operators. Sketched 3D models can be queried visually or quantitatively to provide insights into heterogeneity distribution, facies connectivity or dynamic model behaviour; this information cannot be obtained by sketching in 2D or on paper.

Supplementary material: Rapid Reservoir Modelling prototype (executable and source code) is available at: https://bitbucket.org/rapidreservoirmodelling/rrm. Supplementary screen recordings for the different case studies showing sketch-based modelling in action are available at https://doi.org/10.6084/m9.figshare.c.5084141 and supplementary figure S1-S4 are available at https://doi.org/10.6084/m9.figshare.c.5303043

Capturing geological heterogeneity is a key goal when constructing numerical models of the Earth's subsurface (e.g. Denver and Phillips 1990; Boggs et al. 1992; Hamilton and Jones 1992; Koltermann and Gorelick 1996; Miller et al. 1998; White et al. 2004; Jackson et al. 2005, 2013a, 2015a; Janković et al. 2006; Ronayne et al. 2008, 2010; Linde et al. 2015; Bentley 2016). Such models are used to predict the resource distribution and extraction potential of groundwater resources, geothermal reservoirs, oil and gas reservoirs, and ore deposits, as well as the behaviour of subsurface targets for nuclear waste or CO2 storage (O'Sullivan et al. 2001; Matthäi et al. 2007a; Sun et al. 2008; Geiger et al. 2009; Sech et al. 2009; Ingebritsen et al. 2010; Deng et al. 2012; Shamshiri and Jafarpour 2012; Sommer et al. 2013; Gershenzon et al. 2015; Graham et al. 2018; March et al. 2018). However, the three-dimensional (3D) geometry and spatial distribution of geological heterogeneity in the subsurface is uncertain, as boreholes sample only a small fraction of the rock volume and geophysical imaging methods lack the spatial resolution required to delineate all heterogeneities of interest (Lemon and Jones 2003; Feyen and Caers 2006; Ronayne et al. 2010; Linde et al. 2015). Expert knowledge is often used to predict the likely geometry and spatial distribution of heterogeneity away from boreholes, but geological interpretations are themselves subject to uncertainty. Numerical models based on different geological interpretations can yield highly varying predictions of system behaviour (e.g. Bond et al. 2007; Bentley and Smith 2008; Refsgaard et al. 2012; Deveugle et al. 2014; Brunetti et al. 2019).

Even though geological models of the subsurface are very commonly used, current modelling tools often require specialist knowledge to capture the features of interest and are slow to create or update. Therefore they are not intuitive to use and are inaccessible to non-specialists or different disciplines. This is not ideal to quickly test an hypothesis or explore different geological concepts before starting a detailed modelling workflow. Here we will present a sketch-based approach that is intuitive to sketch geological models in 3D and is meant to rapidly capture geological ideas and concepts in 3D models. This enables us to communicate, learn, quantify or rank geology and its impact on flow in 3D.

Approaches to computer-aided design (CAD) and computational fluid dynamics (CFD) outside of the geological modelling domain often include a significant element of prototyping: a number of simplified models are created to test different design concepts, before detailed models are created of the final agreed design (e.g. Shah et al. 2001; Cherlin et al. 2005; Arisoy and Kara 2014). In geological modelling, analysis of different geological concepts is sometimes termed scenario testing, to differentiate it from the probabilistic modelling approaches that are used to create numerous model realizations around a single geological concept or scenario (Bentley and Smith 2008). However, prototyping in geological modelling is rare, in part because there are no methods or software tools that allow rapid creation and testing of geological concepts. Most previous studies have instead focused on the development of methods and tools for probabilistic modelling around a single concept (e.g. MacDonald and Aasen 1994; Yao et al. 2005; Strebelle 2006; Linde et al. 2015; Brunetti et al. 2019).

Sketch-based interface and modelling (SBIM) is an approach for rapid model creation that is used for prototyping in many non-geological CAD and CFD applications (e.g. Olsen et al. 2009, 2011; Pereira et al. 2011). SBIM allows model concepts to be rapidly sketched and tested. A number of studies have also proposed the use of SBIM in geological modelling (e.g. Amorim et al. 2012, 2014; Lidal et al. 2013; Natali et al. 2014a, b; Jackson et al. 2015a; Costa Sousa et al. 2020). SBIM in this context is based on the concept of using surfaces to capture geological architecture and heterogeneity (e.g. Denver and Phillips 1990; Hamilton and Jones 1992; MacDonald et al. 1998; Willis and White 2000; White et al. 2004; Huysmans and Dassargues 2011). The surfaces define and bound geological domains (Pyrcz et al. 2005; Caumon et al. 2009; Zhang et al. 2009; Jackson et al. 2013b; Ruiu et al. 2016; Jacquemyn et al. 2019). Surfaces may represent faults, fractures, stratigraphic surfaces, facies boundaries, lithological boundaries, diagenetic boundaries or any other type of geological boundary (e.g. Willis and White 2000; White et al. 2004; Pyrcz et al. 2005; Matthai et al. 2007b; Caumon et al. 2009; Geiger et al. 2009; Sech et al. 2009; Zhang et al. 2009; Graham et al. 2015; Jackson et al. 2015b; Massart et al. 2016a, b; Jacquemyn et al. 2019). Geological domains, bounded by surfaces, can be defined without reference to an underlying grid or mesh (e.g. Jacquemyn et al. 2019), although a mesh may be created that conforms to the modelled surfaces when a numerical calculation is required (Jackson et al. 2015a; Zhang et al. 2018). The emphasis on surfaces and surface-bounded domains closely matches how geoscientists conceptualize and represent geological interpretations in traditional tools such as maps, cross-sections and block diagrams. The aim of SBIM here is not to simulate underlying geological processes, but to capture their effects on resulting geometries.

Previous applications of SBIM for geological modelling have been limited in their scope and applicability. Natali et al. (2014a, b) applied a SBIM methodology to sedimentary geology, with the goal of creating figures or animations for teaching or discussion. In their approach, geology must be sketched in depositional order, thus requiring that a geological history be interpreted prior to sketching. This is too restrictive for use as a prototyping SBIM tool, because it does not allow surfaces to be sketched in any order such that different concepts can be tested. Amorim et al. (2012) focused on interpretation of horizons from 3D seismic data using sketched surfaces. In contrast, Amorim et al. (2014) started with a blank screen and allowed the user to sketch lines, representing stratigraphic boundaries or faults, in map-view. The user sketched geological symbols onto the map which were used to control extrapolation of the sketched lines into 3D to form surfaces. Although both approaches constitute advanced examples of SBIM for geological application, they were too limited to be used for realistic geological modelling because, in their worked examples, the sketched surfaces did not interact, and no operations were outlined for how interacting surfaces should be treated. Surface intersections in natural systems are ubiquitous, and a practical SBIM methodology for geology must have a robust and geologically intuitive approach to handle intersections. Furthermore, sketching only in map-view is too limited to replicate the geological architectures observed in natural systems. Geoscientists use sketches in both cross-section and map-view to capture these.

Here we report, for the first time, the use of SBIM to create 3D models of complex stratigraphy and structure, with direct relevance to rapid prototyping of subsurface reservoirs. We demonstrate application of our new SBIM approach using two typical example reservoir modelling scenarios, in which sketches are made at different scales, of different geological heterogeneity, and are constrained by different input data.

The SBIM methods presented here are implemented in an open source research code (Rapid Reservoir Modelling, RRM; Fig. 1; https://bitbucket.org/rapidreservoirmodelling/rrm), developed collectively by the authors, which links a sketching interface with geological operators to build 3D models. The code also includes a flow diagnostics module to provide direct quantitative evaluation of such models. The aim of this paper is not to describe the numerical implementation in detail; a summary is provided, and other publications elsewhere describe the numerical algorithms and underpinning SBIM design concepts (Jackson et al. 2015a). Rather, we focus on how SBIM has been adapted and implemented for geological modelling, and on demonstrating its utility in rapidly testing different interpretations for subsurface reservoirs.

We first describe how 3D surfaces are created from 2D sketches and outline a suite of operators that control surface interactions such that the resulting sketched models are geologically realistic. In this first demonstration of SBIM for geological modelling, we explore its flexibility in stratigraphic and structural case studies. Both case studies illustrate scenarios that are common and important for subsurface reservoir characterization.

Methods

Our geological SBIM framework consists of two main components: a sketching component, which enables sketching on 2D planes from which 3D surfaces are generated, and a surface interaction component, which handles the intersections between the 3D surfaces and guarantees the representational validity of the sketched model. Screen recordings of the real-time sketching process and surface interactions for the models in this paper (Case studies I and II) are available in the supplementary materials.

Creating 3D surfaces from 2D sketches

Geoscientists commonly create 2D sketches using pencil and paper. SBIM adopts the same approach but allows the geoscientist to create 2D sketches in a digital environment using an input device, such as a mouse or stylus, to draw directly on screen. Sketches are made in cross-section and/or map-view; however, the objective of SBIM as presented here is to create a 3D model. We have implemented three different approaches to create 3D surfaces from 2D sketches.

  1. Sketch the intersection of the geological surface with one or multiple vertical 2D cross-sections and create a 3D surface by interpolation between the sketches.

  2. Sketch depth and elevation contours on multiple horizontal 2D planes in map-view and create a 3D surface by interpolation between the sketches.

  3. Sketch the intersection of the geological surface with a vertical 2D cross-section, and then sketch a trajectory in map-view. The cross-section is extruded along the trajectory to create a 3D surface.

Interactions of 3D surfaces

A key innovation in our geological SBIM approach is the definition of operators that dictate how sketched surfaces interact in 3D. These operators are based on fundamental geological rules for the interaction of stratigraphic surfaces that ensure any geological model is valid (White and Barton 1999; Jackson et al. 2005; Caumon et al. 2009).

  1. Surfaces cannot cross.

  2. Surfaces cannot end within a domain.

  3. Surfaces can either terminate against (truncate or conform) or remove (erode), existing surfaces.

The operators are designed to capture the relationships between stratigraphic surfaces that are required for model construction, but not to mimic geological processes. This is important because it allows the users to ‘sketch what they see’ without having to reproduce the underlying processes that created the observed stratigraphy. Moreover, it ensures that the operators are universally applicable to all stratigraphic settings and scales. The selected operator(s) are applied immediately after each new 3D surface is created, which means that surfaces do not need to be created in stratigraphic or hierarchical order. This is important for the prototyping applications we propose here, in which a complete geological interpretation may not be available at the onset of sketching; indeed, sketching and creation of 3D models may be a part of the interpretation process. The operators ensure that models contain watertight volumes (cf. Caumon et al. 2004) such that no gaps exist between truncating surfaces. We have identified two sets of operators, outlined below; visual representations of the operators are available in supplementary material S1 to S4, and screen recordings showing the use of operators in action are available in the supplementary materials.
  1. Operators that define how new sketched surfaces modify existing surfaces. The new sketched surface crops all or some of the existing surfaces.

        1.A ‘Remove Above (RA)’ – the RA operator removes the parts of any existing surfaces that lie above the new sketched surface. Existing surfaces that lie below the new sketched surface remain unchanged. A simple practical example of when this operator could be used is when sketching an erosional surface.

        1.B ‘Remove Above Intersection (RAI)’ – the RAI operator removes the parts of existing surfaces that are intersected by and lie above the new sketched surface. Existing surfaces that are not intersected by the new sketched surface remain unchanged. A simple practical example of when this operator could be used is when sketching an erosional surface after sketching overlying deposits.

        1.C ‘Remove Below (RB)’ – the RB operator removes the parts of any existing surfaces that lie below the new sketched surface. Existing surfaces that lie above the new sketched surface remain unchanged. A simple practical example of when this operator could be used is when sketching basement contacts.

        1.D ‘Remove Below Intersection (RBI)’ – the RBI operator removes the parts of existing surfaces that are intersected by and lie below the new sketched surface. Existing surfaces that are not intersected by the new sketched surface remain unchanged. A simple practical example of when this operator could be used is when sketching lateral accretion surfaces in a channel fill while removing previously interpreted channel-fill deposits.

  2. Operators that set existing surfaces as boundaries for new sketched surfaces. The new sketched surface is cropped outside the assigned boundaries.

        2.A ‘Preserve Above (PA)’ – The PA operator is applied to any existing surface such that it acts as a lower boundary for the next sketched surface. Only the parts of the new sketched surface that lie above this lower boundary are preserved and interact with existing surfaces. A simple practical example of when this operator could be used is when sketching strata that onlap or downlap onto an existing surface.

        2.B ‘Preserve Below (PB)’ – The PB operator is applied to any existing surface, such that it acts as an upper boundary for the next sketched surface. Only the parts of the new sketched surface that lie below this upper boundary are preserved and interact with existing surfaces. A simple practical example of when this operator could be used is when sketching strata that underlie a previously sketched unconformity surface.

        2.C ‘Preserve Between (PBW)’ – The PBW operator combines the PA and PB operators. Only the parts of a new sketched surface that lie above a selected lower boundary and below a selected upper boundary are preserved and interact with existing surfaces. A simple practical example of when this operator could be used is when sketching stratigraphic surfaces within a channel fill by preserving surfaces sketched between the channel-fill top and base, or adding further heterogeneity within layers and regions.

Operators can be combined to provide the flexibility to sketch a model in any order, depending on data type, updating of data, and user preference: different users may favour a different order of sketching or interpretation. For example, some users may prefer to sketch hierarchically, starting with the large-scale features and progressively adding smaller-scale features. Some users may prefer to sketch in stratigraphic order, while others may sketch in no set order. All these approaches are possible in our SBIM implementation; irrespective of sketching order, the choice of appropriate operator allows the user to add or modify previously sketched surfaces, facilitating interpretation while sketching.

Stratigraphic surfaces are typically non-multivalued surfaces (ordinary or single-value surfaces), such that they do not recur in a vertical column unless deformed by later folding and faulting. Hence, the concepts of ‘above’ and ‘below’ can be unambiguously defined. This is not the case for multivalued surfaces resulting from deformation. Moreover, stratigraphic surfaces must be continuous across a model unless they terminate against another surface, whereas faults can end within a domain. Representing the full effects of faulting and folding is beyond the scope of this paper and is the subject of ongoing research. However, folding and faulting that does not result in multivalued surfaces can nonetheless be sketched using the current operators, so long as the fault surfaces do not end within the 3D model volume. We demonstrate preliminary application of SBIM for faulted reservoirs in Case study II.

Numerical implementation

The details of the numerical implementation of the 2D sketching, 3D surface creation and sketch operators outlined above in our SBIM framework are beyond the scope of this geoscience-focused paper; instead, we provide a brief summary.

The 2D sketches are created using a sketching engine that implements well-established SBIM tools (see Olsen et al. 2009 for a review). The engine also implements some bespoke modifications required for our geological modelling application, including (i) vertical exaggeration while sketching, to account for the very high aspect ratio of most geological models, and (ii) the ability to sketch over data in a given sketch plane, such as borehole data, seismic data or outcrop images, as demonstrated in our example applications.

The sketches are converted to a set of points defined in Cartesian space. The point data, which may have been created from multiple sketches on different cross-sections, are then passed to the surface modelling engine which creates a 3D surface from the data in one of two ways. In the first approach, the surface is created by interpolation using thin-plate splines (Duchon 1977) as kernels for surface reconstruction. Overfitting is avoided by reconstructing surfaces using approximate interpolation (Wendland and Rieger 2005), which finds the surfaces that best fit the sketch data within a prescribed margin of error. The balance between accuracy and stability in surface reconstruction is mediated by the resolution of the model being built and the complexity of the 3D surface. As with any interpolation-based surface reconstruction, the resulting surfaces may not preserve geological realism away from the sketched data. Further sketches can be added to constrain the surface reconstruction if needed.

The second approach can better preserve geological realism by extruding a cross-section sketch along a trajectory sketched in map-view. The extrusion method depends on the geometry of the cross-section and trajectory. If the curvature of the guiding trajectory is small enough that a cross-section can be extruded along it without the resulting surface intersecting with itself, then an actual extrusion is computed by interpolating the cross-section and trajectory independently and creating a surface given as a tensorial product of the input curves. To account for the more general case in which a cross-section extruded along a trajectory may self-intersect, the use of a tensorial product is excluded. In this case, the cross-section is guided using a vector field that coincides with the trajectory's tangent vectors at the trajectory location. To minimize distortion it is required that the vector field be divergence-free and that the vectors’ magnitudes near the trajectory be close to one; these restrictions are enforced by using matrix-valued conditionally positive definite functions as kernels for the vector-field interpolation (Narcowich and Ward 1994). A vector field is required that must be integrated while ensuring it is finely sampled where the guided surface has high gradient, and this more general approach for extrusion is thus more CPU intensive than using the tensorial product of cross-section and trajectory.

Once a new 3D surface is created, it is passed to the operator engine. Consistent and efficient implementation of the operators described in the previous section has been one of the most challenging aspects of our numerical implementation of SBIM. In essence, the operators are mapped onto an order-theoretic lattice defined on continuous surfaces to represent the geological models. The operator engine continuously determines surface intersections and performs additional queries, such as ordering the surfaces according to the relative geological age defined by surface interactions. Efficient implementation means that the computational cost of the operator engine scales linearly with the number of surfaces in a model.

Results

Case study I: shallow-marine strata – well correlation

Context and geological setting

Reservoirs deposited in shallow-marine environments are important groundwater resources in many areas of the world (e.g. Taylor et al. 2001; Mayo et al. 2003; Huysmans et al. 2008; Fahad Al-Ajmi et al. 2015) as well as oil and gas reservoirs and targets for geological CO2 storage (e.g. Jackson et al. 2009; Sech et al. 2009; Ashraf et al. 2013). However, stratigraphic sections can be correlated between boreholes in a variety of ways, depending on the number and spacing of boreholes, the data types available to constrain the correlation, and the depositional concept(s) applied. The aim of the modelling in this case study is to test the impact of different borehole correlations and depositional concepts on 3D facies architecture in marginal-marine and shallow-marine strata of the Spring Canyon Member, Blackhawk Formation in the Book Cliffs, Utah, USA. The Spring Canyon Member is interpreted as a series of vertically stacked, regressive–transgressive tongues (parasequences) that were deposited by wave-dominated deltas and laterally adjacent strandplains during regression and wave-dominated, tide-influenced barrier island systems (containing tidal inlets and flood tidal deltas) during transgression (Kamola and Van Wagoner 1995; Hampson and Howell 2005; Campion et al. 2010).

We use SBIM to create the 3D models, honouring data from five measured outcrop sections and boreholes (Fig. 1) (Kamola and Van Wagoner 1995; Campion et al. 2010). These data and associated outcrops are widely used in education and training to develop geological interpretation skills, but their educational value can be significantly enhanced by translating the correlation panel between measured outcrop sections and boreholes into 3D models that demonstrate the impact of different interpretations on predicted aquifer or reservoir character and behaviour. SBIM allows this in a rapid and intuitive manner and could be executed in the field using tablets equipped with touch-sensitive screens. The data used here were vertically exaggerated by a factor of 150 and we chose to sketch in this aspect ratio, but the resulting models can easily be rescaled to true aspect ratio for quantitative calculations.

Variations exist in how interpreted sequence stratigraphic surfaces, facies belts and geobodies can be correlated between the boreholes, and in the map-view geometry of the surfaces, belts and geobodies. The following five interpretation scenarios are explored. These scenarios represent different interpretations of stratigraphic architecture in regions that lie between the outcrop sections and boreholes in the 2D cross-section, and variations in 3D stratigraphic architecture away from the cross-section. The scenarios tested are as follows.

  • Scenario SM1: initial interpretation

  • Scenario SM2: variation of down-dip extent of shoreface facies belts (as recorded by sketched cross-section relationships)

  • Scenario SM3: variation of shoreline orientation and rugosity (as recorded by sketched map-view trajectory)

  • Scenario SM4: variation of the extent and geometry of flood tidal delta deposits (as recorded by sketched cross-section relationships and map-view trajectory)

  • Scenario SM5: variation of the extent and geometry of tide-influenced channelized sandbodies (as recorded by sketched cross-section relationships and map-view trajectory).

These scenarios are broadly consistent with previously published interpretations but differ in their detail to show the flexibility of SBIM to prototype 3D interpretations and rapidly create 3D geological models from sparse datasets.

Initial interpretation scenario SM1

The initial interpretation (SM1; Fig. 2) consists of three parasequences separated by sketched flooding surfaces. Each parasequence is subdivided into facies belts by sketched facies-bounding surfaces, consistent with the outcrop sections and borehole data (Fig. 1). The middle parasequence contains a flood tidal delta, consistent with observations at the Sowbelly Gulch section (Fig. 1), and this is captured by sketching surfaces that represent the top and base of the flood tidal delta deposits, which are considered to form a continuous facies belt in this initial interpretation. The upper parasequence contains a tide-influenced channelized sandbody, consistent with observations at the Gentile Wash section (Fig. 1), and this is captured by sketching surfaces that represent the top and base of the sandbody. In this initial interpretation, the channelized sandbody has been interpreted to comprise a single, tide-influenced deltaic distributary channel that was active during regression of a nearly linear shoreface and associated strandplain. The cross-section view in Figure 2 (top) shows the vertical (2D) sketched correlation. However, in the subsurface, the 3D geometry of the parasequences, facies belts and tidal deposits could not be deduced from the available data; even in this outcrop section, the relative lack of 3D control means there is significant uncertainty.

Here, we exploit the convenience and speed of SBIM to extend the 2D correlation into 3D, based on conceptual understanding of the depositional system. We extrude each surface sketched in cross-section along a trajectory sketched in map-view (Fig. 2, bottom right). We assume a close-to-linear shoreline in this initial interpretation. The trajectories sketched for the parasequence-bounding flooding surfaces represent the interpreted shoreline palaeogeography of each parasequence; these trajectories are re-used for the facies-bounding surfaces, consistent with the interpretation that the facies belts, including the flood-tidal delta deposits, are continuous in the shoreline-parallel direction. Also shown is the sketched trajectory of the tide-influenced channelized sandbody (red), which corresponds to the channel thalweg and is interpreted to be oriented approximately perpendicular to the shoreline. The 3D perspective view in Figure 2 (bottom left) shows the resulting 3D model, created entirely using SBIM. SM1 is the only scenario that needed to be sketched from ‘blank screen’; the other scenarios all modify this initial 3D model, as outlined below.

The sketching workflow for this example simply involves sketching the correlation on the cross-section view (top) and adding a map-view trajectory. All boundaries (flooding surfaces, facies boundaries) are sketched in a single east–west cross-section, except for the boundary of the tide-influenced channelised sandbody, which is sketched in a north–south cross-section.

We initially choose to sketch in hierarchical order, from the largest to the smallest scale. Each 2D sketch is extruded along the associated sketched trajectory to create a 3D surface before the next surface is sketched; we focus here on how the surfaces are sketched in cross-section. The parasequence boundaries (flooding surfaces) are sketched first from base to top, using operator RA or RAI to remove overlapping portions of earlier sketched surfaces (surfaces 1–3 in Figure 3a; screen recording available in the supplementary materials). Next, we add the surfaces denoting the boundaries between facies belts, using the operator PBW to ensure the sketched surfaces are truncated against the parasequence-bounding surfaces.

In the lower parasequence, the PBW operator is activated between its bounding flooding surfaces, and facies boundaries are sketched from bottom to top using operator RA or RAI (surfaces 5–9 in Fig. 3b). In the middle parasequence, the lowest three facies boundaries are also sketched from bottom to top combining PBW with RA or RAI (surfaces 11–13 in Fig. 3c). Next, the flood tidal delta is sketched, again using operator RA or RAI (surface 14 in Fig. 3c). A new sketching region within which to apply PBW is then selected above the flood tidal delta and below the upper parasequence boundary (Fig. 3c). The final two facies boundaries are then sketched within this region using RA/RAI (surfaces 16–17 in Fig. 3c), with use of PBW ensuring they terminate against the top of the flood tidal delta.

In the upper parasequence, the lower three facies boundaries are sketched from bottom to top using RA/RAI, followed by the upper facies boundary (surface 22 in Fig. 3d). This latter surface also represents the top of the channelized sandbody; the operator PB is applied to this surface and the channelized sandbody is sketched using RA/RAI: the resulting erosional channel base truncates earlier sketched surfaces but is itself truncated against the channel top surface (surface 27 in Fig. 3d). A new sketching region within which to apply PBW is then chosen, below the base of the channelized sandbody and above the underlying facies boundary. The final two facies boundaries are then sketched in this region using RA/RAI.

The sketching approach described above is just one of many that can be followed to create the initial interpretation model and we show three alternative approaches that yield the same final 3D model in Figure 3. It is this flexibility in sketching approaches that makes SBIM combined with geological operators so powerful. In general, when sketching surfaces in order from bottom to top, the operators RA and RAI are commonly used (e.g. Fig. 3e–h). The same outcome can be achieved when sketching from top to bottom using operators RB and RBI (e.g. Fig. 3i–l). PBW is a useful operator when sketching in hierarchical order, as surfaces at lower levels of the hierarchy can be constrained to lie within surfaces at higher levels. However, there is no need to sketch in hierarchical order. Some users might prefer to sketch in stratigraphic order (e.g. Fig. 3e–h). Some may prefer to sketch in any order, with no regard to hierarchy or stratigraphy (e.g. Fig. 3m–p). Regardless of the approach used, the 3D model can be created in less than 10 minutes.

Perspective views of different parts of the initial interpretation model (Fig. 4) show the spatial relationships and connectivity of facies belts and geobodies in 3D. Every parasequence consists of a repeated pattern (Fig. 4a–f), with facies belts stacked vertically and laterally from offshore transition (grey) in the deepest and most distal position (east) overlain by lower shoreface (distal: brown; proximal: orange), upper shoreface and foreshore (yellow/pale yellow) to coastal plain facies (green) in the shallowest and most proximal position (west). The linear belt of flood tidal delta deposits in the middle parasequence (Fig. 4d and e) lies at the up-dip termination of the foreshore and upper shoreface facies belts.

Offshore transition (grey) and coastal plain (green) facies are predominantly composed of very low-permeability mudstones, and may form significant aquitards or barriers where they are laterally extensive and continuous. Conversely, upper shoreface (yellow), foreshore (pale yellow), channel-fill and flood-tidal delta (red) facies are high-permeability sandstones, with best aquifer or reservoir potential. As a result, each parasequence is composed of an aquifer/reservoir in its upper part and an aquitard/barrier in its lower part. Parasequences pinch out to the west (palaeo-landward) into a coastal plain aquitard. In this initial interpretation, the lateral continuity of offshore transition aquitards limits the vertical connectivity of aquifer sandstones in each parasequence, with hydraulic communication only possible through heterolithic distal lower shoreface facies (brown) and low-permeability proximal lower shoreface facies (orange). The large dip extent of shoreface facies belts and strike continuity of flood-tidal delta facies in this initial model result in large aquifer volume and high lateral aquifer continuity (Fig. 4d and e). The channelized sandbody in the upper parasequence increases aquifer/reservoir volume and local lateral continuity, and provides the potential to locally erode through the lower part of the parasequence, thus providing vertical connectivity with aquifer sandstones in the underlying parasequence (Fig. 4b).

Alternative interpretation scenarios

Unlike this outcrop example, subsurface datasets typically contain insufficient information to constrain the dip extent of facies belts between boreholes, and this is an important uncertainty to consider. In scenario SM2, the down-dip extents of shoreface facies belts are modified, which results in a different parasequence stacking pattern (Fig. 5a and b). The initial model exhibited a progradational-to-aggradational parasequence stacking pattern (Kamola and Huntoon 1995), whereas this second scenario reflects a uniformly aggradational stacking pattern (Fig. 5b).

Depending on the sketching order used to create the initial model SM1, the sketched parasequence bounding surfaces can be re-used, simply by undoing sketches of facies-bounding surfaces (steps 4–20 in procedure #1 of Fig. 3a–d) before sketching the new surfaces. Alternatively, a combination of PA and PB can be used to select, respectively, the base and top bounding surface of a parasequence, and RA or RB used to remove all facies-bounding surfaces within the selected parasequence before sketching the new surfaces. The sketching process for this second scenario is then identical to that used for the initial interpretation model, including re-use of the same plan-view trajectories (screen recording available in the supplementary materials), except that we vary the sketches of the facies-belt boundaries so that their down-dip terminations occur in different locations compared to SM1.

The modelled dip extent of facies belts has a quantifiable impact on relative facies proportions (Fig. 5b), representing a reduced extent and volume of aquifer sandstones in each parasequence. Offshore transition and coastal plain aquitards retain their lateral extent and continuity, such that aquifers remain poorly connected vertically.

Scenario SM3 explores variations in the map-view trajectory of shoreface facies belts, which reflects shoreline orientation and rugosity (Fig. 5c and d). The cross-sectional correlation between the measured outcrop sections and boreholes is maintained from the initial model SM1, but the orientation and rugosity of shoreface facies belts in each parasequence are sketched using a different map-view trajectory. The same trajectory is used for the upper bounding surface of each parasequence and the facies-bounding surfaces within that parasequence. Here we modify the map-view trajectories for all three parasequences. However, an individual parasequence could be modified while the others remain unchanged. The sketching process for SM3 is identical to the initial interpretation, differing only in the map-view trajectories used (screen recording available in the supplementary materials).

The facies proportions in model SM3 change only marginally relative to the initial model SM1 (Fig. 5d), but the distribution of high-quality aquifer sandstones is modified in each parasequence. The area of foreshore and upper shoreface sandstones that are overlain by coastal plain mudstones at the top of the middle parasequence is increased (Fig. 5d), such that the vertical connectivity between these sandstones and overlying low-quality proximal lower shoreface sandstones in the upper parasequence is decreased.

In scenario SM4, a different geometry is considered for the flood tidal delta deposits in the middle parasequence (Fig. 5e and f). In the previous scenarios, these deposits were interpreted as a linear facies belt aligned parallel to the shoreface (Fig. 2). In scenario SM4, the flood tidal delta deposits are instead interpreted as landward protruding lobes that extend from the shoreface (Hampson and Howell 2005). The sketching procedure is identical to SM1, except a different approach is used to sketch the lobate geometry of the flood tidal delta deposits. The top of each flood tidal delta is constructed by sketching contours on horizontal planes, one at the top of the delta and one at the base (Fig. 6a and c). An additional sketching plane can be added in the middle for extra control on the final geometry (Fig. 6b).

The sketching process is shown in the supplementary screen recordings, starting from the initial interpretation, removing the middle and upper parasequence and finally sketching the new geometry for the flood tidal delta deposits. The SBIM approach is, however, sufficiently flexible to allow an alternative approach, in which these surfaces are sketched in multiple vertical cross-sections, either in dip section or strike section (Fig. 6d–f). For example, on strike-parallel cross-sections, one sketch to delineate the highest point (Fig. 6e), and two cross-sections at the distal and proximal termination (Fig. 6d and f) are sufficient to generate the top surface of the delta.

Facies proportions in model SM4 (Fig. 5f) change slightly from those of the initial model, with an increased proportion of coastal plain facies and a reduced up-dip extent and volume of high-quality aquifer sandstones in the middle parasequence. This results in a reduced volume and less continuous distribution of high-quality flood-tidal delta sandstones that are no longer connected.

Finally, scenario SM5 tests a different interpretation of the tide-influenced channelized sandbody that cuts into the top of the upper parasequence (Fig. 5g and h). In scenario SM5, the sandbody is instead interpreted to form part of a tide-influenced deltaic distributary channel network (Fig. 5h); the geometry of the underlying shoreline and shoreface facies belts is also modified to represent the arcuate geometry of a wave-dominated delta rather than a linear strandplain. To sketch scenario SM5 from the initial interpretation, the upper parasequence is removed (undo last eight sketched surfaces) and a new version sketched. In this example, the same order of sketching is used as in SM1. The curved geometry of the shoreface facies belts in the upper parasequence is constructed simply by using a different map-view trajectory. The distributary channel network is sketched as contours on two horizontal planes, one for the channel bases and one for the channel tops. A screen recording showing the sketching workflow for this model is available in the supplementary materials.

Facies proportions in model SM5 vary slightly from the initial interpretation SM1 (Fig. 5h), with an increased proportion of tide-influenced channel-fill deposits. The area of proximal lower shoreface sandstones in the upper parasequence that overlie foreshore and upper shoreface sandstones at the top of the middle parasequence is reduced, such that the vertical connectivity between these sandstones is decreased. Increased volume of channelized sandbodies in the upper parasequence enhances their potential to locally erode through the lower part of this parasequence, and thus provide vertical connectivity with the aquifer sandstones in the middle parasequence.

Case study II: normal fault array – outcrop cross-section

Context and geological setting

It is well known that faults can significantly affect fluid flow in aquifers and reservoirs, and impact seal integrity for CO2 storage (e.g. Aydin 2000; Bense and Person 2006; Dockrill and Shipton 2010; Frery et al. 2015). In this second case study, we demonstrate that faults can also be captured in a 3D model using SBIM (Fig. 7). The case study models faulted, interbedded sandstones, mudstones and limestones in an exposed part of the footwall damage zone of the Moab Fault (Foxford et al. 1996; Doelling et al. 2002). The exposure contains an array of normal faults in a nearly vertical road cut next to Highway 191, opposite the main entrance to Arches National Park, Utah, USA. The beds in this outcrop belong to the upper part of the Honaker Trail Formation (Late Pennsylvanian), deposited in shallow-marine shelf and nearshore environments (Doelling et al. 2002).

The aim of the modelling is to test different fault configurations that are prominent at the scale of the outcrop (110 × 40 m), ignoring smaller scale fractures and deformation bands. Faults in this outcrop (Fig. 7) show a synthetic or antithetic orientation relative to the main Moab Fault (Torabi et al. 2019).

We use SBIM to delineate bed boundaries (from bottom to top: purple/red/yellow/green) as well as faults (west-dipping: dark blue; east-dipping: light blue) by sketching over the outcrop photograph (Fig. 7 top). Although the features in this outcrop are easy to recognize, there is uncertainty in how the faults may be interpreted and modelled away from the outcrop face, and we use SBIM to rapidly create 3D models that explore and test different scenarios. In some of these scenarios, we deviate from the published interpretation (Torabi et al. 2019), in order to demonstrate the flexibility of the SBIM approach, sketching a conjugate fault set in which antithetic faults exhibit a different strike orientation relative to the main Moab Fault. The scenarios tested are as follows.

  • Scenario NF1: initial interpretation; antithetic faults with strike parallel to Moab fault

  • Scenario NF2: antithetic faults forming intersecting conjugate fault sets

  • Scenario NF3: antithetic faults forming intersecting conjugate fault sets with two additional fault blocks.

Initial interpretation scenario NF1

In the initial interpretation model, strike-parallel synthetic and antithetic faults are sketched (respectively dark and light blue in Fig. 7), as reported by Torabi et al. (2019). The plane of the outcrop face is assumed to be oriented perpendicular to the strike direction of all faults.

This initial model is constructed by first sketching all the faults, and then sketching the bed boundaries within each fault block, in the outcrop cross-section (Fig. 8a–d). Sketching of the faults can begin with any fault (Fig. 8a and b). For each fault, the sketching operator to apply depends on whether it is in the footwall or hanging wall of previously sketched faults. When sketching a new fault in the footwall of a previously sketched fault, RB or RBI is most appropriate. When sketching a new fault in the hanging wall of a previously sketched fault, RA or RAI is most appropriate. Alternatively, PBW can be used to select two faults and sketch within all fault blocks that they bound. Sketching the bed boundaries within a fault block also relies on PBW, to constrain the boundaries within the selected block (Fig. 8c and d). Within a fault block, bed boundaries can be added in any order, similar to Case study I (Fig. 3). A step-by-step screen recording of this sketching approach is available in the supplementary materials.

An alternative approach is to mix sketching of stratigraphy and faults (Fig. 8e–h). Some bed boundaries are sketched first (Fig. 8e), and are subsequently truncated by a sketched fault (Fig. 8f). If the sketched bed boundaries lie in the footwall of the new fault surface, the RA or RAI operator is used to truncate the bed boundaries correctly, and RB or RBI if the previously sketched boundaries lie in the hanging wall. This will generate two fault blocks, one with bed boundaries already sketched, and another that is empty. Bed boundaries (or additional faults) are then sketched within the empty fault blocks (Fig. 8g and h). The combination of sketch operators allows faults and fault blocks to be added after stratigraphy has been sketched.

Once a fault or bed boundary has been sketched in the outcrop cross-section (Fig. 8), it is then extruded perpendicular to the sketching plane, from front to back of the model. Consequently, all cross-sections parallel to the outcrop face are identical (Fig. 9a–c). The perspective views in Figure 9 show how the continuity of some beds is affected by faulting. Assuming that transmissibility of the faults is not impacted by clay smear or cementation (i.e. there is no membrane sealing), most beds are still connected across most faults. However, in this initial interpretation model, no connectivity is preserved with the westernmost fault block (highlighted by the blue rectangle in Fig. 9b). Thinner beds (e.g. yellow bed in Fig. 9a) are generally disconnected. If these thin beds were conduits to flow, the faults would create flow barriers by juxtaposition sealing. Conversely, if the thin beds were barriers to flow, the faults would disconnect them laterally such that they form baffles around which fluids can pass.

Alternative interpretation scenarios

In scenario NF2 (Fig. 9d–f), different strike orientations are sketched for the faults interpreted as synthetic (Fig. 7; dark blue) and antithetic (Fig. 7; light blue) to represent intersecting conjugate fault sets. In the outcrop cross-section, the interpretation does not change and the sketching workflows in Figure 8 are applied, but the synthetic faults are extruded along a NNW–SSE trajectory, whereas the antithetic faults are extruded along a NNE–SSW trajectory. The 3D geometry of the model shows the more complex relationships between fault blocks (Fig. 9d–f). The differences in fault orientation are most obvious in the lowest beds of the model (Fig. 9e). An important difference between this scenario and NF1 is visible on the right-hand (NW) corner of the model (Fig. 9e; highlighted by the blue dashed square), where the geometry of the fault block is changed and yellow and light green beds are exposed, in contrast to the initial model (Fig. 9b; highlighted by the blue dashed square). Because the model is no longer a perpendicular extrusion of the sketched cross-section, not every volume is intersected and visible in a single cross-section. Moving to different cross-sections towards the front or back of the model exposes the lack of some beds, which must be sketched on different cross-sections. A step-by-step sketching video is available in the supplementary materials.

The change in strike of the faults to a conjugate set rather than parallel faults has an impact on bed connectivity. In comparison to scenario NF1, most beds remain connected from west to east, but now also including connectivity with the westernmost fault block (Fig. 9b and e; highlighted by the blue dashed squares); here, the orange bed is only just connected with the westernmost fault block. East–west fluid flow through the orange bed would be possible in this scenario. There is no change in connectivity of the thin beds, as they are not juxtaposed across any fault.

In scenario NF3 (Fig. 9g–i) additional fault blocks are modelled by sketching additional antithetic faults in model NF2 (highlighted by the blue arrows in Fig. 9d and g and blue circles in Fig. 9f and i). Because the newly sketched faults and fault blocks do not intersect the outcrop cross-section, this is an equally possible scenario as NF2. To sketch the additional fault blocks, PBW is used to constrain new sketched fault surfaces within an existing fault block(s). RA or RB is used depending whether the existing bed boundaries or faults are respectively in the footwall or hanging wall of the new sketched fault. The existing fault blocks are truncated by this new fault, and new bed boundaries are sketched within the footwall of the new fault (steps 37–42 in workflow #2 of Fig. 8h). Structure and stratigraphy do not have to be sketched in a prescribed order but can be sketched in different combinations. Adding a new fault and sketching bed boundaries inside the new fault block was completed in less than four minutes. A step-by-step sketching video is available in the supplementary materials.

The addition of two new antithetic faults has no impact on the connectivity of thicker beds (e.g. red and orange beds in Fig. 9h and i are still connected). However, the additional faults have caused further disconnection of thin beds (e.g. yellow bed highlighted by the blue arrow in Fig. 9d and g).

Discussion

Our new SBIM modelling approach allows, for the first time, geoscientists to create 3D geological models using an intuitive approach that is similar to sketching in 2D on paper. Such sketches are familiar to geoscientists as they are often used to conceptualize and communicate geology in maps, cross-sections and block diagrams. Our geological operators allow surfaces to be sketched in any order so geoscientists can interpret as they sketch, prototyping geological concepts and testing alternative conceptual models and scenarios. 3D models can be created very rapidly compared to most currently available expert software and with minimal training. Sketching a model from a blank screen typically takes minutes to tens of minutes, and models can be sketched in the field using tablets with touch-sensitive screens, as well as in the office.

The resulting models are quantitative and, although a detailed description of such calculations lies outside the scope of this paper, they can be used to determine key aquifer and reservoir properties. Volumetric calculations, such as total or connected volumes and fluid storage capacity of a given facies, can be obtained directly from the sketched model. However, we have also integrated our implementation of SBIM for 3D reservoir modelling with computationally cheap flow diagnostics. Such flow diagnostics allow key flow properties and behaviour to be assessed using a single pressure solution (Zhang et al. 2017, 2020). Together, SBIM and flow diagnostics allow rapid, quantitative assessment of the impact of different geological concepts and scenarios on resource volumes and flow properties. Figure 10 shows an example of our flow diagnostics module in action on sketched model SM5. Here, the pressure field obtained for a given combination of injection and offtake boreholes allows the fluid ‘time of flight’ to be calculated, highlighting key flow paths through connected aquifer facies. Petrophysical properties are assigned to facies or other geologically defined rock units using a simple interface, and boreholes can be interactively added and moved around the model with flow properties updated essentially instantaneously.

Comparison of sketch-based modelling, as implemented here, to other modelling tools reveals a number of advantages and disadvantages. The advantages of sketch-based modelling include ease of use and flexibility in the features that can be modelled. Sketching a geological model is intuitive: it does not require users to learn complicated workflows and their underlying methods. The operators introduced here allow users to sketch surfaces in any order, rather than being forced to follow stratigraphic ordering, a length-scale based hierarchy, or order imposed by a pre-defined workflow. Intricate geometries, which are hard to replicate using existing modelling approaches, can be sketched in 3D. The disadvantages of sketch-based modelling are inherent in its design as a propotyping tool. Despite its capability to generate models quickly, sketch-based modelling is not designed to generate hundreds or thousands of realizations of the same base case as in conventional stochastic modelling tools and workflows, but rather to explore scenarios that are conceptually different. We emphasize that our aim is not to replace conventional workflows, but rather to complement them. Finally, sketching is easiest on touchscreen-enabled devices, such as laptops and tablets, either using a stylus or fingers. Sketching with a computer mouse is less natural and requires the user to find a balance between precision and speed.

Potential applications of our implementation of SBIM for 3D geological modelling are numerous and include characterization of subsurface groundwater and hydrocarbon resources and potential sites for CO2 storage. There are also broad potential applications in geoscience education and training. 3D models can easily be sketched in the classroom or the field, so that qualitative observations and quantitative calculations can be assessed directly and communicated with peers. The quick turnaround time from conceptual sketch to model output encourages geological thinking in 3D, and integration with flow diagnostics demonstrates the impact of geological interpretation on predicting subsurface flow behaviour.

The current set of geological operators are sufficiently robust for any stratigraphic model and a range of simple structural scenarios. Ongoing development of RRM is focused on extending the sketching of structure to include more complex scenarios including faults that tip-out within the model domain, and on further simplifying the user experience, for example by allowing users to edit existing models and sketch template surfaces that can be re-used to avoid repetitive sketching. We are also expanding the geological operators to allow a wider range of structural settings and diagenetic alteration to be modelled.

Conclusions

In this paper we present the first application of SBIM for rapid and intuitive creation of realistic 3D geological models. We integrate SBIM with geological operators to allow a flexible approach to sketching. The result is a fast, geologically robust prototyping tool that leverages traditional techniques to conceptualize geological interpretations such as maps, cross-sections and block diagrams, without requiring specialist modelling expertise. The resulting models can be assessed visually and their volumetric properties and dynamic flow behaviour evaluated quantitatively. In addition to a wide range of applications for resource estimation, there are educational benefits of our modelling approach in developing 3D geological interpretation and visualization skills.

Acknowledgments

The authors want to thank the members of Phase 1 of the Rapid Reservoir Modelling Consortium (Equinor, ExxonMobil, Petrobras, Shell and IBM Research Brazil/IBM Centre for Advanced Studies, Alberta, Canada) and Phase 2 of the Rapid Reservoir Modelling Consortium (Equinor, ExxonMobil, Petrobras, Petronas and Shell) for granting permission to publish this paper. Geiger acknowledges partial funding for his Chair from Energi Simulation. The authors thank Subject Editor Sarah Boulton, reviewer Tore Grane Klausen and an anonymous reviewer for their insightful and constructive comments on the manuscript.

Author contributions

CJ: conceptualization (equal), data curation (lead), formal analysis (equal), funding acquisition (equal), investigation (equal), methodology (equal), software (equal), validation (equal), visualization (lead), writing – original draft (equal), writing – review & editing (lead); MEHP: conceptualization (equal), formal analysis (equal), funding acquisition (equal), investigation (equal), methodology (equal), software (equal), validation (equal), visualization (equal), writing – original draft (lead), writing – review & editing (equal); GJH: conceptualization (equal), formal analysis (equal), funding acquisition (equal), investigation (equal), methodology (equal), software (equal), validation (equal), visualization (equal), writing – original draft (equal), writing – review & editing (equal); MDJ: formal analysis (equal), funding acquisition (equal), investigation (equal), methodology (equal), software (equal), validation (equal), visualization (equal), writing – original draft (equal), writing – review & editing (equal); DP: formal analysis (equal), funding acquisition (equal), investigation (equal), methodology (equal), software (equal), validation (equal), visualization (equal), writing – original draft (supporting), writing – review & editing (equal); SG: conceptualization (equal), formal analysis (equal), funding acquisition (equal), investigation (equal), methodology (equal), software (equal), validation (equal), visualization (equal), writing – original draft (supporting), writing – review & editing (equal); CCM: conceptualization (equal), formal analysis (equal), funding acquisition (equal), investigation (equal), methodology (equal), software (equal), validation (equal), visualization (equal); JDMS: conceptualization (equal), formal analysis (equal), funding acquisition (equal), investigation (equal), methodology (equal), software (equal), validation (equal), visualization (equal), writing – original draft (supporting), writing – review & editing (equal); SJ: conceptualization (equal), formal analysis (equal), funding acquisition (equal), investigation (equal), methodology (equal), software (lead), validation (equal), visualization (equal), writing – original draft (supporting), writing – review & editing (equal); FR: conceptualization (equal), formal analysis (equal), funding acquisition (equal), investigation (equal), methodology (equal), software (lead), validation (equal), visualization (equal), writing – original draft (supporting), writing – review & editing (equal); MCS: conceptualization (equal), formal analysis (equal), funding acquisition (lead), investigation (equal), methodology (equal), software (equal), validation (equal), visualization (equal), writing – original draft (supporting), writing – review & editing (equal)

Funding

Funding was received from Phase 1 of the Rapid Reservoir Modelling Consortium (Equinor, ExxonMobil, Petrobras, Shell and IBM Research Brazil/IBM Centre for Advanced Studies, Alberta, Canada) and Phase 2 of the Rapid Reservoir Modelling Consortium (Equinor, ExxonMobil, Petrobras, Petronas and Shell).

Data availability

All data generated or analysed during this study are included in this published article (and its supplementary information files). Screen recordings are available on https://doi.org/10.6084/m9.figshare.c.5084141. Rapid Reservoir Modelling prototype software and source code are available here: https://bitbucket.org/rapidreservoirmodelling/rrm

Scientific editing by Sarah J Boulton

This is an Open Access article distributed under the terms of the Creative Commons Attribution 4.0 License (http://creativecommons.org/licenses/by/4.0/)