Realistic modelling of faults in LoopStructural 1.0

Without properly accounting for both fault kinematics and faulted surface observations, it is challenging to create 3D geological models of faulted geological units that are seen in all tectonic settings. Geometries where multiple faults interact, where the faulted surface geometry significantly deviate from a flat plane and where the geological interfaces are poorly characterised by sparse data sets are particular challenges. There are two existing approaches for incorporating faults into geological surface modelling: one approach incorporates the fault displacement into the surface description but does not 5 incorporate fault kinematics and in most cases will produce geologically unexpected results such as shrinking intrusions, fold hinges without offset and layer thickness growth in flat oblique faults. Another approach builds a continuous surface without faulting and then applies a kinematic fault operator to the continuous surface to create the displacement. Both approaches have their strengths, however neither approach can capture the interaction of faults within complicated fault networks e.g fault duplexes, flower structures and listric faults because they either (1) impose an incorrect (not defined by data) fault slip direction; 10 or (2) require an over sampled data set that describes the faulted surface location. In this study we integrate the fault kinematics into the implicit surface by using the fault kinematic model to restore observations and the model domain prior to interpolating the faulted surface. This approach can build models that are consistent with observations of the faulted surface and fault kinematics. Integrating fault kinematics directly into the implicit surface description allows for complex fault stratigraphy and fault-fault interactions to be modelled. Our approaches show significant improvement in capturing faulted surface geometries 15 especially where the intersection angle between the faulted surface geometry and the fault surface varies (e.g. intrusions, fold series) and when modelling interacting faults (fault duplex).


Introduction
Building 3D geological models of complex geometries (e.g. fault duplexes, flower structures, folds, fold interference patterns) 20 from only field observations is still a challenge for most 3D modelling approaches. This challenge exists because many of the inverse problem. In this probabilistic framework model data points become parameters for the forward model and fault displacement can be incorporated by including a likelihood function where the fault displacement is expected to come from a prior probability distribution, that may be guided by observations of the fault displacement (de la Varga and Wellmann, 2016;Wellmann et al., 2017;Godefroy et al., 2018). However, this approach is still limited by the inability for fault kinematics to be 60 incorporated into the step function.
The alternative approach is to apply the fault displacement to an already interpolated surface. A first possibility is to apply restoration of the geometric structures surrounding the faults and find the displacement field based on some mechanical likelihood criterion (Thibaut et al., 1996;Maerten and Maerten, 2015). However, finding an objective mechanical criterion is difficult and computationally challenging, an alternative is to construct the geometry of the continuous surface prior to faulting 65 and then to apply the fault operator to the faulted surface (Godefroy et al., 2018;Laurent et al., 2013;Georgsen et al., 2012).
In this approach, both the slip direction and the magnitude are used to constrain the fault displacement. Since a continuous surface needs to be interpolated prior to applying the fault, the surface observations within the fault ellipsoid are discarded when interpolating the initial surface. The observations within the fault ellipsoid can then be used to determine the fault displacement magnitude by minimising the misfit between the surface and observations. 70 Using step functions in the implicit surface description and adding faults to existing surfaces are both useful in basin settings where the pre-fault geometry is reasonably simple and most geometrical complexities are associated with the fault geometry.
However, it is unclear how these approaches can be applied to fault networks with a large number of interconnected discrete faults or when modelling complex poly-deformed geometries. When modelling multiple interconnected faults (e.g. splay faults, crosscutting faults and abutting faults), it can be challenging to properly account for the relative kinematics between faults, 75 especially where the observations of the faulted surfaces are limited to a map view.
In this contribution, we adapt the approach presented by Laurent et al. (2013) and Godefroy et al. (2018), by incorporating the fault operator directly into the implicit function. This uses the fault kinematics to restore the model domain and data points to a pre-fault geometry. The restoration function is defined within a structural frame (Grose et al., 2020) where the coordinates define fault surface, fault slip direction and the fault extent. This restoration function allows for the kinematics of the fault to 80 be directly incorporated into the modelling work flow. Faults are modelled starting with the most recent and the kinematics are added to every consecutive geological feature added to the model. This application of time aware geological modelling follows on from previous work for modelling folds by (Grose et al., 2017;Laurent et al., 2016). In this paper we outline the framework for integrating faults used in LoopStructural (Grose et al., 2020). We show case studies ranging in complexity from a simple faulted intrusion where we demonstrate the limitation of step functions, a thrust duplex system, and a faulted fold series using 85 the fold modelling framework from Laurent et al. (2016) and Grose et al. (2017).

Surface representation
Three dimensional structural geological models are a representation of subsurface geology where the geological units are either represented using boundary surfaces (upper and lower contacts) (Wellmann and Caumon, 2018) or a discretisation on a 90 predefined support. There are two approaches that are commonly used for representing surfaces in 3D geological models. The first approach uses explicit surface representation: the geometry of surfaces are contained using a support that is collocated with the surface geometry. The surfaces are represented using discrete objects such as triangulated surfaces, two-dimensional grids or parametric surfaces. The surface geometry is usually built by either directly triangulating data points or using interpolation algorithms to create a smooth surface fitting the data (Chilès and Delfiner, 1999;de Kemp, 1999;Mallet, 1992). The explicit 95 surface representation means that the geometry of the surface is only represented where the surface is located.
The second approach, implicit surface representation, uses isovalues of one or several scalar fields in 3D space to represent the geometry of geological surfaces (e.g. stratigraphic horizons and fault surfaces) (Jessell, 1981;Calcagno et al., 2008;Caumon et al., 2013;Cowan et al., 2003;Frank et al., 2007;Hillier et al., 2014;Lajaunie et al., 1997;Mallet, 2002Mallet, , 2014Maxelon et al., 2009;Moyen et al., 2004;Yang et al., 2019;Renaudeau et al., 2019;Manchuk and Deutsch, 2019;Gonçalves et al., 2017;100 De La Varga et al., 2019). The value of the scalar field represents the distance away from a reference horizon. Alternatively, if the geological interfaces are conformable, the surfaces can be represented by isovalues of the scalar field representing the relative thicknesses between the interfaces. The gradient of the scalar field is a representation of the orientation of the surface being modelled. These scalar fields can be constructed using a variety of different interpolation methods: e.g. co-kriging (Calcagno et al., 2008;De La Varga et al., 2019;Gonçalves et al., 2017;Lajaunie et al., 1997), radial basis functions (Cowan et al., 2003;105 Hillier et al., 2014) or using discrete interpolation on a predefined support (Irakarama et al., 2020;Caumon et al., 2013;Frank et al., 2007;Moyen et al., 2004;Mallet, 1992). Most modelling approaches use the same interpolation algorithms for all types of geological surfaces even though the physical properties and expected geometries can vary.
These interpolation methods mean models can be created with less bias and are more reproducible. Implicit methods also make it possible to generate a set of model realisations reflecting geological uncertainty (Manchuk and Deutsch, 2019;110 Cherpeau et al., 2010;Henrion et al., 2010;Wellmann et al., 2010;Lindsay et al., 2012;Yang et al., 2019).

Faults in implicit modelling
There are a number of existing techniques available for incorporating faults into the implicit approach with methods allowing for topology (Calcagno et al., 2008;Moyen et al., 2004), kinematics (Laurent et al., 2013;Godefroy et al., 2018) and uncertainty (Cherpeau and Caumon, 2015;Cherpeau et al., 2012) to be incorporated into the modelling workflow. Incorporating faults into 115 implicit surface representation is difficult because faults represent a discontinuity in geological feature being modelled and require a discontinuity in the implicit function. The geometry of the faulted geological feature (stratigraphic interfaces, faults or foliations) should be consistent with the kinematics of the fault.
There are three possible approaches that can be used to incorporate faults into implicit surface description: (1) interpolate fault domains using independent implicit functions; (2) incorporate the fault into the domain discretization or; (3) apply a fault 120 operator to an already interpolated surface.
The first approach separates the faulted feature using two implicit functions, one on the hanging wall and another on the foot wall. A boolean operation using the fault surface geometry is used to define which implicit function is used to represent the faulted feature (Caumon et al., 2013;Cowan et al., 2003;Frank et al., 2007). When the hanging wall and foot wall are represented using separate implicit functions, the displacement of the fault is not explicitly defined and will be a result of The second approach directly incorporates the location of the fault into the polynomial trend of the implicit function using a discontinuous step function (Calcagno et al., 2008;De La Varga et al., 2019;Marechal, 1984). The fault displacement magnitude is not used as an observation for building the model but is found by solving the interpolation system for the cokriging matrix. de la Varga and Wellmann (2016) applied Bayesian inference to incorporate the fault displacement into the 135 implicit modelling work flow by using the observations as parameters for the inverse problem. In their approach models with different displacements are generated by perturbing the observations of the faulted surface. The objective function that is incorporated into their Bayesian model as a likelihood, only considers the magnitude of the displacement and not the direction of displacement. This approach provides an interesting idea of using Bayesian inference as a method for optimising the model parameters in order to try and generate models which fit both geological knowledge and geological observations. An alternative 140 to this approach for discrete implicit methods is to introduce a discontinuity in the regularisation term, by either cutting the interpolation support by the fault geometry (Caumon et al., 2013;Frank et al., 2007) or introducing ghost nodes into the discrete elements (Irakarama et al., 2020).
The third approach integrates the fault kinematics into implicit 3D modelling work flows by using a numerical operator that defines the distribution of fault slip on and close to the fault surface (Godefroy et al., 2017(Godefroy et al., , 2018Laurent et al., 2013). In these 145 approaches a fault operator models fault displacement within an ellipsoidal fault volume. The displacement is discontinuous across the fault surface and decreases smoothly away from the fault. In this approach, a fault volume is defined where the fault displacement is non-zero. The geological interfaces are defined by first interpolating a continuous surface across the fault. Displacement profiles can be identified, either manually defined or using numerical optimisation techniques (Godefroy et al., 2018). This approach applies the fault to a surface and is not integrated with the implicit function meaning it would not 150 directly integrate with approaches where the implicit function represents a continuous features such as foliations, conformable stratigraphy or grade shells.
Alternative kinematic models for describing fault displacements exist such as the trishear model (Erslev, 1991;Allmendinger, 1998). Trishear is a kinematic model describing the geometry of fault propagation folds and has been presented in pseudo 3D (Cristallini and Allmendinger, 2001) and in true 3D (Cardozo, 2008). Trishear has been used to understand 155 the structural evolution that has resulting in the geometries observed in the field and as part of seismic interpretations. In the interpretation of reflection seismic data, kinematic fault restoration approaches based on vector fields have been proposed (Hale, 2013;Wu et al., 2015). For example, Wu et al. (2015) determine the fault tip by seismic image processing, then interpolate fault shifts in the volume by least-squares smoothing. However, this model is most directly applicable in contexts where observations are sparse, and the smoothing assumes that fault shifts depend on the stratigraphic orientation. It also 160 restores all faults at once whereas sequential restorations allow for the fault topology to be includes.

A kinematic framework for modelling faults
Faults are described by geologists by the movement of rock packages and not only by the location of the discontinuity. Structural geologists describe faults with the fault slip and the displacement magnitude. In a lot of cases measuring the fault slip direction is impossible, the geologist will have a conceptual model to describe the fault displacement direction, for example 175 using their conceptual model based on geological knowledge. This means that to properly capture the displacement kinematics in structural models the additional kinematic constraints need to be encoded into the modelling algorithms.
In this contribution, the work flows presented by Laurent et al. (2013) and Godefroy et al. (2018) are modified by incorporating the fault operator into the implicit representation of the geological feature in the model (Fig. 2). The faulted stratigraphy can then be interpolated as a continuous layer where the data points are restored using the fault kinematics. The fault kinematics 180 are integrated into the implicit surface description, meaning that the scalar field when evaluated in the model domain captures the fault kinematics. Integrating the fault kinematics into the implicit modelling system is achieved by: 1. Building the fault frame, a curvilinear coordinate system representing the fault geometry 2. Defining the fault displacement within the model domain  3. Adding the fault kinematics to the implicit surface description 185 3.1 Fault frame A curvilinear coordinate system (the fault frame, described empirically for natural faults by Walsh and Watterson (1987)) represent different geometrical elements of the fault geometry and is built by interpolating structural observations. The first coordinate ( f 0 ) represents the distance to the fault surface and the isovalue of 0 will represent the fault surface and can be interpolated from any structural observations of the fault surface e.g. strike and dip controls and observations of the fault 190 surface location. The second coordinate ( f 1 ) measures the displacement direction of the fault and the normal to the gradient of this field will be parallel to any kinematic indicators for the fault (e.g. slickenslides, stretching lineations) and parallel to the fault surface. The third coordinate ( f 2 ) measures the distance in the direction of the fault extent and the gradient of this field will be orthogonal to the gradient norm of the fault surface ( f 0 ) and also orthogonal to the fault displacement field ( f 1 ).
Three local direction vectors are implicitly defined by the normalised gradient of the fault frame for any location: The fault frame can then be queried for any location within the model returning the distance to the fault centre and the fault frame vectors.

Fault displacement
The displacement direction of the fault is represented by f 1 . The fault displacement magnitude can be either defined by a 205 constant displacement on the hanging wall or by a function of the fault frame coordinates allowing for the fault displacement to be represented relative to a fault origin using a combination of attenuation profiles as described by Laurent et al. (2013) and Godefroy et al. (2018). We define the local fault displacement vector for any location as where d is the scalar value of the fault displacement magnitude at that location. should be related back to the geological topology, e.g. the overprinting relationships between geological events.
To represent a faulted geological feature we first define the geological feature in the model using the associated structural data. A fault can then be applied to this geological feature which produces a faulted geological feature. The fault restores the observations attached to the faulted geological feature to their location using the fault kinematics from the fault operator.
The faulted geological feature can then be interpolated using the restored locations. To evaluate the scalar field of the faulted 220 geological feature the fault transformation must be added to the implicit surface representation. This is done by applying the fault restoration kinematics to the data points prior to evaluating the scalar field values.
Using this work flow it is possible to stack multiple fault operators and apply, backwards in time, the appropriate displacements to restore the location to its position prior to faulting. This approach can be applied to any structural elements in the model, 225 including lineations, fold vergence, tectonic foliations and fold frames (Grose et al., 2017).

Implementation
The generalised work flow presented above does not depend on the specifics of the interpolation scheme. We present an overview of how to build the fault frame using discrete interpolation algorithms using LoopStructural (Grose et al., 2020), a open source 3D modelling library 1 . Using discrete interpolation algorithms, the complexity of the interpolation problem is 230 defined by the resolution of the support and not by the number of constraints.

Building fault frame
The fault frame is built by interpolating the scalar fields using the following procedure: 1. Interpolate coordinate 0 to represent the geometry of the fault surface so that the isosurface of 0 contains the fault trace and the field is parallel to the orientation of the fault surface. 3. Interpolate coordinate 2 so that the direction of its gradient norm is orthogonal to the direction of the gradient norm of the fault surface field and to the fault slip direction field.
In this study we use a discrete implicit modelling technique where the implicit function is represented on a discrete support 240 by either a piece-wise linear function on a tetrahedral mesh (Mallet, 2002;Frank et al., 2007;Caumon et al., 2013;Grose et al., 2020) or by using finite differences on a Cartesian grid where an objective function such as the bending energy (Renaudeau et al., 2019) or the second derivative is minimised (Irakarama et al., 2020). field (∇ f j ) so that the gradient norm of the new field is orthogonal to the gradient norm of the existing field.
The scalar field can then be interpolated by minimising the misfit between the orthogonality constraint and any observations of the fault slip direction.

250
The element orthogonality constraint can be applied to two fields and the resulting scalar field will be close to orthogonal to both. Frequently it is difficult to measure or observe the fault slip direction. In these cases the geologists would need to use qualitative geological knowledge to define the fault slip direction using their knowledge about the geology of the area. It is also possible to build the model with different fault slip directions (and magnitudes) in order to explore the range in possible geometries and compare these with observations.  δ(X) = D 0 ( f 0 (X)) · D 1 ( f 1 (X)) · D 2 ( f 2 (X)) (5)

Splay faults
Faults are rarely single surfaces and are often defined as fault zones with complex linkage between individual fault segments.
During the development of a fault zone, an individual surface can splay meaning that the displacement of the fault is separated onto multiple slip surfaces (Huggins et al., 1995;Marchal et al., 2003). The resulting geometries can be increasingly complex and includes duplexes and flower structures.

270
An advantage of our approach, is that multiple faults can be applied sequentially where the total displacement is defined by the accumulation of the local displacement fields. Modelling faults in a time aware approach is similar to how folds are modelled (Laurent et al., 2016;Grose et al., 2017Grose et al., , 2018 where the most recent fold is modelled first and this defines the geometry of the folded foliations.
A fault zone can be represented by a piece-wise combination of fault segments. To model multiple connected faults the fault 275 frames need to equal along the branch line. In the following section, hard constraints are incorporated into the interpolation system using using Lagrange multipliers.
Discrete interpolation approaches define the interpolation within a volume of interest by a linear system where the unknowns are the property values on the mesh nodes (Frank et al., 2007;Caumon et al., 2013). The unknowns are found by solving the linear system: Where A is an interpolation matrix that contains the different geological constraints (contact observations, orientation constraints, orthogonality constraints) and a regularisation term that controls the smoothness of the resulting scalar field and x are the values of the scalar field at the interpolation nodes. The linear system defined in equation 6 is over determined and can be solved in a least squares sense by minimising the L-2 norm ||Ax − b|| 2 2 .

285
Splay faults (Fig. 3a) share the same geometry as the parent fault at the point of intersection. To ensure the faults share the same geometry an additional constraint needs to be added to the linear system so that the scalar field value of the shared nodes of the fault are equal. In Fig. 3b the main fault is interpolated using all interpolation nodes. To interpolate the splay faults the hard constraints need to be identified by finding the values of the fault extent field ( f 2 ) where the fault is active. Equality constraints are identified by finding the elements that contain the isosurface of the fault extent field containing the intersection, 290 point e.g., black nodes in Fig. 3c. The hard constraints form an an additional linear system, where C is the interpolation matrix for the equality constraints.
To solve this system ensuring that the hard constraints are strictly honoured and the other interpolation constraints are satisfied in a least squares sense we define the Lagrangian for the least squares problem, where λ are the Lagrange multipliers. To apply this constraint we need to know where the two faults share geometries so that the values of the main fault can be used as equality constraints to interpolate the splay fault. Where the geometry of the two surfaces is shared, the property value of the nodes inside the shared region are fixed so that both surfaces have the same value ( Fig. 3c). This can be implemented where C i, j = 1 where i is index of hard constraint and j is the node global index. d j is equal 300 to the shared node value between the two scalar fields. Alternatively, the hard constraints could be defined using the discrete support shape functions (e.g. linear tetrahedron or trilinear interpolation), however in this case the line of intersection would need to be known, which may be complicated or impossible to quantify.
In Fig. 4B and C, two splay faults are interpolated, in B) the same constraints are used to solve both surfaces and the resulting faults overlap but do not share the same geometry in the region where they are expected to splay. In B) equality constraints 305 are used to constrain the surface geometry to be the same before the fault splay. When modelling complex fault networks for example splay faults or other faults with shared geometries the hard constraints ensures that the faults have the same geometry at intersection lines. For splay faults this will also ensure that the interpolation of the splay fault also includes the parent fault geometry. In the following sections we will present three synthetic case studies that demonstrate the application of our approach to both extensional and contractional systems. All of the case studies have been implemented in the open source 3D structual modelling library, LoopStructural (Grose et al., 2020). The case studies are provided as interactive Jupyter notebooks in Appendix A.

Faulted intrusion
The first case study is a synthetic intrusive body that has been faulted by a planar normal fault. The intrusion has been randomly 315 sampled with one set of data points describing the distance to the intrusion surface, where the value is 0 and another set of data points where the value is 1 which describes the distance to the intrusion surface. In Fig. 5A  the fault displacement magnitude. In this example a fault displacement of 500m is used. Fig. 5B shows the fault displacement vector field applied to the data points on the hanging wall of the fault indicating the direction of movement to restore the data points to their location prior to faulting. The data points that are identified as being on the footwall of the fault are not affected by the fault restoration and are grey coloured without a vector field. Fig. 5C shows the data points from Fig. 5A in the restored locations after applying the fault operator. The intrusion shape is interpolated within the restored space Fig. 5C,D using a 325 discrete interpolation approach with a weak regularisation constraint to allow for the propagation of a high curvature geometry.
A weak regularisation constraint is applied because modelling intrusions with implicit interpolation approaches has similar  (Laurent et al., 2016), where high curvature is expected but the interpolation actually minimises the curvature. A kinematic operator is added to the interpolated scalar field so that when evaluating the scalar field at a location, the location is restored using the fault displacement field. Fig. 5E shows the model support for a regular grid where the points in 330 the hanging wall are restored using the restoration vector field. In practice, the points are not actually moved and the restoration vector field is actually added to the implicit surface description resulting in a faulted intrusion Fig. 5F.
In this example the intrusion volume is maintained when the fault is applied. The volume of an intrusion when modelled using step functions is likely to change as shown in Fig. 1.

350
Different profiles can be chosen either to fit observations or to enforce a particular conceptual model.
The resulting 3D displacement scalar field is shown in Fig. 6A, where displacement scales from 0 to ±1 close to the fault.
This displacement field and the fault slip direction vector field can be added to the implicit surface scalar field (shown in Fig. 6C). Fig. 6D shows the faulted stratigraphic surfaces where the implicit function is evaluated using the fault restoration function. In Fig. 6B and D the displacement of the fault decreases away from the fault and at depth with the displacement not

Faulted fold series
In the following example we model a faulted doubly plunging fold. The reference model was generated using Noddy (Jessell and Valenta, 1996) and the orientation of the fault surface, fold axial surfaces and folded surface were sampled at irregular locations shown in the geological map (Fig. 8).

380
To model a folded and refolded structures we use the framework introduced by Laurent et al. (2016) and Grose et al. (2018).
To model the faulted fold series the most recent structure needs to be modelled first. In this case this is the fault as it overprints the fold series. We The fold frame for F1 is interpolated in a similar way to the fault frame, by first interpolating the axial foliation and then interpolating the fold axis direction field (Laurent et al., 2016). The fault is added to the fold frame meaning that the interpolation of the fold frame occurs after the observations have been restored. The scalar fields for the fold frame can be seen in 390 Fig. 10A and B.
Within the fold frame, the intersection lineation is calculated by finding the cross product between the normal to observations of bedding and the gradient of the interpolated axial foliation. The fold axis rotation angle is calculated by finding the angle between the intersection lineation and the fold axis direction field. The fold axis rotation angle is interpolated throughout the model by fitting a Fourier series to the calculated fold rotation angles. We use the Levenberg-Marquardt algorithm for fitting the wavelength is slightly larger (blue curve in Fig. 9), corresponding to a more consistent alignment of the fold hinges across A single fault is challenging to incorporate into 3D geological modelling because it defines a discontinuity in the older geological surfaces. This challenge is further complicated because the fault is not only defined by the location of the discontinuity 420 but by the kinematics associated with the displacement of the faulted objects. This means that to properly model faults, the kinematics need to be directly incorporated into the modelling process. These complications are further exaggerated when modelling a fault network where there are multiple overprinting and interacting faults.
In this contribution faults are modelled using a curvilinear coordinate system and fault operator adapted from Laurent et al. (2013) and Godefroy et al. (2018) to directly incorporate fault kinematics into the surface description. Our approach 425 applies the fault operator in reverse prior to interpolating the surface, meaning that the fault description is incorporated into the implicit surface representation rather than applied to an already discretised surface. Both approaches require building a reference coordinate system for the fault that defines the fault surface, fault slip direction and the fault extent. The coordinate system can be easily constrained from observations where the fault slip direction must always be orthogonal to the fault surface and the fault extent is always orthogonal to the fault slip direction and the fault surface. We show that this coordinate system 430 can be built allowing for multiple overprinting and interacting faults to be modelled using a serial approach where the scalar fields representing the coordinates are modelled sequentially. While this approach produces an acceptable coordinate system for the faulting, it can be costly to evaluate especially on high resolution models. The following strategies could be applied for modelling large fault networks: (1) Lower the resolution of the interpolation support for the fault network, as the geometry of the fault should be less complicated than stratigraphy.
(2) Only model the fault geometry within the fault volume.

435
One of the key advantages of our approach is the ability to make use of kinematic observations for constraining the geometry of faulted geological surfaces. However, in a lot of cases (e.g. geological survey maps, potential field geophysical interpretations), the required kinematic indicators are unknown. One approach for solving this problem would be to build the fault frame by aligning the fault surface with observations of the fault surface (this step is no different to any other implicit modelling approach) and second to align the fault extent direction with the map extent of the fault. This defines the fault volume as an 440 ellipsoid where the major axis is aligned with the map pattern of the fault. The fault slip direction is then orthogonal to both of these interpolated fields. This approach would provide a good first pass approach for modelling the fault displacement and is a similar outcome to the inferred slip direction used by step functions, albeit using our approach the fault kinematics are modelled appropriate to the dip of the fault.
One of the challenges for the presented approach is determining the correct displacement magnitude and profile. Although 445 in some cases these parameters are observations recorded on a geological map they are not always present. Godefroy et al.
(2018) applies particle swarm optimisation for these parameters by optimising the displacement to the data points within the fault volume. In a similar approach, an objective function could be applied to the restored continuous surface minimising some metric, for example for folding this could be incorporated into the probabilistic geological inversion presented by Grose et al. (2018Grose et al. ( , 2019. In this framework, fitting a fold model to the S-Plot could be used to determine whether a displacement parameter 450 and fault profile are acceptable. Fig. 9C shows the fold limb rotation angle for the faulted fold series where the fault has been added (blue dots) and where the fault hasn't been added (red dots). Adding the fault changes the shape of the S-Plot, and could be added as a likelihood function in an inversion framework. As with all geological problems this would be highly uncertain and require additional domain knowledge to result at a solution.

455
The geometry of fault networks and the displacement of faults on other geological surfaces is challenging to integrate into implicit 3D modelling algorithms. Directly integrating faults into the implicit description of surfaces using step functions does not allow for kinematic observations to be incorporated and will result in fault geometries with inconsistent kinematics. In order to mask these flaws, it is possible to constraint fault geometries using interpretive constraints. In this contribution the kinematics are applied in reverse to restore the observation into a pre-fault location. Within the restored model space a continuous surface