Surface-Based Geological Reservoir Modelling Using Grid-Free NURBS Curves and Surfaces

Building geometrically realistic representations of geological heterogeneity in reservoir models is a challenging task that is limited by the inflexibility of pre-defined pillar or cornerpoint grids. The surface-based modelling workflow uses grid-free surfaces that allows efficient creation of geological models without the limitations of pre-defined grids. Surface-based reservoir modelling uses a boundary representation approach in which all heterogeneity of interest (structural, stratigraphic, sedimentological, diagenetic) is modelled by its bounding surfaces, independent of any grid. Volumes bounded by these surfaces are internally homogeneous and, thus, no additional facies or petrophysical modelling is performed within these geological domains and no grid or mesh discretisation is needed during modelling. Any heterogeneity to be modelled within such volumes is incorporated by adding surfaces. Surfaces and curves are modelled using a parametric non-uniform rational B-splines (NURBS) description. These surfaces are efficient to generate and manipulate, and allow fast creation of multiple realisations of geometrically realistic reservoir models. Multiple levels of surface hierarchy are introduced to allow modelling of all features of interest at the required level of detail; surfaces at one hierarchical level are constructed so as to truncate or conform to surfaces of a higher hierarchical level. This procedure requires the joining, terminating and stacking of surfaces to ensure that models contain “watertight” surface-bounded volumes. NURBS curves are used to represent well trajectories accurately, including multi-laterals or side-tracks. Once all surfaces and wells have been generated, they are combined into a reservoir model that takes into account geological relationships between surfaces and preserves realistic geometries.


Introduction
Surface-based reservoir modelling and simulation is an approach to modelling subsurface reservoirs (water, geothermal, petroleum) in which reservoir heterogeneity that impacts fluid flow is represented explicitly by their bounding surfaces (Caumon et al. 2009;Deveugle et al. 2011;Graham et al. 2015a, b;Hassanpour et al. 2013;Jackson et al. 2009Jackson et al. , 2013Jackson et al. , 2015Massart et al. 2016a, b;Pyrcz et al. 2005;Sech et al. 2009;Zhang et al. 2009). In this approach, all reservoir heterogeneity is represented by the use of surfaces, such that petrophysical properties within surface-bounded volumes are homogeneous (Jackson et al. 2013(Jackson et al. , 2015. The approach is capable of preserving detailed, complex input geometries over the entire modelling workflow, from geological representation to fluid flow simulation. The inputs for surface-based modelling are the bounding surfaces that define structural features, sediment bodies, facies and/or lithology distributions, and other types of geological heterogeneity. This is very different from conventional approaches in which geostatistical or object-based methods are used to distribute properties over pre-defined grid cells, and from hybrid surfaceand grid-based approaches in which some scales of heterogeneity are represented by surfaces that are then discretised on a cornerpoint grid on which further facies or petrophysical modelling is performed (Hassanpour et al. 2013;Pyrcz et al. 2005Pyrcz et al. , 2009Ruiu et al. 2016;Zhang et al. 2009).
Most conventional reservoir modelling and simulation workflows represent geological heterogeneity on (typically pillar) grids, which are defined early in the modelling process. However, it is challenging to represent many common geological features using pillar grids, including faulted domains and non-monotonic surfaces (e.g. recumbent folds, diapir flanks and the margins of intruded or injected bodies). Rock types with diverse petrophysical characteristics are "averaged" within grid cells of arbitrary size and shape, and the continuity and connectivity of low-permeability baffles and barriers or high-permeability zones is often lost. Stair-stepping effects are introduced by cornerpoint grids around any feature that is not aligned to the grid orientation. Such effects are common around faults, but also impact, for example, the continuity and connectivity of sinuous channels ( Fig. 1) (Jackson et al. 2013). Note the difference in connectivity of the pink and purple bodies in the two models. The effects become more pronounced as the grid element size is increased. It is also challenging to represent features across a wide range of spatial scales, even though important geological heterogeneity may occur across these scales, because the same number of pillars must be present across all layers. Consequently, there is little flexibility to adjust the areal grid resolution and large features are over-resolved and small features are under-resolved or omitted. For example, refining the areal definition of grid cells in order to capture a small-scale geological feature that is present in only one grid layer will result in refined grid resolution in all overlying and underlying grid layers, even though properties in these layers might be homogeneous. In this example, the grid refinement provides local geometrical accuracy (i.e. small grid cells) at the expense of computational efficiency (i.e. a large number of grid cells). Pyrcz et al. (2005) and Zhang et al. (2009) Jackson et al. 2013) recognise that these limitations also apply to surface-based modelling if surfaces are represented by gridded surfaces, such that adequate element resolution is needed to capture geological detail. The representation of wells in reservoir models is subject to the same issues. Discretising wells on a pre-defined grid introduces geometrical errors, since grid cells are several orders of magnitude larger in area than the wellbore, even if the grid is locally refined across all layers. Deviated wells are represented by stair-stepping geometries that may compromise their connectivity with geological features, and wells located close to inclined faults can be represented in the wrong structural domain. Some hybrid approaches model heterogeneity as a combination of surface-based volumes that are then infilled by pillar grids (Pyrcz et al. 2005;Ruiu et al. 2016;Xie et al. 2001;Zhang et al. 2009) or rastered onto locally refined grid cells (Hassanpour et al. 2013;Pyrcz et al. 2009), which are then used for further (geostatistical) facies or petrophysical modelling. These methods result in realistic heterogeneity distributions, but they require a high-resolution grid in order to capture the modelled heterogeneity, and, thus, an additional upscaling step from geological grid to simulation grid, which re-introduces the grid-related issues described above. In the modelling approach advocated in this paper, heterogeneity is represented by surface-bounded volumes that are internally homogeneous and, thus, no extra petrophysical modelling step is introduced.
A surface-based heterogeneity modelling approach is presented that is based on using non-uniform rational B-splines (NURBS) to generate surfaces that represent all the geological heterogeneity of interest without imposing it on a pre-defined grid or imposing a resolution. The resolution at which heterogeneity is represented depends on the scale of surfaces and surface-bounded volumes that are constructed, rather than the spatial resolution of grid cells. Our focus is on the generation and interaction of different surfaces across all length scales of their hierarchical arrangement, which are the essential elements of surface-based modelling without reference to a grid. These models are only discretised when required for numerical analysis; using geometryadaptive fully unstructured grids avoids potential grid-related and upscaling problems, and allows full advantage to be taken of adaptive, unstructured-grid simulators (Jackson et al. 2015;Salinas et al. 2017). This paper does not present new stochastic methods to generate surfaces or to condition models to input data, but a grid-free surface-based modelling approach that can leverage existing deterministic, stochastic and eventbased object-based modelling methods in these regards (Georgsen and Omre 1993;Hassanpour et al. 2013;Parquer et al. 2017;Pyrcz et al. 2009;Ruiu et al. 2016;Wang et al. 2018).
NURBS interpolation (Piegl and Tiller 1997) is an efficient way to represent curves and surfaces, and is commonly used in other modelling disciplines (e.g. computeraided design, civil and mechanical engineering, aerodynamics). The main advantage of NURBS surfaces is that complex non-monotonic shapes can be represented using only a few control points, which allows both geometrical accuracy and computational efficiency. NURBS surfaces may be used to represent many common geological heterogeneity types across a wide range of length scales. Wells may also be represented by NURBS curves, which allow geometrically realistic well trajectories that preserve connectivity and well locations.
Because reservoir modelling is conventionally started from a pre-defined grid, NURBS surfaces have not been used widely in this context. Smooth surface-based modelling based on discrete smooth interpolation methods was introduced by Wietzerbin and Mallet (1994), but the resulting surfaces differ from NURBS in that the control points lie on the surfaces. To date, NURBS surfaces have been used to delineate the large-scale subdivision of a reservoir Caumon et al. 2009;Florez et al. 2014;Jones et al. 2002;Zehner et al. 2015) and for fracture modelling (Belayneh et al. 2006;Corbett et al. 2012;Geiger and Matthäi 2012;Paluszny et al. 2007). The former generally use NURBS to construct a deterministic model. However surface-based stochastic approaches have been proposed to represent meso-scale architectural elements, such as channelised sediment bodies prior to cornerpoint gridding (Parquer et al. 2017;Ruiu et al. 2016), rasterising channel geometries onto a grid (Hassanpour et al. 2013) or to represent lobe element geometries that constrain heterogeneous cornerpoint grids (Pyrcz et al. 2005;Zhang et al. 2009).

Fig. 2
Non-uniform rational B-splines (NURBS) curve representation of five wells (dark blue to light blue) that extend into the subsurface from the same surface wellsite. Some control points are shared between different wells to ensure that different wells share part of a common trajectory The approach presented here uses NURBS surfaces in a similar way to this previous work, but across all length scales, both deterministic and stochastic, and to represent all heterogeneity of interest, without modelling on a grid. There is no scale or resolution limitation to the NURBS curves and surfaces used in the approach presented herein (e.g. none of the figures feature a scale bar).

Non-Uniform Rational B-Splines (NURBS)
NURBS curves and surfaces are parametric functions over, respectively, one or two parameters that result in a smooth curve or surface defined by a set of weighted control points, their interpolation degree and a parametric knot vector (Piegl and Tiller 1997).

Curves
NURBS curves are defined by Eq. (1), which yields an interpolation between a series of control points (Fig. 2). The shape of the resulting curve is a function of the location in x, y and z of the control points, the individual weights of the control points, the parametric influence ranges of the control points [knot vector, Eq. (2)] and the degree of the curve.
A NURBS curve C of degree p along parameter u is defined by in which P i are the n control points with coordinates (x i , y i , z i ), w i are the weights of the respective control points and N i, p (u) are the non-rational B-spline basis functions defined on a non-uniform knot vector where r n + p + 1 (Piegl and Tiller 1997). An important geometrical property of NURBS curves for modelling applications is local modification. If a control point P i is moved or its weight w i is changed, it only affects the portion of the curve in the parametric interval u ∈ [u i , u i+ p+1 ). This is useful when detail is added to a portion of the curve because it does not affect the complexity of the rest of the curve.

Surfaces
NURBS surfaces are defined by Eqs. (3) and (4), which yields an interpolation between a net of control points. The shape of the resulting surface is a function of the location in x, y and z of the control points, the individual weights of the control points, the parametric influence ranges of the control points [knot vectors, Eqs. (5) and (6)] and the degree of the surface. It should be noted that the degree of the surface can be different for u and v parameters.
A NURBS surface S of degree p along parameter u with n control points and degree q along parameter v with m control points, respectively, is the bivariate vector-valued piecewise rational function of the form (Piegl and Tiller 1997) in which P i, j are the n × m control point positions (x, y, z) i, j and rational basis functions in which w i, j are the weights of the respective control points and N i, p (u) and N j,q (v) are the non-rational B-spline basis functions defined on knot vectors where r n + p + 1 and s m + q + 1. Local modification for NURBS surfaces is analogous to that for curves. If a control point P i, j is moved or its weight w i, j is changed, it only affects the portion of the curve in the parametric intervals u ∈ [u i , u i+ p+1 ) and v ∈ [v j , v j+q+1 ).
The example of a NURBS surface illustrated in Fig. 3 shows how the control point position and degree of the surface affect the surface shape. The same net of control points (Fig. 3a) is used to generate a surface of degree 2 and 3, respectively (Fig. 3b, c). The order of a surface (degree + 1) determines the number of control points shaping the surface at any point on that surface. NURBS are a computationally efficient representation of curves and surfaces. For example, the surface in Fig. 3 only requires 24 control points with x, y, z coordinates and weights, and two knot vectors of, respectively, 7 and 9 elements (for a degree of 2) in the u and v directions.
The degree of a curve or surface determines the parametric smoothness of the respective curve or surface and the continuity at curve and surface connections. In practice, a degree of 1, 2 and 3 generates linear, quadratic and cubic curves and surfaces, respectively. Assigning a degree higher than 3 has few practical applications.

Geological Surface Construction
Geological domains are closed volumes that are bounded by at least one surface. However, most volumes will be bounded by several interacting surfaces. The construction of surface-based models, therefore, requires different NURBS surfaces to interact with each other, in order to form closed volumes bounded by surfaces. Surfaces that bound closed volumes have to satisfy two criteria: (I) they must have an appropriate geometry that represents the shape of (or part of) the geological heterogeneity to be modelled and (II) they must have "watertight" interaction with all neighbouring surfaces to form closed volumes.
Representing geological heterogeneity across a wide range of length scales requires a clear understanding of the spatial arrangement of heterogeneity, which is often hierarchical in character (Miall 1988(Miall , 1991Vakarelov and Ainsworth 2013;Van Wagoner et al. 1990). Such heterogeneity can be represented in surface-based models using a series of hierarchically arranged surfaces (Graham et al. 2015b;Sech et al. 2009). The spatial distribution of, and interaction relationships between, surfaces is dependent on their hierarchical level relative to other surfaces (Fig. 4). The appropriate geometry of such surfaces is defined by: (I) the shape of (part of) the geological heterogeneity to be modelled and (II) the hierarchical level(s) at which the modelled geological heterogeneity occurs, which determines its interaction with other geological heterogeneities. Surfaces are built directly to match both constraints or, instead, they can be built to represent the shape of the selected geological heterogeneity first and be then transformed ("warped") to match geometrical criteria imposed by surfaces at higher hierarchical levels. Depending on the data available, deterministic surfaces can Fig. 3 a An example of 24 control points, organised in six rows of four points (or four rows of six points), that determine the overall geometry of a surface. b Surface of degree 2 defined by control points in a, with equal weights for all control points. The right-hand side of the surface is non-monotonic. c A different surface based on the same set of control points, with the same weight, but with a degree of 3. This higher degree results in a smoother surface be created to represent one specific case, or surface geometry and distribution can be sampled from databases to produce stochastic surface-based models.
There are three types of interaction between surfaces that need to be considered in representing geologically realistic geometries: (I) joining surfaces at their edges to create closed volumes; (II) termination of hierarchically arranged surfaces; and (III) warping of surfaces and surface sets. These types of interaction are documented below in Sect. 2.5.  a Channelised geometry, constrained by its thickness, width, asymmetry and side angles. b Dipping clinoform template constructed by height, length and dip angle. c Channel top and twochannel lag and dipping accretion surfaces are constructed using variables relative to the channel geometry (a). d Facies boundaries inside a clinothem are constructed using variables relative to its bounding clinoforms (e.g. b)

Surfaces Extruded from Cross Sections
Many geological features that form heterogeneities can be described and quantified by their cross section, such as channelised sediment bodies, erosional scours and fault planes. The three-dimensional geometry of the feature is often determined by a trajectory along which this cross section is translated (extruded) and varied in shape. This information can be conveniently translated to a NURBS description by considering parameters u and v as the cross section and its extrusion trajectory, respectively. This approach is particularly useful when a well-defined geometrical model exists that defines the full three-dimensional geometry of the feature (e.g. Ruiu et al. 2016).
The definition of a NURBS cross section depends on the geometry to be modelled, and a new cross-section template is constructed for every type of geometry. Figure 4 shows geometrical templates for cross sections of channels and clinoforms. For most geometries, a degree of 3 is used for both the NURBS cross-section template and extrusion trajectory, in order to have the most control over the smoothness of geometry (Fig. 3). Weights are initiated at 1 for all control points and knot vectors are distributed uniformly. A few geometries require the weights and knot vectors of control points to be lower; for example, a degree of 2 and weights below 1 are used to generate ellipses (and ellipsoids).
The number and position of control points of the template depend on the information and constraints available for the geometry of the cross section (Fig. 4). For cross sections of degree 3, at least four control points are needed: one at each end of the cross section and two or more to shape the cross-section template. For each location that is determined exactly by the input geometrical information (e.g. channel thalweg), three collinear control points are added: a central control point exactly on the location and two control points on either side of and at the same distance from the central control point. This ensures that the template passes through the specified location defined by the central control point. The distance between these three control points depends on other constraints for the template. For every cross-section template, it should then be decided where the extrusion trajectory intersects the cross section (grey circles in Fig. 4a, b). The relationship between the cross-section template and extrusion trajectory may be different for each geometrical feature to be modelled, and should, thus, be chosen for each feature individually.
The construction of a surface now requires the positioning of the cross-section template along an extrusion trajectory (Fig. 5). The trajectory can be defined as a parametric function or as a list of coordinates through which it passes. The number of cross sections required to define the NURBS surface depends on the definition and complexity of the trajectory. For a linear trajectory, two cross sections are sufficient. However, most trajectories are likely to be more complex. For curved trajectories, an upper angle limit is used to sequentially position control points. Where the orientation of the trajectory deviates more than the chosen angle limit, a new cross section is inserted, until the end of the trajectory is reached. Cross-section templates are then translated and rotated such that they are perpendicular to the extrusion trajectory (Figs. 3 and 5). For trajectories that are defined by a list of coordinates, a curve construction procedure analogous to that used in modelling wells (see below) can be used. The number of cross sections and trajectory knot vector is then equal to the number of control points and knot vector of the NURBS curve. Every cross section of a surface uses the same template; however, the positioning of the control points within each cross section is also controlled by the input parameters of the cross-section template, and, thus, the final shape of the cross section can change along a trajectory (Fig. 6).

Joining Surfaces at Edges
If two surfaces have the exact same edge geometry, then they form a "watertight" connection (Fig. 7). The construction of a surface of which an edge corresponds to the edge of another surface is achieved by positioning the control points along the shared edge (CP1 and CP2 in Fig. 4c). Linking surfaces that have the same degree in the direction (u or v) of their shared edge is as straightforward as choosing the same coordinates for the control points of both surfaces along the shared edge. Surface edges of different degrees can be joined by a degree elevation procedure (Huang et al. 2005;Piegl and Tiller 1994;Prautzsch 1984), where the edge of the surface of higher degree is a degree-elevated copy of the edge of the lower-degree surface. If the tangents  (Fig. 4a). The geometrical parameters describing the channel shape (depth, width, asymmetry, side angles) are varied along the extrusion trajectory (green) between surfaces should be continuous, more control points may need to be added along the surface joints. For non-abrupt changes in surface orientation across a shared edge, the adjacent control points on either side of the edge need to be positioned collinear with the control points on the edge.

Termination of Hierarchically Arranged Surfaces
The construction of a surface that terminates on a different surface can be defined for a general case and also for the specific case where the terminated and terminating surfaces have a shared trajectory. The general case of intersecting/terminating surfaces is complex; however, determining the edge between surfaces becomes much simpler when the position of the terminated edge can be defined relative to the terminating surface hierarchically (Fig. 4c, d). Transformation of control points can then be used to define the control points for the edge. Existing geological volumes, which represent large-scale units at the higher level(s) of a hierarchy describing geological heterogeneity, can be subdivided into successively smaller sub-volumes, which represent smaller-scale units at lower levels of the hierarchy, by the addition of internal surfaces (Fig. 8). To ensure that these new sub-volumes are closed, the internal surfaces must connect to the existing (outer) bounding surfaces. This implies that surface-based models should be designed from the top down, starting at a higher hierarchical level and then moving to progressively lower hierarchical levels. However, the application of surface terminations is not limited to representing different hierarchical levels of geological heterogeneity, but is also essential to attach surfaces of the same hierarchical level. In general, termination operations are used to sculpt surfaces and volumes that have a clear genetic link to their parent surfaces in terms of geometry, orientation or location. Surfaces of a lower hierarchical level are Fig. 7 Two surfaces (green channel base; yellow channel top) are joined together to form a "watertight" connection. The surfaces are of different degrees in their cross section (degrees 3 and 1, respectively) but share the same trajectory, and, thus, the connection can be achieved by positioning the control points of the top surface (blue) coincident with the control points of the basal surface edge (red) constructed in relative coordinates to surface(s) of the higher hierarchical level with which they interact.
The termination of surfaces that share a cross section or extrusion trajectory with another surface (e.g. bar macroform boundaries within channelised bodies in Fig. 4c) is an operation that is similar to generating "watertight" shared edges. The edge of one surface needs to correspond to or intersect another surface. These terminated surfaces are generated with the same degree along the edge that is touching the terminating surface. The position of the control points on the terminated edge are positioned to correspond to the terminating surface. Ideally, the cross sections of terminated surfaces are positioned within the cross sections of terminating surfaces (Fig. 4c, d).
For the termination of two surfaces that do not share a cross section or extrusion trajectory (e.g. tributary channels, levees and lateral splays adjacent to trunk channels), the terminated surface must extend across the terminating surface. A "watertight" connection is then achieved by a Boolean operation that cuts away the terminated surface where it extends beyond the terminating surface. This operation ensures that there are no gaps between surfaces.

Warping of Surfaces
Warping surfaces to generate geometrical relationships, for example to superimpose surfaces onto pre-existing depositional topography or to mimic their deformation by post-depositional structure, is achieved by a two-step approach (Fig. 9). Firstly, the surfaces to be warped are constructed in their appropriate areal position and orientation, but using relative vertical coordinates with respect to a horizontal reference plane.

Fig. 9
Warping of surfaces to conform to the geometry of a reference surface. A downstream-widening and downstream-shallowing scour (yellow) running down a clinoform surface (grey). Similarly, a concretionary body along part of the clinoform is generated by constructing a bounding surface (blue) with reference to a horizontal surface that represents the depth of cementation below the clinoform Subsequently, the control points describing the surfaces to be warped are repositioned vertically according to their parent surface(s). The vertical position of each control point on the surfaces to be warped is adapted from the horizontal reference plane to its corresponding vertical position relative to the parent surface.
Surfaces and sets of surfaces are also commonly distorted by post-depositional faulting and folding. In principle, such distortion may be represented by the warping of surfaces, as outlined above, in which the control points describing the surfaces to be warped are translated both vertically and horizontally to define faults and folds (Fig. 10). Such translation would represent the forward modelling equivalent of (inverse) structural restoration, for which there are several well-defined procedures that are appropriate for particular structural settings (Caumon et al. 2009;Laurent et al. 2013Laurent et al. , 2016.

Well Construction and Interaction
The number, location, trajectory and completion design of wells that are required to inject fluids into or produce fluids from a subsurface reservoir are the subject of intense planning because they represent a key aspect of reservoir development and reservoir management. Defining these well-related parameters so that they allow optimal reservoir development is a common goal of many reservoir modelling studies. Well trajectories can be represented in a spatially accurate and computationally cheap manner as NURBS curves (Fig. 2).
Well trajectories are constructed as NURBS curves following the procedure outlined below (which may also be adapted to extrude a cross section along a complex Fig. 10 Three common structural configurations that can be accurately represented using NURBS surfaces, but not with conventional pillar grids: a faulted recumbent fold; b conjugate normal fault sets with unstructured mesh for flow simulations; c listric growth fault trajectory) (Fig. 2). Well trajectories are deterministic, and their input data (azimuth, dip and distance) define n points Q k in space through which the well should pass. The coordinates of n + 1 control points P i and parametric coordinates of the NURBS curve representing a well should, thus, be calculated to meet this requirement. This is achieved by first calculating the parametric coordinatesū k and knot vector {u 0 , . . . , u m } using the averaging technique described by Eq. (9.8) of Piegl and Tiller (1997) u 0 · · · u p 0, u m− p · · · u m 1, with either chord length (a 1) or centripetal (a 2) spacing methods [Eq. (8)] for parametric coordinates. Chord length spacing is the most widely used; however, the centripetal method provides better results for curves exhibiting sharper turns (Floater 2008;Lee 1989). Equations (7) and (8) are used to formulate a system of n + 1 linear equations to compute the knot vector and parametric coordinates. The coordinates of the control points can then be simply calculated by Eq. (9). This method results in a smooth curve along the well trajectory, Interaction between wells, associated with generating multi-lateral wells or sidetracks, is modelled by combining different well sections that share control points where they are joined (Fig. 2). The methodology above for one well is repeated for every lateral or side-track.

Model Construction
The NURBS surfaces and curves described above form the building blocks for surfacebased reservoir models. All surfaces and wells are constructed and then exported for assembly into a full reservoir model (Fig. 11). Because surfaces are not discretised, there are several approaches that can be used to create and mesh these reservoir models. Our approach exploits the parametric representation of geometry that is inherent in CAD systems. This exports surfaces in a standard, CAD-compatible format (.STEP) and employs the functionality of CAD engines to assemble the model, either automatically or manually (Melnikova et al. 2016). Alternatively, the NURBS surfaces can also be discretised and exported as triangulated surfaces (e.g. .TS, .OBJ, .PLY) and as gridded surfaces (e.g. .ZMAP) if they are monotonic, for integration into conventional reservoir modelling and visualisation packages.
Model assembly consists of three steps that are applied sequentially: (I) surface interactions (joining, terminating and stacking surfaces) follow rules based on the logic of superposition and cross-cutting relationships that arise from rock-forming processes (Jackson et al. 2015); (II) wells are merged with "watertight" volumes bounded by surfaces; and (III) "watertight" volumes are meshed (Melnikova et al. 2016). Each surface is assigned metadata that define its hierarchical level, relationship with other volumes and surfaces, such as removing older surfaces or being removed by younger surfaces, and the rock properties within the volumes it bounds. During model assembly, all geometrical modifications are performed through Boolean operations (intersection, subtraction, union) (Melnikova et al. 2016). Once all geometrical modifications have been completed, the model is meshed when required for numerical analysis (e.g. flow simulation or flow diagnostics, Fig. 10b), for example using a fully unstructured grid that ensures all geometries are preserved (Gomes et al. 2017;Jackson et al. 2015;Pellerin et al. 2017;Salinas et al. 2017;Zhang et al. 2017). Figure 6 shows three examples of NURBS surfaces that represent commonly observed geometries in sedimentary rocks: a channel, a downstream-widening and downstreamshallowing erosional scour, and six inclined surfaces (e.g. clinoforms).

Surface Extrusion
The cross-sectional geometry of a sinuous channel (Fig. 6a) can be reproduced by a NURBS curve of degree 3 using just five control points. This cross-section curve is then extruded along a sinuous plan-view trajectory. Based on the distribution of geometrical parameters that describe channel depth, width, cross-sectional curvature and sinuosity, control points are repositioned to construct similar cross sections along the trajectory. The extrusion trajectory of the channel and variability in cross-sectional channel geometry along the trajectory can also be stochastically generated by sampling from the geometrical parameter distributions. Sinuous channel trajectories can be generated by a number of methods (Parquer et al. 2017;Rongier et al. 2017;Ruiu et al. 2016;Viseur 2004).
More complicated channelised geometries are exhibited by scours that widen and become shallower downstream, as described and quantified by Jacquemyn et al. (2018) (Fig. 6b). In cross section, the scours exhibit pit-and-wing geometries of greater complexity than a simple symmetrical, concave-upward curve. Such geometrical complexity can be represented by increasing the number of control points along every cross section. The local modification property of NURBS surfaces ensures that additional detail is only added where needed and does not affect the geometry and complexity of the rest of the surface. The downstream widening and shallowing of the scours is represented by NURBS surfaces that have the same level of complexity as for a scour of constant cross section.
Clinoforms are inclined depositional surfaces (Rich 1951) whose geometries can be described and quantified by their cross section and plan-view trajectory (Graham et al. 2015a). In the NURBS-based approach described herein, clinoform surfaces can be modelled as cross sections (u-section) extruded along a trajectory (v-section) (Fig. 6c).

Shared Edges
The most evident way to form closed volumes is to join individual surfaces at their edges. For example, channel-top and channel-base surfaces should close a channel volume (Figs. 4 and 7). This type of interaction is achieved by generating the edges of the top surface in such a way that they conform to the edge of the base surface. Within every cross section, control points at the edge of the top and base surfaces are coincident (e.g. CP2 in Fig. 4c is coincident with CP7 in Fig. 4a). The coincidence in control point positions at the edge of the base and top surfaces is also applied in three dimensions, in order to ensure a "watertight" connection between the two surfaces (Fig. 7).

Surface Termination
The termination of hierarchically arranged surfaces is demonstrated here using two examples. The first example is a channelised sediment body deposited by a meandering river (Fig. 8a) (cf. Miall 1988). The base of the channelised sediment body (Figs. 6a and 8a) represents the highest hierarchical level of heterogeneity to be modelled. Five surfaces representing geological heterogeneity developed at lower hierarchical levels are terminated against the channel-base surface: the top of a channel lag and four lateral accretion surfaces that represent deposition on each inner bend of the channel during the deposition of point bars (Fig. 8a). The second example is a regressive shallow-marine deposit (parasequence sensu Van Wagoner et al. 1990) (Fig. 8b). The deposit contains a set of clinoform surfaces (Fig. 6c), which represents the highest hierarchical level of heterogeneity shown in Fig. 8b. Three surfaces that represent the boundaries of facies belts occur at a lower hierarchical level, and, thus, terminate at each clinoform surface (cf. Sech et al. 2009).

Surface Warping
Stratigraphic surfaces are constructed by deposition or erosion on top of pre-existing ("inherited") topography. The geometry of the resulting stratigraphic surfaces typically reflects the geometry of the underlying topographic surface, at least to some extent (e.g. Fig. 9, in which the geometry of the downstream-widening, downstream-shallowing scour shown in yellow conforms to the geometry of the grey-coloured clinoform into which it is eroded). Similarly, diagenetic surfaces and volumes may conform to the depositional and/or compactional geometries of underlying or overlying surfaces (e.g. Fig. 9, in which the geometry of the blue-coloured concretion follows the geometry of the overlying grey-coloured clinoform surface).

Case Study: Model of Net-Transgressive Deltaic Deposits
The full NURBS-surface-based model construction workflow is illustrated by building a model of net-transgressive deltaic deposits that is inspired by the SPE10 bench-marking study (Christie and Blunt 2001) and the Upper Brent Group reservoirs of the North Sea on which the SPE10 study was based (Fig. 11). The upper part of the Brent Group comprises a coastal plain succession that contains channelised fluvial sandbodies (Ness Formation) overlain by retrogradationally stacked, wavedominated shoreface parasequences (Tarbert Formation) (Livera 1989;Morris et al. 2003;Rønning and Steel 1987). The geometry and dimensions of channelised fluvial sandbodies in the coastal plain succession, and of clinoforms and facies belts in shoreface parasequences, are constrained by outcrop analogue measurements (Deveugle et al. 2011;Flood and Hampson 2015;Graham et al. 2015a;Villamizar et al. 2015), which have been modified to fit the dimensions of the SPE10 model and to demonstrate the capability of the NURBS-surface-based model construction workflow (Table 1). Multiple model realisations have been generated by sampling from the same input data (Fig. 11). The model shows that facies relationships and geometries are preserved correctly using the workflow, and that heterogeneities at scales across several orders of magnitude can be combined and represented explicitly within one model.

Discussion
In the following section, the advantages of surface-based models constructed using NURBS over conventional grid-based models are discussed. The three main advantages are: (I) stronger linkage with geological concepts and interpretations; (II) improved geometrical accuracy; and (III) greater computational efficiency.

Direct Conceptual Link Between Geological Interpretation and Reservoir Model
Surface-based modelling using NURBS surfaces preserves the direct conceptual link between geological interpretation of heterogeneity and a resulting reservoir model, as exemplified by the model shown in Fig. 6. Geoscientists typically conceptualise and communicate their interpretations in cross section(s) and map view. Consequently, the quantitative input data required to construct surface-based models are often available from studies of outcrop and modern analogues, geometrical databases, process-based models, and interpreted seismic lines and volumes. When the surfaces generated do not provide sufficient detail to represent the heterogeneity of interest, additional surfaces at lower hierarchical levels should be inserted. For example, this approach may be useful in representing gradational trends in petrophysical properties. Additional surfaces can be inserted to subdivide an existing volume into several smaller, sequential volumes that each contain incrementally higher or lower petrophysical properties; these sequential volumes together represent the trend in petrophysical properties at a level of accuracy that is determined by the number of inserted surfaces.
The models shown (Figs. 8,9,10 and 11) demonstrate the possibilities of surface-based modelling representing geological heterogeneity in a more geometrically realistic way than conventional modelling approaches. Further geometri-

Geometrical Accuracy
A NURBS surface representation preserves geometrical accuracy to whatever level of detail is needed. Explicitly representing heterogeneities that are several orders of magnitude smaller than the full model volume only increases the complexity of the model locally, where surfaces that capture the small-scale heterogeneities are added, and the rest of the model remains unaffected. In contrast, such heterogeneities cannot be mapped directly onto conventional pillar grids because they are too small or geometrically complex to be discretised on grid cells of a given size and shape, even in grid-based implementations of surface-based modelling (Pyrcz et al. 2005;Zhang et al. 2009). Incorporating such heterogeneities in pillar grids would require upscaling or complex grid manipulations. The resulting grid-based models typically compromise the underlying geological concepts if large grid cells are used, or are computationally expensive and inefficient if small grid cells are used. A higher degree of geometrical accuracy, which captures intended facies proportions and connectivity, is achieved by using NURBS surfaces without having to increase model resolution.
The advantage of using parametric surfaces such as NURBS to represent geological geometries in an accurate way can arguably be illustrated most clearly by considering common, but complex, structural configurations that are difficult or impossible to represent using conventional pillar grid approaches (Georgsen et al. 2012;Gringarten et al. 2008;Hoffman and Neave 2007;Holden et al. 2003;Qu et al. 2015). Figure 10 shows NURBS-surface-based models of three examples of such structural configurations. In these examples, structural deformation is represented by surfaces and surface sets that have not been warped (e.g. using the methods documented by Caumon et al. 2009;Laurent et al. 2013Laurent et al. , 2016, but instead have been constructed in their final structural configuration. Recumbent folds (Fig. 10a), which may also be faulted, are common in foreland fold-and-thrust belts, and adjacent to salt and shale diapirs. Beds in such settings may be locally overhanging (i.e. non-monotonic) and, thus, difficult to model with conventional pillar grids, which do not permit multi-value surfaces. Conjugate sets of intersecting normal faults (Fig. 10b) are also problematic for conventional pillar grids, even though all beds are monotonic (e.g. cases in Hoffman and Neave 2007). Growth faults, which have a listric geometry (Fig. 10c), are also problematic for pillar grids, especially in hanging wall regions where rotated growth strata are juxtaposed against the fault surface. Additionally, faults are not surfaces but have volumetric properties, which are often translated into transmissibility mapped on a modelled fault surface in conventional grid-based models. Surface-based modelling allows fault-core and fault-damage zones to be represented by volumes bounded by surfaces, similar to volumetric fault modelling (Fredman et al. 2008;Manzocchi et al. 2010), but without the representation of these zones on a cornerpoint grid.

Computational Efficiency
The examples described herein (Figs. 6,7 and 8) show how little information is needed to construct a NURBS surface of relatively complicated geometry, and NURBS surfaces are, thus, very efficient to construct, store and manipulate. This computational efficiency is a significant advantage over conventional reservoir modelling approaches when creating multiple model realisations that are sampled from the same input data distributions. Surfaces can be built quickly and cheaply because they do not require discretisation or gridding.
Using a NURBS representation for curves and surfaces allows efficient and accurate modelling of complex geometries, in which the detail of one surface does not affect the detail of other surfaces. The hierarchical arrangement of surfaces governs the interaction between surfaces, and allows the modeller to switch on and off specific heterogeneities or length scales, because no modification of an underlying grid is needed. Consequently, the computational efficiency of NURBS representation makes it possible to generate quickly multiple realisations of surface-based models, and to test the impact of different geological scenarios.
Because NURBS surfaces are computationally efficient, generating multiple realisations is not time-consuming. For example, the shoreface parasequences and channelised fluvial sandbodies in the upper and lower parts of the model shown in Fig. 11 take, respectively, 12 and 23 s to build on a standard workstation. Exporting .STEP files requires less than 2 s. In contrast, exporting cornerpoint-gridded or triangulated versions of the same NURBS surfaces is at least two orders of magnitude slower, depending on the resolution of discretisation.

Compatibility with Conventional and Next-Generation Reservoir Modelling Workflows
NURBS-surface-based models can be exported in different formats, depending on the workflow design. Our workflow uses an unstructured tetrahedral mesh to preserve the geometry of the modelled surfaces (Melnikova et al. 2016), which can be used directly, without upscaling, for flow simulation using a next-generation, adaptivemesh simulator (e.g. IC-FERST; Jackson et al. 2015;Salinas et al. 2017) or to conduct flow diagnostics (Zhang et al. 2017). However, NURBS surfaces can be exported in discretised formats, such as triangulated surfaces for visualisation and quality control (.TS, .OBJ, .PLY), or as gridded surfaces for compatibility with most commercial reservoir modelling packages (.ZMAP), effectively enabling a hybrid surface-and grid-based approach. Alternatively, surface-based models can be used to generate three-dimensional training images for modelling workflows that use multiple-point statistics (Strebelle and Levy 2008).

Conclusions
Reservoir modelling using parametric non-uniform rational B-splines (NURBS) surfaces and curves has many advantages over conventional grid-based reservoir modelling approaches. A NURBS-surface-based modelling approach allows all geological heterogeneity of interest and wells of any trajectory to be represented without reference to a pre-defined grid. With a limited number of control points, geometrically complex NURBS surfaces and curves can be created, with a high level of detail where needed. As a result, NURBS provide a computationally efficient way to represent geological heterogeneity associated with complex geometries (e.g. non-monotonic surfaces, volumetric pinchouts, non-alignment to pre-existing grid) and occurring over a range of length scales, which are difficult or impossible to represent in conventional grid-based models. Because NURBS form smooth surfaces, no stair-stepping effects related to grid resolution are introduced. The interaction of different NURBS surfaces by joining, terminating and stacking creates "watertight" volumes that are the constituents of a surface-based reservoir model that is ready for meshing and numerical analysis.