pynoddy 1.0: an experimental platform for automated 3-D kinematic and potential field modelling

Introduction Conclusions References


Introduction
A wide range of methods exists for the computational synthesis of geological models as interpreta-20 tions about the structure of the subsurface(see, for example, Jessell et al., 2014, for a recent overview of methods). Each modelling method focusses on different aspects of geological data and concepts, but they can be broadly classified in terms of: (1) surface-or volume-based interpolation techniques, (2) pure geophysical inversions, and (3) mechanical or kinematic modelling approaches. We present here a set of open-source Python modules for the efficient, flexible and reproducible construction of 25 kinematic structural models to enable the analysis of uncertainties in geological models.
Structural geological models are generally produced by combining information from direct observations (e.g. measurements in outcrops or boreholes) and indirect data, for example interpreted from geophysical data. Additional aspects of the conceptual geological model or the structural setting are, in the general case, only indirectly taken into account. Computational methods, which are able to 30 capture several or all of the previous considerations, are then used to produce the model.
Regardless of the approach taken, the resulting models always contain uncertainties. These uncertainties are increasingly recognised (Bond et al., 2007;Caers, 2011;Bond, 2015) and addressed with novel methods for uncertainty analysis and visualisation (e.g. Bistacchi et al., 2008;Suzuki et al., 2008;Jessell et al., 2010;Polson and Curtis, 2010;Wellmann et al., 2010;Lindsay et al., 2012;35 Cherpeau et al., 2012;Lindsay et al., 2013;Laurent et al., 2015). So far, the analysis of uncertainties in kinematic models has been performed for balanced cross sections (e.g. Judge and Allmendinger, 2011) and detailed fault displacement models (Laurent et al., 2013). We contribute here with a methods to analyse and visualise uncertainties in automatically constructed kinematic forward models.
To enable this functionality, we extend the capability of an existing kinematic modelling method, 40 implemented in the software Noddy (Jessell, 1981;Jessell and Valenta, 1996), with a flexible set of dedicated scripting modules developed in the programming language Python. Our aim is to provide high-level access to the underlying model construction methods, enabling: (1) flexible and rapid construction of kinematic models; (2) the definition of fully reproducible modelling experiments, and (3) a framework for automatic model generation, to enable experiments and analyses that require

Materials and Methods
Because we extend the functionality of an existing kinematic modelling package, Noddy (Jessell, 1981;Jessell and Valenta, 1996), we briefly describe its functionality here, and then provide details about the implementation of the Python package we have developed, referred to hereafter as 60 pynoddy. Finally, in order to describe the main difference between our approach and other commonly used structural interpolation methods, we also briefly review the relevant approaches in this direction.

Structural geological modelling concepts
Structural geological models can be constructed with different approaches, and the choice of a spe-65 cific modelling method directly depends on the model applications and the available input information.
The approach that we apply here is based on kinematic modelling concepts. The distinction between interpolation and kinematic methods is most apparent when considering the types of data and geological constraints that are honoured. The most common approach to construct structural 70 models is based on surface and volume interpolation methods (Mallet, 1992;Lajaunie et al., 1997;Sprague et al., 2006;Caumon et al., 2009;Hillier et al., 2014;Jessell et al., 2014). An example of the general interpolation function is presented in Fig. 1a. Structural interpolations focus on honouring parameterised surface contact points (Caumon et al., 2009), although secondary data like orientation measurements can also be taken into account (Lajaunie et al., 1997;Calcagno et al., 2008;Hillier 75 et al., 2014). Constraints on the shape of geological surfaces, or the interaction with other units or faults, are then assigned to different surfaces, according to observations in the field or the expected geological settings. While these considerations are clearly based on geological reasoning, it is not guaranteed that an interpolated structural model matches all the known aspects of the geological evolution of an area. For example, it is easily possible that constraints on thickness of geological 80 units are not consistent, for example across a fault, leading to a violation of mass conservation. Additionally, a wide range of structures observed in multiply deformed terranes, such as complex fault networks or refolded folds, are difficult to construct consistently using current interpolation methods.
Another end member in the evaluation of the structural setting are simulations of physical processes (e.g. Gerya and Yuen, 2007;Moresi et al., 2007;Kaus et al., 2008;Regenauer-Lieb et al., 85 2013). Instead of starting with geological observations, these methods are based on mathematical models capturing relevant physical processes that led to the formation of specific structures (Fig. 1b).
For realistic simulations, meaningful constitutive models and boundary conditions are required. Multiple different methods exist which capture different aspects of the mechanical deformation, and more and more commonly also the effect of coupled Thermo-,  simulations. However, these types of simulations are not yet commonly applied to model the entire  (c) kinematic models (modified from Jessell and Valenta, 1996) complexity of multiply deformed geological regions as simulations are computationally demanding and rock properties, and boundary conditions are not always perfectly known. Furthermore, they require an initial distribution of rock properties in space as initial conditions, often determined from an explicit or implicit interpolation approach.

95
Kinematic modelling methods focus on major tectonic and metamorphic events in geological history (Jessell and Valenta, 1996) and are therefore conceptually located between the previously described end members (Fig. 1c). In these models, the complexity of deformation is greatly reduced and captured in simplified kinematic functions as surrogate models. This means that direct geological observations of surface contacts and orientation measurements are not taken into account in 100 the simulation step. However, the simulations are very fast and enable therefore a quick testing of different deformational scenarios, and the interaction of multiple events in geological history. Furthermore, rapid simulation makes direct (and ideally automated) comparisons between the model and observed structures feasible, allowing the indirect incorporation of geological observations. We will present several examples in which this trade-off between physical realism and geological ob-105 servations can lead to useful insights into the interaction and relevance of deformational events in geological history.

Kinematic structural modelling with Noddy
Noddy models begin as a layer cake stratigraphy, for which the heights of the stratigraphic contacts and geophysical rock properties are defined. A history of relevant events that affected the model 110 region is then developed from a predefined set of events, including: folds, faults, and shear zones; unconformities, dykes and igneous plugs; regional tilting and homogeneous strain. In addition to modifying the initial stratigraphy, each event can define (geophysical) alteration halos, penetrative cleavages and lineations.
Each Noddy event is defined by four classes of properties: form, position, orientation, and scale 115 (Jessell, 1981;Jessell and Valenta, 1996). For example, a fault is defined by its dip and dip direction, the pitch and magnitude of the slip vector, and the position of one point on its surface (note that more complex definitions of the fault plane are also possible, c.f. Jessell, 1981). The use of geological descriptions provides a natural and intuitive framework for geologists to build a model. Even though the structural events themselves are relatively simple, complex geometries quickly develop as two Displacement equations are stored as a "history", which provides parameterised definitions of the model kinematics and rock properties. A voxel model of any 3D rectangular volume of interest can be calculated from this history by considering each voxel independently using the Eulerian (inverse) form of the defining Lagrangian displacement equations, and applying them in reverse chronological 125 order (i.e. starting with the most recent deformation event). This operation transforms the x, y, and z position of each voxel into the x, y, z position at the time the associated volume of rock was created.
The properties of this voxel can then be calculated directly from the base stratigraphy.
New lithologies can also be created during three specific event types: unconformities, dykes, and plugs. These events are assumed to be instantaneous, and are ordered relative to other events. In 130 order to simplify the underlying kinematic equations, they are all defined in a standard reference frame that is orthogonal to the symmetry of the deformation event. The real world reference frame is rotated into the standard reference frame prior to the calculation of each event, and then subsequently rotated back to the real world reference frame using the variations in the z-values as a continuous implicit field that can be iso-surfaced to produce stratigraphic horizons.

135
As well as the initial position of the point, a binary "discontinuity code" is stored, that records each time a voxel is affected by an event described by a discontinuous displacement equation (faults, unconformities, dykes, and plugs) but ignores events described by continuous displacement equations (folds, shear zones, strain, rotation, foliations and lineations). This discontinuity code allows the accurate transformation of the voxel data set to a vector data set, since only voxels which have 140 exactly the same sequence of discontinuity codes are part of the same contiguous volume of rock. If two adjacent voxels have different codes, the difference in the discontinuity code that occurred most recently defines the discontinuity which separates them.  Similarly, the orientation of a linear feature is calculated from the intersection of two planes.
Thus both the Eulerian (inverse) and Lagrangian (forward) descriptions of the displacements must be available for a new deformation event to be included in the modelling scheme. For this reason, the displacement equations governing each Noddy deformation event are kept as simple as possible, and superimposed deformation events are combined to produce structural complexity. A full description 155 of the Noddy implementation is presented in Jessell and Valenta (1996).

Basic concept
Noddy has the capability to calculate potential field responses of gravity and magnetics based on the spatial distribution of rock properties in the subsurface, and this functionality is exposed in pynoddy. 160 The petrophysical rock properties of a specific volume are defined by their original stratigraphic value, unless a specific deformation event (faults, unconformities, plugs and dykes) has an associated alteration/metamorphic character, with allows the modification or replacement of pre-existing properties based on that locations distance from the structural feature at the time of the activity of the event. A further complication is possible if a model with an anisotropic magnetic susceptibility 165 (a tensor property) or magnetic remanence (a vector property) is defined, in which case there is the possibility of calculating the voxel-level reorientation of these properties as a result of deformation (for example having a remanence vector deformed during a folding event).
For all surveys the rock property of a cube is defined as the value at the centre of the cube, and for grid surveys (that is, not arbitrary surveys or borehole surveys) the field strength is calculated 170 at the x,y location above the centre of each cube. The Total Magnetic Intensity value calculated for all schemes is actually the value projected onto the Earth's field, following the convention of many modelling schemes. The gravity field calculated is for the z component only.
Three geophysical computational schemes are available in Noddy to calculate magnetic and gravity potential fields. The criteria as to which scheme should be used depends on required accuracy, 175 speed and the various geological situations being modelled. A brief description of each scheme is provided below.

Spatial convolution scheme
The spatial convolution scheme works by calculating the summed response of all the cubes within a cylinder centred on the sensor, with a radius defined by the spatial range term. The calculation 180 for each cube is based on the analytical solution for a dipping prism presented by Hjelt (1972) and Hjelt (1974). In order to calculate solutions near the edge of a block, extra geology is used to produce a padding zone around the block equal in width to the spatial range, so that there are no edge effects in this scheme. The scheme only provides exact solutions when the range is larger than the length of the model. For reasonably complex geology this limitation does not result in inaccurate 185 models, however for idealised geometries using a range that is too small results in a kink in resultant profiles. The spatial convolution scheme is slower than the Spectral scheme for medium ranges (10-20 cube ranges), but generally much faster than the Full Spatial Calculation. As long as the range is greater than the spacing between high density/susceptibility features, the inaccuracies associated with truncating the calculation is probably not evident. The draped survey and down-hole surveys 190 have not been implemented for this scheme.

Spectral scheme
This scheme, based on pioneering work by Parker (1972) works by transforming the rock property distributions into the Fourier domain, applying a transformed convolution, and then transforming this result back into the Spatial Domain. The calculation is performed for each horizontal slice through 195 the geology, and the results are summed vertically. The Spectral scheme produces a different result than the other two schemes in terms of absolute numbers for three reasons: 1. The Fourier transform implies that the geology is infinitely repeating outside the calculation area. This produces edge effects when high susceptibility or density bodies are found near the edges of the survey area. This effect can be lessened by the choice of a suitable padding around 200 the block, including over specified areas of interest, however it cannot be totally removed.
2. The calculation loses the absolute base line of the gravity or magnetic field, so even when comparisons are made for well-padded Spectral and large range Spatial models, an overall offset is apparent between the two schemes. When trying to model real data this offset is not a problem as any regional is removed before the modelling process. 3. There is a high frequency component to the calculated field that is of the same wavelength as the cube size and especially apparent when there are steep gradients in the values of the rock properties.

Full spatial scheme
This is similar to the Spatial Convolution scheme except that all the cubes in the model are summed 210 using the Hjelt schemes in order to calculate the response at any point. It generally takes significantly longer to apply this calculation scheme than either of the other schemes. The only exception is when there is a relatively sparse geological model, in which case contiguous blocks with identical petrophysical properties are aggregated to form rectangular blocks, which reduces computation time.
In the extreme case where only one cube has non-zero values for both density and susceptibility, any 215 cubes which have both zero density and susceptibility are ignored. This is the only scheme that can accurately calculate draped surveys, down-hole surveys and arbitrarily located airborne surveys.

Creating input files for kinematic modelling with Noddy
Noddy histories are stored as ASCII files with a simple keyword-value ordering. These files can be written or adapted with any text editor, and the kinematic modelling result computed with a compiled 220 command line version of the program and results visualised with other software.
A graphical user interface (GUI) has previously been created to simplify this model setup, combining convenient input file generation directly with computation and visualisation of the results.
This GUI is freely available (http://tinyurl.com/noddy-site), though currently only runs on Windows operating systems. The GUI version of Noddy is also limited to user-driven workflows, restricting 225 further automation or extension of the methods for scientific experiments.
In order to overcome the problem of either having to work with a direct text input file, or being restricted by the limitations of a GUI, we have developed flexible modules in the programming language Python that enable scripted access the kinematic modelling functionality and to enable the extension to uncertainty estimations.

Implementation of pynoddy
Python is an object-oriented scripting language that is widely used in scientific computation (e.g. Langtangen, 2008). It is highly flexible language, and contains a variety of programming and visualisation libraries ideal for scientific purposes. Python also runs on virtually every operating system that is available, meaning that python wrappers retain the platform independence of C applications.

235
The pynoddy module described here contains a set of classes and functions for managing Noddy input files, passing them to the Noddy command-line application, and processing the results. This approach has many advantages, as it allows automatic generation and analysis of kinematic models in a Python environment, while retaining the performance of Noddy itself (which is written in C).

Overall module structure 240
The package pynoddy contains three main modules: pynoddy.history, pynoddy.output and pynoddy.experiment. The pynoddy.history and pynoddy.output modules provide interfaces for managing Noddy inputs and outputs, while classes defined in pynoddy.experiment provides methods for implementing and performing repeatable modelling experiments. The output of Noddy simulations can be processed and analysed with classes in pynoddy.output, and exported in VTK file formats for

Experiments combining Noddy input and output
If all steps of a pynoddy experiment are automated properly, they can be integrated into one script for model set-up and analysis. This method is leading to a possible reproduction of results (as an example: see the scripts that generate the figures in this manuscript, see appendix B for availability).
This method is often used successfully to ensure reproducibility. It does, however, have one signifi-270 cant drawback: intermediate results or adapted simulation settings have to be stored in separate files and all of those files have to be available to continue with an experiment at a given state.
In order to overcome this limitation, we follow here the aim of including an entire experiment, from the definition of input parameter of the model, to parameters that are specific to an experiment, to the post-processing of results, within a single Python object. Specific experiments can then be 275 defined as child classes inheriting a set of useful base methods. This object can then be stored (for example with a serialisation using the Python pickle package) and retrieved exactly in the state that it was used and defined for a complete reproduction of results, or the adaptation of model parameters to test different model outputs.
The core of the pynoddy.experiment module is the Experiment class, which inherits methods 280 from both the NoddyHistory and NoddyOutput classes, combining and extending their functionality into a single interface that allows a flexible modelling procedure were the Noddy computations are automatically executed when required and outputs directly updated. In addition, methods are provided to encapsulate relevant parameters of an experiment in the most efficient and flexible way.
We consider this last point essential to ensure a full reproducibility of scientific experiments with In order to generate a specific type of experiment, new child classes can then be defined, inheriting from the Experiment base class. Several classes for specific types of experiments are already implemented in the pynoddy package, and we show below the application of one such child class, the UncertaintyAnalysis, applied to a Monte Carlo error propagation experiment.

290
For more details on the implementation and the structure of the modules in pynoddy, please see the documentation and associated source code at the pynoddy GitHub directory (see appendix ?? and B).

Applications
This section outlines the functionality and utility of our pynoddy implementation using a variety of 295 case studies. Firstly, the structural effect of multiple faulting events is investigated, serving mainly as an introduction to the generation of event histories in pynoddy and the visualisation of results.
Then, a model from the Atlas of Structural Geophysics is used to evaluate the sensitivities of calculated gravity potential-field values to changes of parameters in kinematic events. Finally, we use the pynoddy framework to evaluate uncertainties in a case study of the Gippsland Basin, Australia.

Analysis of fault interactions
We start here with an example that is conceptually simple, but can quickly lead to complex structural settings: the interaction of a sequence of fault events on a predefined stratigraphy (Fig. 3a). A more detailed description and interactive version is available as an IPython notebook as part of the repository and as supplementary online material for this manuscript (see Appendix ?? and B).

305
This model is constructed from a stratigraphic sequence containing five units, each 1000 m thick.
We consider a model domain of 10,000 x 7,000 x 5,000 m in x, y, and z-directions. In the following descriptions, we define points with respect to an origin in the model at the top, SW corner (i.e.: the point (0,0,-1000) is at a depth of 1000 m at the SW corner). A representation of the model in a (x,z)-section is given in Fig. 3a.  Fig. 3b. The third event is also a fault, defined with a surface at position (8000, 3500, 0), dipping 60 → 270 and a slip of 1,000 m (Fig. 3c).

315
In terms of this definition of kinematic equations, the two fault events are symmetrical. However, the combination of both events leads, as can be expected, to a non-symmetrical interaction pattern, here clearly visible in the central part of the model (Fig. 3d). The previous example is included to present the possibilities for the simple construction of a kinematic model from start. The model itself is mostly interesting from an instruction or teaching 320 perspective and we will move to more complex models in the following.

Potential field modelling and the Atlas of Structural Geophysics
One motivation for the development of Noddy was to provide a method to explain and teach the effect of subsequent geological events, as we presented an example above. The capability of Noddy to calculate geophysical fields can furthermore be used to provide insights for the interpretation of 325 geophysical potential field data. We can, for example, quickly evaluate how changing the properties of a geological event (for example the dip angle of a fault) influences a simulated potential field.
In fact, this capability of Noddy has been a main driver to develop the "Atlas of Structural Geophysics", an online collection of geological models with their simulated corresponding potential fields for a wide variety of typical structural geological settings (http://tectonique.net/asg). 330 We provide in pynoddy the functionality to directly load models from this atlas into python objects, for further testing and manipulation. In addition, the pynoddy.output module also contains a class definition to read in the calculated potential field responses (NoddyGeophysics). In combina-tion, these methods enable us to quickly test the effect of different event properties on calculated potential fields.

335
As an example, we evaluate here how changing properties of deformational events affects the forward calculated gravity field with a model of a fold and thrust belt (Fig. 4). The required commands to download a model from the web page, to adjust cube size (for better representation), to write it to a file, and to run the model, are combined in a tutorial notebook for detailed reference (see Appendix B). The 3-D visualisation in Fig. 3c was generated through the pynoddy export to VTK and  We calculate the gravity field for this model with the spectral scheme (Sec. 2.3.3) by calling pynoddy.compute in the geophysics simulation mode. The resulting z-component of the gravity field is visualised in Fig. 5a.
As a next step, we evaluate how the effect of a different wavelength in the folding event, as the 345 latest event in the model history, affects the calculated gravity field. This adaptation, as well as the recalculation and visualisation of the geophysical field (Fig. 5b), can be performed with a few lines of Python code (see tutorial notebook for details). In addition, we use simple Python commands to calculate and visualise the difference between the gravity fields of the original and the changed model (Fig. 5c).  With the following example, we now want to highlight an essential advantage of our new implementation in pynoddy: the high-level definition of scientific experiments with kinematic models.

Reproducible Experiments with pynoddy
One main motivation for the definition of a python package to access the functionality of kinematic If all steps of a pynoddy experiment are automated properly, they can be integrated into one script for model set-up and analysis. If implemented properly, this method enables a complete reproduction of results. As described in section 2.5.4, we provide a high-level object-oriented method for classes 370 of full kinematic experiments, combining Noddy input and output, automatic computation when required, and the additional integration of further methods from external Python packages.
In the following example, we show how we use the pynoddy.experiment methods to investigate error propagation with a Monte Carlo experiment for a complex geological model of the Gippsland Basin. The tectonic history input to Noddy is shown in Fig. 6a. This simplified, but representative ge-375 ological history has been primarily derived from Rahmanian et al. (1990), Norvik and Smith (2001), Moore and Wong (2002) and Lindsay et al. (2012). Each event shown in Fig. 6a corresponds to an event interpreted from the Gippsland Basin, a Mesozoic to Cenozoic oil and gas field in southeastern Australia (Cook, 2006;Rahmanian et al., 1990). Our model basement is Ordovician rocks and the cover sequences include the Oligocene Seaspray and Pliocene Angler sequences. Of particular inter-380 est for oil and gas prospectivity is the Paleocene to Late Miocene Latrobe Group, which includes the Cobia, Golden Beach and Emperor Subgroups (Bernecker et al., 2001 km wavelength) folding is observed, however the basin retains an overall layer-cake stratigraphy. We now want to evaluate how uncertainties in the kinematic parameters of the different tectonic events (Fig. 6a) propagate to the final constructed model (Fig. 6b). The general procedure is briefly outlined here, for more details please see the IPython notebook with the complete example and more thorough descriptions (See documentation and tutorial, Sec. B). 390 We use here the class UncertaintyAnalysis, which contains methods for Monte Carlo-type error propagation and subsequent uncertainty analyses. As a first step, we consider relevant kinematic modelling parameters now as random variables, instead of deterministic variables. The properties of these random variables can be described as probability distributions in several ways. We use here a simple definition in a table, stored in a comma separated file, that can be loaded directly into the 395 object.
We assign normal distributions to location points and layer thicknesses, with a mean value accord-

t a i n t y A n a l y s i s ( h i s t o r y _ f i l e , p a r a m s )
We can now directly generate n random samples from this model with: ua . e s t i m a t e _ u n c e r t a i n t y ( n ) The set of results is, by default, saved directly within the object, and can be extracted in the form of 410 Python numpy arrays for further processing. In addition, A set of standard post-processing methods and utility functions is already implemented in the class definition. For example, it is directly possible to generate analyses and visualisations for the probability of outcomes for a specific geological lithology per voxel (Wellmann et al., 2010;Lindsay et al., 2012), and for the analysis of voxel-based information entropy measures (Wellmann and Regenauer-Lieb, 2012;Wellmann, 2013).

415
In this example of the Gippsland Basin, we perform Monte Carlo error propagation for a set of 32 parameters of all kinematic events in the model, and generate 100 random realisations of the model (see tutorial notebook in documentation). For post-processing, we analyse and visualise results in a 3-D plot of cell information entropies (Fig. 6c). The estimated uncertain areas in the model are clearly visible, and highest uncertainties exist in areas where the effect of uncertainties in different 420 events overlaps (see Fig. 6c).
The previous experiment is a typical example of Monte Carlo sampling methods (Metropolis and Ulam, 1949). One characteristic of the sampling is that all realisations are drawn independently.
Therefore, a parallel implementation of the sampling is directly possible. As one possibility, we pro-vide a parallel sampling scheme implemented in the pynoddy.experiment.monte_carlo.MonteCarlo class, based on the Python threading module, and we used this scheme successfully on a supercomputer. For more information on this possibility, see documentation (Appendix B) and the source code of the monte_carlo.py module.

Discussion
We have presented a newly developed python module for performing scientific experiments with forward. In essence, the combination of input and output generation with on-demand computation allows a high flexibility, as well as an integration of essential aspects of entire kinematic experiments in a single object. As the random state is stored, this encapsulation facilitates easy reproduction of entire scientific experiments with kinematic models.
Limitations of modelling approach presented here are related to (a) the specific kinematic func-450 tions available in Noddy, (b) the conceptual simplification of representing complex dynamical evolutions with purely kinematic functions in general (see Fig. 1), and (c) the limited consideration of surface contact information and measurements.
The first limitation is based on the fact that we did not extend the basic functionality of the kinematic equations implemented in Noddy, and these equations may not be complex enough for specific 455 modelling requirements. For example, the fold model in Noddy is based on a simple fold concept, and this may be a limitation when other fold mechanisms need to be modelled. For a full definition of possibilities and limitations, please see Jessell (1981) and Jessell and Valenta (1996). However, as presented in this manuscript (Fig. 6), in the examples of the Atlas of Structural Geophysics (Sec. A6), or even in recent publications (Armit et al., 2012), complex models can easily evolve from the interaction of multiple kinematic events.
The second limitation, that is the use of kinematic equations instead of a full dynamic simulation, is a significant conceptual simplification and has to be kept in mind when constructing and interpreting results of kinematic modelling, to ensure that the methods are used in the scope where they are valid.
With the examples presented in this manuscript, we wanted to highlight such applications; 465 in addition to the instructive aspect of using kinematic models to teach and visualise the effect of interacting deformational and magmatic events, we believe that main advantages come from the potential to automatically generate multiple model realisations. These methods are facilitated by the fact that the generation of a single kinematic model is typically very fast (in the order of seconds to minutes on a single core) compared to full dynamic simulations. This possibility therefore enables 470 investigation of interaction between simplified deformational events, but with the consideration of uncertainties in event parameters, orders, and types.
The final limitation is that kinematic modelling only allows indirect consideration of actual observations and measurements in the models. An encouraging avenue of investigation is the inclusion of observations facilitated by combining kinematic modelling with interpolation methods (Fig. 1a). 475 We note at this point the similarity between the kinematic modelling methods described in our work, and object modelling methods in geostatistics (Pyrcz and Deutsch, 2014), which are widely and successfully use in reservoir modelling. We envisage that experience from applications of these object modelling methods can be transferred to kinematic modelling concepts based on the flexible methods presented in this work.

480
The methods we have implemented are platform independent, as they are completely implemented in Python, and Noddy itself in C. It is therefore possible to port developed experiments and code easily to other computational environments. We have, for example, tested numerical experiments on supercomputers, a possibility that is especially important for the generation of multiple (i.e. thousands or more) high-resolution model realisations, or the combination with complex post-processing 485 methods. In addition, the platform independence circumvents a limitation of the current GUI for Noddy which is restricted to one operating system. One of the main motivations for the original development of Noddy, for use a teaching tool, is therefore also ensured.
Geological modelling is most often not an end in itself, but the input to further modelling and simulation methods. For example, structural geological models are often used as an input for sub- Future extensions of the developed code will include an optimised application in parallel environments, including a better storage of results (e.g. in HDF5 formats), and a better link to geological 500 data sets and parameters (e.g. through the use of GeoSciML, see Sen and Duffy, 2005;Simons et al., 2006). In addition, we are actively working on developments of additional experiment classes, for example for detailed topological analyses of structural models, and further post-processing and uncertainty quantification methods. Another path of future research is to investigate the possibility to integrate kinematic modelling with Noddy into inference frameworks, for example to test the pos-505 sible inversion of kinematic parameters from observations and geophysical measurements. We hope to include functionality developed by other external users into the main package, and encourage an active participation with successfully developed extensions.

Code Availability
The information provided here is relfecting the current state of the repository at the time of manuscript 510 preparation. In case you find information outdated, please contact the corresponding author.
-For detailed information on the license, see the agreement in the LICENSE file of the repository.

515
-Documentation is available as part of the package and online:

A1 Notes on installation
A successful installation of pynoddy requires two steps: 1. An installation of the python modules in the package pynoddy; 2. The existance of an executable Noddy(.exe) program.
Currently, pynoddy and Noddy can be installed in two alternative ways: (a) directly from the 530 source code with the full repository, or (b) with a direct installation from the Python Package Index and pre-compiled executables. We suggest to use option (a) for the most recent and most complete version of the code. Version (b) is suggested for less experienced users who would like to quickly test and apply kinematic modelling methods. We describe the installation the alternatives in the following.

535
Note: for clarity, we denote command line prompts with a > symbol below: > command to be executed Which should produce the general output: Arguments < h i s t o r y f i l e > < o u t p u t f i l e > < calc_mode > : Note: if the executable is correctly placed in a folder which is recognised by the (Environment) path variable, then you should be able to run Noddy from any directory. If this is not the case, please see Sec. A3.3.

A4.2 Testing pynoddy
The pynoddy package contains a set of tests which can be executed in the standard Python testing 630 environment. If you cloned or downloaded the repository, then these tests can directly be performed through the setup script: The original graphical user interface for Noddy and the compiled executable program for Windows can be obtained from http://tinyurl.com/noddy-site. This site also contains the source code, as well as extensive documentation and tutorial material concerning the original implementation of the software, as well as more technical details on the modelling method itself.

645
The Atlas of Structural Geophysics contains a collection of structural models, together with their expression as geophysical potential fields (gravity and magnetics), with a focus on guiding the interpretation of observed features in potential-field maps.
The atlas is currently available on: http://tectonique.net/asg. The structural models are created with Noddy and the history files can be downloaded from the atlas. Together with the Python package The most convenient way to get started with pynoddy is to experiment with the interactive IPython notebooks, for example to reproduce and adapt the examples given in this manuscript. These notebooks are a part of the repository. The only requirement is to have a running Jupyter installation, 660 see http://jupyter.org for more information. We furthermore plan to have these interactive notebooks available for web-based experiments with pynoddy in the future.
Appendix C: Additional information on Models and Results in this Publication

C1 Example models in this manuscript
All example models presented in this manuscript, respectively the python and pynoddy code to gen-665 erate them, are available as part of the repository. The experiments are directly accessible as Jupyter notebooks to re-generate the presented experiments, or to test different parameters (in docs/notebooks).

C2 Gippsland Basin uncertainty study
The Gippsland Basin model was inspired by previous work of the authors in this region (Lindsay 670 et al., 2012), and further references to the geological setting can be found there. For the purpose of this work, the kinematic parameters for the geological events, as well as the probability distributions consideration of these parameters as random variables, are given in