Keywords: Databases 2, Results relevance 3, Results 4, Generalization

Table of contents

I) Introduction 1. Central nervous system diseases and glia

Central nervous system (CNS) diseases are neurological disorders affecting the function or the structure of the brain and the spinal cord. There is a great number of different types of CNS diseases, spanning from addiction to neurodegenerative disorders. The treatments for those diseases often consist in medication. However, drugs affecting the brain usually come with important and frequent adverse effects or lack of response, varying between individuals [START_REF] Health | Pharmacogenomic Testing for Psychotropic Medication Selection: A Systematic Review of the Assurex GeneSight Psychotropic Test[END_REF].

A combination of a psychotropic and a secondary drug meant to either remove the adverse effects or enhance the psychotropics' efficacy can be considered as a potential solution for this problem. The choice of the secondary drug requires to find the target that is responsible for the psychotropic's limitations. Over recent years, some psychotropic have been determined to have secondary effects on the glia [START_REF] Liu | The Psychostimulant Modafinil Enhances Gap Junctional Communication in Cortical Astrocytes[END_REF], potentially outlining a new type of combination between a neuronal target drug and a glial target drug. Glial cells are the neuron's environmental cells. They constitute around 50% of the brain volume. They are mainly divided in 3 different types: microglia, astrocytes and oligodendrocytes [START_REF] Jäkel | Glial Cells and Their Function in the Adult Brain: A Journey through the History of Their Ablation[END_REF]. While they were originally thought to only support neurons passively, the recent decades have proven that assumption to be false. The glial cells have been found to control multiple aspects of the development and wiring of the brain, as well as providing the metabolic support to neurons and regulating cell-cell interactions [START_REF] Zuchero | Glia in Mammalian Development and Disease[END_REF]. Despite these many discoveries, many of the mechanisms behind their function have yet to be discovered.

Drug Repurposing/Repositioning

The development of a new drug is an expensive, complex and lengthy process costing up to billions of dollars and taking up to decades for a single drug [START_REF] Morgan | The Cost of Drug Development: A Systematic Review[END_REF]. Furthermore, developing a drug is always a risky venture as the success rate for the FDA (USA) evaluations is only around 10% [START_REF] Weng | Pharmacogenetics and pharmacogenomics: a bridge to individualized cancer therapy[END_REF] and once the drug is fully developed and marketable, there is still a risk of commercial failure.

Drug development, while it has the potential to be highly profitable, can fail on many steps of the way and is never guaranteed success [START_REF] Cha | Drug Repurposing from the Perspective of Pharmaceutical Companies[END_REF]. In an effort to reduce the risks and the costs, drug repurposing (also called drug re-tasking, drug repositioning or drug reprofiling) offers an alternative to the classic approach of drug development [START_REF] Pushpakom | Drug Repurposing: Progress, Challenges and Recommendations[END_REF]. Drug repurposing is the process in which already registered and studied drugs will be addressed to another indication. This method allows drug companies to bypass most of development costs as well as having more confidence for safety studies as those have already been conducted for the first commercialization of the compound [START_REF] Morgan | The Cost of Drug Development: A Systematic Review[END_REF]. This process started by empirical discoveries of drugs effects on pathologies they were not supposed to affect. It has now evolved into the use of bioinformatics tools to determine new potential substances usages [START_REF] Talevi | Drug repositioning: current approaches and their implications in the precision medicine era[END_REF].

Those bioinformatics tools use computational approaches to predict novel uses by exploiting various data such as gene expression profiles, protein targets, pharmaceutical mechanisms, etc [START_REF] Ozdemir | Methods for Discovering and Targeting Druggable Protein-Protein Interfaces and Their Application to Repurposing[END_REF] [START_REF] Chou | Theoretical Basis, Experimental Design, and Computerized Simulation of Synergism and Antagonism in Drug Combination Studies[END_REF]. The information is exploited by creating networks of similarities between drugs thus determining which drugs could potentially find new uses [START_REF] Li | Building the Drug-GO Function Network to Screen Significant Candidate Drugs for Myasthenia Gravis[END_REF]. However, those tools are not foolproof. Most of the data generated will be predicted but not proven, therefore the results always have to be considered with caution and be considered only as potential ways to explore via future experiments.

Drugs combination

The combination of multiple drugs is an old practice that can be traced back to millennia ago.

Those drugs were developed to target a single disease at first but it has recently begun to be used to treat multiple diseases at once [START_REF] Foucquier | Analysis of drug combinations: current methodological landscape[END_REF].

The many discoveries in biology these past decades allow for a better understanding of diseases and of the drug used as treatment, allowing science to progress from empiric discoveries to theoretical possibilities. It is now possible to determine if a disease can be better treated by targeting one or multiple targets and the drugs secondary effects can be modulated as well [START_REF] Huang | Systematic prediction of drug combinations based on clinical side-effects[END_REF]. This opens the opportunity for drugs combination which offers multiple benefits compared to simple drugs such as:

-Increased efficacy: one drug enhancing another one by blocking the elimination process.

-Increased efficacy: one drug enhancing the properties of the other without modulating its PK parameters. The only interaction would be pharmaco-dynamic (additional effect).

-Decreased toxicity: one drug counter-acting the negative secondary effects of another drug.

Rationale and Objective

The aim of this internship is to develop a tool to predict potential drug combinations for brain pathologies.

An important proportion of the drugs used to treat central nervous system diseases target neurotransmission while also impacting the glia. We make the hypothesis that this secondary effect is be detrimental but can be modulated by a secondary drug that could be combined with the psychotropic. This tool would look for a selected neuronal drug's secondary effects on the glia by exploring the multiple databases available online and extracting the recorded drug interactions. From those interactions their pathways' contents are selected and filtered to obtain as output the name of the potential glial targets that could be modulated through a second drug.

The first part of the internship consisted in finding databases containing potentially relevant data to our objective to then gather data from the different sources chosen in RDF format. For some of those databases, the data are freely accessible online and need no further work but other databases require to download their data and re-create the database locally, sometimes needing to be converted to the appropriate RDF format as well. The second part consisted in exploring the gathered databases with SPARQL queries and creating an outline of all the information that needs to be extracted to make an appropriate selection. This selection will then remove irrelevant data to our study to retain only glia results.

We focused on Modafinil, a drug used to treat excessive daytime sleepiness in narcoleptic patients, as a proof of concept. An interaction between Modafinil and two glial targets, the connexins 30 and 43, has been discovered in the past years [START_REF] Duchêne | Impact of Astroglial Connexins on Modafinil Pharmacological Properties[END_REF] but is not indicated in drugs or targets databases. As such, finding the connexin in the results, while not a foolproof validation, would serve to add credit to this method.

II) Materials and methods

RDF/SPARQL

Resource Description Framework (RDF) is a model of graph designed to store data and metadata. It was developed by the World Wide Web Consortium (W3C) and is the Semantic Web's main language. It is not formally considered a format but a data model. There are multiple RDF formats such as RDF/XML, Turtle, JSON-LD, OWL, etc.

RDF is based on the concept of a triple: subject, predicate and object. The predicate declares a relationship between the subject and the object thus creating a statement. A collection of these statements represents a directed multi-graph. The data stored under the RDF model can be accessed in two ways: locally or remotely. The remote access is done through endpoints, these endpoints are links to access the data stored within.

SPARQL is a query language used to retrieve and manipulate data stored in an RDF format.

The name is a recursive acronym: SPARQL Protocol And RDF Query Language. It is one of the major technologies of the Semantic Web which is constituted of hundreds of SPARQL services offering access to huge amounts of data, growing continuously [START_REF] Faye | A survey of RDF storage approaches[END_REF].

Virtuoso is a query service allowing for SPARQL queries on multiple remote data endpoints at once. This type of query going through multiple databases is called a federated query.

However, for those federated queries to work, it is necessary for the multiple endpoints used to possess common identifiers for compounds as to be able to connect them.

Unfortunately, this is not the case for all the databases that were used for this internship. The solution consisted in using Python to convert queries' results to make them useable by other databases with which there is no common identifiers.

Two Python packages were used to execute the SPARQL queries:

-RDFLib [START_REF] Boettiger | rdflib: A high level wrapper around the redland package for common rdf applications (Version 0[END_REF] is a library containing an RDF/XML parser and is used on local RDF files that could not be stored in endpoints.

-SPARQLWrapper is a wrapper around a SPARQL service. It is used to send queries on endpoints directly from Python and return the results in a easily useable format.

AskOmics [START_REF] Bettembourg | AskOmics : Intégration et interrogation de réseaux de régulation génomique et post-génomique. IN OVIVE[END_REF] is a visual SPARQL query interface supporting both intuitive data integration and querying while shielding the user from most of the technical difficulties underlying RDF and SPARQL. In this project, it has been be used to convert databases that are not available in RDF but that are available in tsv/csv format. This software will convert files in such format into a RDF format.

AskOmics also offer the possibility to load RDF data into a local private endpoint.

Data used

We looked mainly for two types of databases:

-Drug databases from which we could extract the drugs recorded interactions as well as the sources documenting them.

-Target databases containing proteins, genes and pathways as to link the drugs' databases interactions and connect them to the genome and proteome. On this figure is the example of a classic RDF structure, in this case the Uniprot endpoint. From the Uniprot accession number, the entire data associated is accessible through the corresponding relations. If the accession number is not known, it is possible to obtain it if enough of the associated data is known as the triples in the RDF model work both ways.

SPARQL data

For example, the entry P17302 is the entry for the connexin 43. You can obtain the mnemonic name with the following relation: VALUES (?entry){up_core:P17302} ?entry up_core:mnemonic ?mnemonic

Subject Predicate Object

If the mnemonic name is known but no the entry, this relation can be used as well.

VALUES (?mnemonic){up_core:CXA1_HUMAN}

?entry up_core:mnemonic ?mnemonic This is however one of the simpler structures as, like biopax3 which is used by the Ensembl endpoint. some can prove more complicated III) Results

Databases a. Drug databases

The input of the pipeline being a drug, the first step is to find drug databases. The IDs from other databases as well as the known interactions (and their sources) will be extracted from those. The interactions are mostly genes and proteins interactions. KEGG is one of the most complete databases available online and could have brought plenty of data on both drugs and targets. Unfortunately, it does not offer programmatic access and is not available to download for free. Therefore, it was taken out of the list of drugs databases.

PubChem gathers an important part of all other drugs databases as well as providing bioassays which could have brought new possibilities for target identification. However, the API offers only partial access to the database. It is available for download in RDF, and it is possible to recreate the database locally but it required a 500 Go SSD with 64 Gb of memory. The possibility of using the GenOuest platform to upload the database in the cloud was considered, however due to the short duration of the internship, this would have been a complicated and timeconsuming process for only a few weeks of use.

It should be noted that this choice of database is not exhaustive, the main point of focus for this selection, after their availability, are the drug's interactions. The objective in gathering the maximum of databases is to have a maximum of known interactions, as those differ from database to database. Since there is no centralized database that gathers all the information from all the available sources, we are left with the only option of exploring as many databases as possible. With more time, more databases would have been selected and integrated.

b. Target databases

Next are the target databases, the objective there is to gather enough data on the most commons targets, being genes and proteins.

Wikipathway

Endpoint freely accessible

PharmGKB and DGIdb were available to download freely, not in RDF but instead in TSV. Some small modifications were necessary, mostly in the columns' names but also in some of the data which contained unreadable symbols, to make the files readable by AskOmics. The software then converted the two databases in RDF. These were then accessible through the local endpoint created by the software.

DrugBank was downloadable in one XML file of 1.7 Gb. I developed a Python script using the minidom package, which contains an XML parser, was used to convert in an RDF format:

Turtle. The parser will scan the entire file and extract a DOM tree on which it will be possible to select elements by name. The major inconvenient with this parsing method is that the file is being loaded up entirely and re-created in DOM format. While this is not a problem with small files, the DrugBank database weights 1.7 Gb which can be a problem depending on the computer's power.

XML RDF (Turtle)

From there, data of interest can be extracted with the function getElementsByTagName().

However, this function is recursive, meaning it will search for the given tag name in every subsection of the tree. The parser cannot make the difference between a drug tag referencing a Drug entity and a Drug-Drug interaction. The workaround here is to use the attributes provided in the tags. For example, a Drug entity tag will add information such as: <drug type="biotech" created="01/01/2000"> Not all of the data provided was kept for the conversion as some of it was not necessary and would have slowed the conversion and led to heavier files which would have slowed the

SPARQL queries

The following data was kept:

Identification: Name / DrugBank ID / CAS ID / ATC ID / External IDs Interactions: Targets / Drug-Drugs interactions / Drugs-Proteins interactions / GO terms /

PubMed References

The database was converted in two ways: one file for the entire database and one file per drug (15 000 files in total). The objective was to find which solution was optimal. The conversion took the same amount of time for both methods. Converting into one file amounted to a smaller size since the prefixes at the beginning of each turtle file was only necessary once but the difference was only of 6 Mb (around 3%).

As for queries, they are much slower when looking for data on a single drug on the whole database rather than on one drug file (around 10 minutes against 6 seconds). This means the multiple files method is much more practical and was thus selected. However, the drawback is that since drugs are within separate files, it is not possible to execute a query on multiple drugs at once, instead requiring multiple queries whereas only one would have been necessary with the whole database in one file.

Harmonizome has an API that is programmable through Python, but the documentation is outdated making the API unusable. However, on each of the database's entity webpage, it is possible to retrieve the full page's data in JSON format. This was used by creating a Python script with the urllib package to query the URI leading to this JSON formatted data which could then be extracted.

At this point, all databases are accessible and connected. The SPARQL queries are carried out through the Python package SPARQLWrapper (with the exception of DrugBank for which the RDFlib package is used since the files are not located in an endpoint). Python is used here for multiple reasons:

Remote Endpoint

AskOmics endpoint Local RDF files Webpages

-It allows to execute multiple SPARQL queries at once -Connecting some of the databases can prove tedious due to differences in the way data is stored. Python allows us to convert results from a query to make them readable to another database. The vast amount of packages also ensures that there is a way to use the data in whatever format it is.

-The access to Harmonizome was done through the designed Python API at first and some elements of this API have been recycled to be used with the URL parser.

Throughout the internship, there was no noticeable difference between the query time of SPARQLWrapper on Virtuoso and Virtuoso directly.

Uniprot, Ensembl, Wikidata DGIdb, PharmGKB Harmonizome DrugBank Here we can observe that our set of databases is well connected, with the exception of the two smaller databases DGIdb and Harmonizome. The pipeline input is the name of the drug, in this case the Modafinil. With this, a query is made on Wikidata to extract the DrugBank ID. This query is necessary as the DrugBank database is the only one of those selected from which data cannot be extracted merely with the drug name as the files are named with the drug's DrugBank ID.

Connections between Databases

Pipeline with Modafinil

Queries are then sent to the drugs databases, the Drugs-Target interactions and the sources documenting those are extracted. Three kind of interactions are obtained: proteins, genes and genes variants, this last category is not used for the rest of the pipeline as there is no gene variants database in our list.

WikiPathway is then used to extract all the pathways which contain the genes and proteins interacting with the drug. The content of those pathways is then stored.

However, all the results returned by WikiPathways are not relevant since we are looking for glial results. In an effort to make the results more relevant and to reduce their number, some selections will be made:

-The pathways that do not take place in the brain are removed. A threshold is set where a given percentage of the pathway's proteins have to be present in the brain. This percentage is not 100% because the databases are not exhaustive and the ontology for some of the proteins might not be complete leading to some pathways being relevant but removed because of a few elements.

-The pathways are given a score. This score is computed on the size of the pathway and the amount of recorded interactions occurring in this pathway.

-The proteins are given a score. This score is computed on the amount of times the protein appears in the results and on the score given to the pathway from where it's coming from.

At the time of writing, the computation is still in progress and will not be displayed in the results.

Output

Out of all the drugs databases, a total of 30 interactions have been extracted. They are constituted of 16 genes, 11 proteins and 3 variants. As it has been said before, there is no variant database in our list at the moment so those interactions are not used.

After searching in the pathways containing the interactions, we extracted a total of 9358 human proteins. A filtration is made with Uniprot to select only the proteins localized in the brain. This filtration is not made on the glia yet because Uniprot does not have that level of precision, the purpose is to reduce the amount of results and gin time for the next step. This selection is not done with one of the databases from our list but with a list of glial genes coming from an external source. After this selection, we are left with 1016 results. In these results, we do find the connexins 30 and 43 that were found to interact with Modafinil. [START_REF] Duchêne | Impact of Astroglial Connexins on Modafinil Pharmacological Properties[END_REF]

IV) Discussion

Databases

There is plenty of information available on the many online databases publicly available.

However, this information is dispersed. The recorded drug interactions for example, differ greatly from one database to the other even if some elements are constant.

Few of those many databases were directly available through RDF, whether it is because they are not stored this way or because they chose not to give full access to their data. Some of the most important databases (Uniprot and Ensembl for example) provide a freely accessible RDF endpoint via Virtuoso which make them very easy to access and eventually connect to other databases. It should be noted that those two endpoints are already linked, with the documentation in both of those indicating how to do so. Unfortunately, there are no major drugs databases offering that same service at the exception of ChEMBL which is included in Ensembl. In spite of this, we have established that it is possible to integrate them in a RDF network of databases but it is a precarious answer. Since the data is downloaded, it is not kept to date with the live database and gets no update, potentially meaning that an error in the data could stay until the next download. It also means it is necessary to re-download each time a database is updated, which is a doable but inconvenient process. Furthermore, the files available for download sometimes present errors or discrepancies such as names miswritten or written differently (chemical name instead of usual name). In TSV files, some of the information is sometimes placed in the wrong column.

What is surprising about this state of things is that plenty of databases either used to be available through a public endpoint created by the database administrators or were part of a project to gather them inside a global endpoint. Those endpoints are not activated anymore but can still be found on Bio2RDF, a biological database designed to gather RDF databases.

The availability of the database is not the only blocking step to establish a network of connections. One important issue is the structure of the information. In our case, there is no problem with DGIdb, DrugBank or PharmGKB because there is a direct access to the files and we can control how to set up the data inside the local AskOmics endpoint or in the Turtle files.

Because of this, we have full knowledge on the way the data is contained and can easily create SPARQL queries. However, this is not the case for the public endpoints of Uniprot and Ensembl. Those are both huge databases containing many different types of information which need to be stored and linked together. This results in a complex data network which design' is either specific to the database or following a set of rules. In our case, Uniprot has it's own specific structure while Ensembl is constructed under the biopax3 model. In both cases, the understanding of the structure can be tedious especially when the documentation provided is incomplete, erroneous or imprecise which can be the case for both those databases. For example, some of the prefixes and relationships between entities indicated in the Uniprot documentation either do not exist or are wrongly flagged. As for Ensembl, the documentation for the biopax3 model is very basic, despite it being deep and complex, which means it is necessary to search for a more thorough documentation elsewhere. Fortunately, example queries are provided for the two endpoints which helps as a good introduction for both simple and more complex queries but there is few of them and they do not reflect the vast amount of information available.

Results relevance

There are several factors that need to be taken in account in order to consider our results. Firstly, every result obtained is entirely dependent of the databases' validity. Because we are gathering data from multiple sources to gather as much information as possible instead of cross-examining it, any erroneous result from the database could not be filtered. If an error in a drug data happened to slip by in a drug database, there would be very little way for us to detect it. The solution here lies solely on the choice of the databases and the careful verification of sources from the results. This verification, however, would only remove false positive results but not false negative ones since they would be automatically removed from our results before we could manually detect the error.

Secondly, the design of the pipeline is debatable as well, exploring the pathways is but one of the many possibilities for exploring the available databases. Multiple options were considered and this one was selected because it was considered pertinent and practical. An interaction between a drug and a gene or a protein could logically impact the levels of a molecule from that pathway, but not necessarily. The results obtained are solely predictions and not confirmed interactions, as those could only be obtained via experiments.

Results

We have a total of 1016 proteins as results. This is a number too big to draw conclusions on the relevance of this method, however the two targets that were supposed to be found (connexins 30 and 43) do come up in our list which is an encouraging sign. The pathway filtration has yet to be implemented and may be the decisive step to determine if those results came up by chance or not. If the proteins are ranked by score, the important amount of proteins is not going to be a problem since we will be able to know which ones are truly relevant and which are not.

Generalization

Aside from the functions currently being implemented (pathway and proteins statistical selection as well as glial filtration), there is a lot of room for improvement of the pipeline. For starters, more databases could be used to potentially gather interactions that were not recorded within our own set of databases. A variant database would also be necessary to utilize the variants interactions.

But overall, there is an almost endless amount of possibilities when accessing such large databases. The list of databases chosen, now that they're gathered and linked together, can be used as basis for an entirely different focus just as it could be used to complement the current work. For example, the PubChem database, which wasn't included because of a lack of time, possess Bioassays which could be used to compute similarities between multiple drugs. This way, we could potentially predict common genes and proteins interactions between drugs. Overall, each database added to the list from this point could bring a new potential type of analysis.

IV) Conclusion

There is a huge amount of publicly available information online, whether it is under a RDF format or not. With the exception of a few databases making the deliberate choice of not providing their data, this information is easily accessible, mostly through downloadable files and some databases even provide a free direct programmatic access.

As for the objective of this internship, we obtained a total of 1016 proteins which is too much to consider it achieved at this point. However, the proteins used as validation are present in our list and a ranking score should be implemented in the coming weeks.

Figure 2 .

 2 Figure 2. Glial cells -neurons interactions. [14]

Figure 1 .

 1 Figure 1. Comparison of Drug Development Costs with Revenue Earned After Approval. [4]

Figure 3 .

 3 Figure 3. Structure of the Uniprot SPARQL endpoint. [21]

Figure 4 .

 4 Figure 4. Example of data relationships using the biopax3 model. [22]

Figure 4 .

 4 Figure 4. Example of XML and Turtle files containing the same information.

Figure 5 .

 5 Figure 5. Structure of a DOM tree.

Figure 6 .

 6 Figure 6. The two ways considered to convert DrugBank from XML to Turtle.

Figure 7 .

 7 Figure 7. How the data was extracted

Figure 8 .

 8 Figure 8. Connections between our list of databases.

Figure 9 .

 9 Figure 9. Pipeline's path through selected databases when extracting the potential proteins interactions of Modafinil

Acknowledgements

I thank Olivier DAMERON and Adeline DUCHENE for their help, advice and supervising during this internship. I thank the GenOuest bioinformatics core facility (http://www.genouest.org) for providing access to their cluster


Gautier Laurent
email: laurent@gocad.org

Guillaume Caumon
email: caumon@gocad.org

Antoine Bouziat
email: antoine.bouziat@cggveritas.com

Mark Jessell
email: mark.jessell@get.obs-mip.fr

A parametric method to model 3D displacements around faults with volumetric vector fields $

Keywords: structural modeling, fault, parameterization, kinematics, time integration

This paper presents a 3D parametric fault representation for modeling the displacement field associated with faults in accordance with their geometry. The displacements are modeled in a canonical fault space where the near-field displacement is defined by a small set of parameters consisting of the maximum displacement amplitude and the profiles of attenuation in the surrounding space. The particular geometry and the orientation of the slip of each fault is then taken into account by mapping the actual fault onto its canonical representation. This mapping is obtained with the help of a curvilinear frame aligned both on the fault surface and slip direction.

This formulation helps us to include more geological concepts in quantitative subsurface models during 3D structural modeling tasks. Its applicability is demonstrated in the framework of forward modeling and stochastic sequential fault simulations, and the results of our model are compared to observations of natural objects described in the literature.

Introduction

Faults dramatically impact fluid flow, mineralization, facies localization and the geometry and connectivity of rock units. In subsurface modeling, an accurate description of faults is therefore paramount in maximizing a model's predictive capabilities. In general, uncertainties associated with faults are significant, particularly concerning their connectivity and the displacement of surrounding rocks. These two aspects are intimately related to the kinematics and mechanics of the fault. Unfortunately, the compatibility with deformation history, tectonics, kinematics and mechanical concepts is generally secondary in current processes of quantitative 3D modeling, the main objective being to fit the data while honoring geometrical constraints such as minimizing the curvature of the structures [START_REF] Caumon | Surface-based 3D modeling of geological structures[END_REF].

Typical geomodeling workflows proceed by interpolating the stratigraphic information contained in data points while faults are taken into account by introducing topological discontinuities [START_REF] Mallet | Geomodeling[END_REF][START_REF] Calgagno | Geological modelling from field data and geological knowledge part i. modelling method coupling 3D potential-field interpolation and geological rules[END_REF][START_REF] Caumon | Surface-based 3D modeling of geological structures[END_REF]. While such approaches have demonstrated their efficiency and are in daily use by the natural resources industry, they lack an explicit control on kinematics and mechanics, which can lead to implausible structures [START_REF] Caumon | Current bottlenecks in geomodeling workflows and ways forward[END_REF]. A validation step is then required for rejecting structurally in-compatible models e.g. by simulating the retrodeformation that restores the geological structures to their supposed initial state [START_REF] Dahlstrom | Balanced cross sections[END_REF][START_REF] Tanner | Threedimensional retro-deformation of the lechtal nappe, northern calcareous alps[END_REF][START_REF] Moretti | Kine3D: A new 3D restoration method based on a mixed approach linking geometry and geomechanics[END_REF][START_REF] Maerten | Chronologic modeling of faulted and fractured reservoirs using geomechanically based restoration: Technique and industry applications[END_REF][START_REF] Durand-Riard | Balanced restoration of geological volumes with relaxed meshing constraints[END_REF]. Even if such approaches have proved their ability to highlight certain structural inconsistencies, they imply a repeated trial-and-error process to achieve a completely kinematically compatible model, which would require prohibitive computational time. For these reasons, the introduction of kinematics and geomechanics remains a major bottleneck of geomodeling workflows, albeit crucial for the predictive capabilities of the models [START_REF] Fletcher | Can we understand structural and tectonic processes and their products without appeal to a complete mechanics[END_REF].

Several approaches aim at producing numerical models of structural surfaces that natively honor some geological rules, for example: developability [START_REF] Thibert | A direct method for modeling and unfolding developable surfaces and its application to the ventura basin (california)[END_REF], thread geometry of fault surface [START_REF] Thibaut | An inverse method for determining three-dimensional fault geometry with thread criterion: application to strike-slip and thrust faults (western alps and california)[END_REF], sedimentation and compaction rates consistency [START_REF] Mallet | Space-time mathematical framework for sedimentary geology[END_REF]) and fold models [START_REF] Kaven | Constraining surface interpolations using elastic plate bending solutions with applications to geologic folding[END_REF][START_REF] Hjelle | A hamilton-jacobi framework for modeling folds in structural geology[END_REF]. Because 3D displacement patterns associated with natural faults play a prominent role in fault characterization, we suggest that taking them explicitly into account while building the faults is a key aspect to increasing their consistency.

The application of such concepts to uncertainty modeling and to inverse problem solving [START_REF] Cherpeau | Method for stochastic inverse modeling of fault geometry and connectivity using flow data[END_REF][START_REF] Jessell | Towards an integrated inversion of geoscientific data: What price of geology?[END_REF][START_REF] Georgsen | Fault displacement modelling using 3D vector fields[END_REF] requires us to express a fault's displacement with appropriate parameters:

• The parameters should preferably correspond to the structural characteristics of geological objects to make their interpretation and use easier.

• They have to be limited in number in order to maintain the dimensionality of these problems to acceptable limits.

This paper presents an approach for parameterizing the faults and their associated displacements in the form of a volumetric vector field, in a spirit similar to [START_REF] Jessell | Structural geophysics: integrated structural and geophysical modelling[END_REF] and [START_REF] Georgsen | Fault displacement modelling using 3D vector fields[END_REF]. It consists of an extension of the fault parameterization presented in Cherpeau et al. (2010b[START_REF] Cherpeau | Method for stochastic inverse modeling of fault geometry and connectivity using flow data[END_REF]. The proposed approach is based on previous work in Computer Graphics (Von Funck et al., 2006), which was first adapted to geological modeling in an earlier conference paper [START_REF] Bouziat | Vector field based fault modelling and stratigraphic horizons deformation[END_REF]. Our paper completes this fault parameterization and introduces significant improvements to the methodology with a particular attention to fault kinematics and partitionning between nearfield and far-field deformation.

Our contributions to a complete fault parameterization are:

• The definition of a 3D curvilinear fault space offering a general and appropriate frame for displacement computation (section 3.2).

• The differentiation of the displacement evolution in the three principal directions of a fault: away from fault, along fault in the slip direction and along fault orthogonal to slip direction (section 3.4).

• The ability to combine different flanking structures at different scales around the same fault (Fig. 10).

• The ability to use complex prior information to characterize the displacement field even when few or no data are available.

• The recourse to a time integration scheme to progressively build the displacement (section 3.5), which allows us to combine displacements with different tectonic origins.

We present synthetic applications of this model in a forward modeling context (section 4.1). Its integration into stochastic sequential simulations of fault networks for improving fault data clustering is also depicted (section 4.3). Finally, adaptation to complex fault cases is presented through the example of a roll-over anticline (section 4.2).

Related work

Numerous approaches have allowed fault-related displacements to be taken into account as an emerging character from either the geometrical interpolation of geological structures due to topological discontinuity introduced by faults [START_REF] Caumon | Surface-based 3D modeling of geological structures[END_REF] or the geometrical and mechanical restoration of horizon structures [START_REF] Egan | Computer modelling and visualisation of the structural deformation caused by movement along geological faults[END_REF][START_REF] Moretti | Kine3D: A new 3D restoration method based on a mixed approach linking geometry and geomechanics[END_REF][START_REF] Maerten | Chronologic modeling of faulted and fractured reservoirs using geomechanically based restoration: Technique and industry applications[END_REF][START_REF] Durand-Riard | Balanced restoration of geological volumes with relaxed meshing constraints[END_REF]. For a complete fault parameterization, the associated displacements have to be explicitly integrated in the description of the fault.

The context of forward modeling requires fault operators capable of representing all possible ranges of displacements along and across faults, even when few or no data are available. To overcome this hurdle, [START_REF] Jessell | Structural geophysics: integrated structural and geophysical modelling[END_REF] present several fault operators analytically modeling the 3D displacement fields associated with different kinds of canonical faults. For regional faults going through a studied domain, the displacement fields in fault blocks are modeled as pure translation or rotation. For faults of shorter extent, an elliptical decrease of slip intensity based on [START_REF] Walsh | Distributions of cumulative displacement and seismic slip on a single normal fault surface[END_REF] is introduced.

Similar concepts of fault operator have also been developed for model editing in fault uncertainty models [START_REF] Hollund | Havana-a fault modeling tool[END_REF][START_REF] Georgsen | Fault displacement modelling using 3D vector fields[END_REF], using two kinds of fault models: a piece-wise planar model for large faults and an elliptic model for small-scale faults. The curvature of fault surfaces can be accounted for by requiring constant distance between the displaced points and the fault surface [START_REF] Jessell | Structural geophysics: integrated structural and geophysical modelling[END_REF][START_REF] Georgsen | Fault displacement modelling using 3D vector fields[END_REF]. Because slip distribution is generally far more complex than depicted by classical elliptical models [START_REF] Barnett | Displacement geometry in the volume containing a single normal fault[END_REF] it is possible to derive the slip field by kriging horizon displacements projected on the medium plane of the fault [START_REF] Georgsen | Fault displacement modelling using 3D vector fields[END_REF].

In this paper, we present a general model based on the observation that, in spite of the wide diversity of faults, some common characteristics emerge and make it possible to describe them in a unified model, based on two main components:

• a 3D curvilinear fault frame, whose axes are oriented according to the orientation of the fault surface and the displacement direction (detailed in section 3.2).

• the profiles of the attenuation of the displacement along the three axes.

The fault frame makes it possible to map complex fault geometries onto canonical cases, and the attenuation profiles enable the modeling of complex, theoretical and/or observed, displacement patterns. Displacement attenuation profiles are controlled by a limited set of geometric parameters, making this model compatible with both forward modeling and inverse approaches [START_REF] Cherpeau | Method for stochastic inverse modeling of fault geometry and connectivity using flow data[END_REF][START_REF] Jessell | Towards an integrated inversion of geoscientific data: What price of geology?[END_REF][START_REF] Georgsen | Fault displacement modelling using 3D vector fields[END_REF].

A parametric model describing faultrelated displacements

General presentation of the model of displacements

The model presented in this section aims at deforming the structures cut by a fault in a kinematically consistent way. It relies on the modeling of the displacement field which represents the effects of the fault on the surrounding structures. Two kinds of displacement fields associated with faults are generally considered [START_REF] Barnett | Displacement geometry in the volume containing a single normal fault[END_REF]: the farfield and the near-field, representing respectively the continuous and discontinuous parts of the displacement (Fig. 1).

The far-field describes the global displacement field in which a fault occurs. At the large scale (with regard to the size of a fault) only the far-field is perceptible. At a smaller scale, a fault localizes deformation which enables the accommodation of part of the stress related to the far-field displacements. This accommodation comes in the form of an additional displacement located around the fault, referred to as near-field. This kind of displacement affects the surrounding rocks and produces flanking structures, including normal and reverse drag (Fig. 2). These two terms designate ductile deformation of geological markers cut by a fault. They are distinguished by the direction of the curvature of the resulting folds:

• Normal drag describes a decrease of apparent displacement near the fault surface producing bending towards the opposite direction of the block displacement (Fig. 2a).

• Reverse drag is the opposite phenomenon. It corresponds to an increase of the apparent displacement near the fault surface bending the horizon towards the direction of displacement (Fig. 2b).

The term "normal drag" comes from the similarity with the geometry produced by friction phenomena [START_REF] Hamblin | Origin of reverse drag on the downthrown side of normal faults[END_REF][START_REF] Billings | Structural geology[END_REF] and is now well established even if it has been recognized to be misleading [START_REF] Hamblin | Origin of reverse drag on the downthrown side of normal faults[END_REF][START_REF] Hobbs | An outline of structural geology[END_REF][START_REF] Grasemann | Reverse and normal drag along a fault[END_REF]. Indeed, the frictional resistance is unable to properly explain normal drag and would limit the development of reverse drag [START_REF] Reches | Drag along faults[END_REF], which seems contradictory to the fact that normal and reverse drag can be observed together on a single fault.

c. Resulting displacement b. Near-field displacement a. Far-field displacement Figure 1: The different displacement fields related to a fault. a. Far-field: pure shear displacement field with diagonal principal direction (white arrows: shortening, black arrows: elongation). Faults preferentially appear with an angle below 45 • from principal stress direction, depending on the friction coefficient of the host rock (dashed lines: two possible orientations for a fault). In this context, a geological marker (double line) would rotate with an angle depending on the orientation relative to the principal directions. b. Near-field: the development of a fault generates an additional displacement decreasing up to the near-field boundary (dashed line). The attenuation of the near-field displacement causes any orthogonal horizon to form a reverse drag. c. Resulting displacements: the combination of far-field and near-field displacements results in a complex displacement field. The isolines show the distribution of the displacement magnitude (orange is high, blue is low).

The interpretation of fault drag is complicated their diverse origins. The curvature of the fault surface is one of the first general effects accounting for flanking structures, e.g. the roll-over anticlines observed in the hanging-wall of listric faults [START_REF] Hamblin | Origin of reverse drag on the downthrown side of normal faults[END_REF]. The attenuation of the near-field displacement around faults of limited extent also plays a role in producing flanking structures as it naturally produces reverse drag [START_REF] Barnett | Displacement geometry in the volume containing a single normal fault[END_REF], even for planar faults. Normal drag is also likely to appear due to a low angle intersection between marker and fault [START_REF] Grasemann | Sense and non-sense of shear in flanking structures[END_REF][START_REF] Grasemann | Reverse and normal drag along a fault[END_REF] or due to ductile deformation occurring before or in relation to fault rupture [START_REF] Reches | Drag along faults[END_REF].

It is important to notice that normal and reverse drag generally act at different scales, normal drag being localized closer to the faults' surface or tips. The ability to superimpose these two effects with appropriate scales, as demonstrated in Fig. 10, is necessary for an accurate representation of flanking structures.

As in most of the literature, we focus on the nearfield displacement for two reasons:

• The modeling of the far-field is demanding because it cannot be directly deduced from the fault geometry. It depends on the tectonic context and is beyond the scope of this study.

• The near-field is generally predominant in the characterization of flanking structures (Fig. 2).

The observation of faults in the field reveals that actual near-field displacement patterns may be quite complex, for example in the case of faults growing by segment linkage [START_REF] Walsh | Formation of segmented normal faults: a 3-d perspective[END_REF][START_REF] Kim | The relationship between displacement and length of faults[END_REF]. Rather than directly taking all possible displacement patterns into account, the model presented here focuses on the description of simple fault patterns, i.e. faults having the following characteristics:

• The displacement intensity is globally decreasing from a point of maximum displacement, referred to as fault center.

• The fault surface and the displacement intensity isolines along the fault surface are roughly circular or elliptical, i.e. convex 2D shapes with major and minor geometrical axes.

• The displacement direction is parallel to one of these axes.

It is then possible to model more complex displacements by combining several simple fault patterns.

In the present model, finite and infinite faults will be considered with the following assumptions:

• For finite faults, the continuity of displacements of the matter beyond the faults' tip line implies the continuity of the far-field displacements. The discontinuous part of the final displacement field is represented by the near-field a. Normal drag b. Reverse drag displacements. This near-field implies a volumetric strain which can be locally non-zero, but we assume that its values globally compensate inside the envelope of the near-field ensuring the global conservation of rock volume. This condition is however relaxed for faults whose near-field envelope intersects a free surface, such as the topographic surface.

• Very large faults can be considered as infinite with regard to the domain of study. For infinite faults, near-and far-fields are entangled and are considered as a single discontinuous far-field. Some local deformation leading to normal drag of limited extent is also likely to be observed. They are interpreted as residuals of pre-or peri-faulting ductile deformation [START_REF] Reches | Drag along faults[END_REF]. In our approach, very large faults are represented by a discontinuous far-field with a possible special near-field representing ductile deformation (section 4.2). No assumption on local or global volumes is made because their conservation is not systematic [START_REF] Means | Stretching faults[END_REF]. These conditions are replaced by assumptions on the extension rate in a given direction.

The computation of the near-field displacement relies on a given fault geometry, a set of structural parameters representing the fault's extent and the profile of attenuation of the displacement in the surrounding rocks. The final displacement is obtained by progressively integrating an instantaneous displacement field v over time (section 3.5), which al-lows to account for the variation of the displacement direction and intensity. Furthermore, the calculation of v is rendered independent of the fault geometry by using a curvilinear fault frame presented in the next section.

3D curvilinear fault space

Despite a large variety of fault types, two principal directions can generally be identified: the direction orthogonal to the fault surface and the direction of slip. Our model implements these observations through the use of a 3D curvilinear frame F, inspired by the one describing chronostratigraphic coordinates [START_REF] Moyen | 3D-parameterization of the 3D geological spacethe GeoChron model[END_REF][START_REF] Mallet | Space-time mathematical framework for sedimentary geology[END_REF]. Each of its coordinates is defined by three scalar fields in 3D space:

• g 0 : represents both the geometry of the fault surface (iso-surface of value 0) and the signed distance to its surface.

• g 1 : is defined by the local normal to the fault and the slip direction. It complements g 0 for representing the displacement direction: along the fault surface, the isolines of g 1 correspond to the fault striae.

• g 2 : is a complementary coordinate orthogonal to both g 0 and g 1 . g 2 is used to describe the lateral evolution of the displacement field on the fault surface.

These three coordinates define an implicit mapping between the actual geometry G and the regular fault space F where the fault appears as an ideal elliptical fault (Fig. 3). To complete this mapping, a fault center C is also required. C is chosen as the point of maximum displacement along the fault surface, which facilitates the computation of the displacement profiles (section 3.4).

The distance D to the fault center is defined in the frame F as the norm of the vector of coordinates (g0, g 1 , g 2) denoted as g:

∀g ∈ F, D(g 0 , g 1 , g 2) = g 2 0 + g 2 1 + g 2 2 (1)
For convenience, g is normalized so that D becomes 1 at the boundary of the near-field influence zone (see Fig. 1b).

The fault surface is defined by the positions where: g 0 = 0 and D g ≤ 1. The tip line of the fault is defined by: g 0 = 0 and D g = 1. An algorithm to obtain adequate coordinates while modeling a fault with a complex geometry is given in section 4.1.

Constraints on displacement direction

The direction and polarity of v are directly derived from the fault frame F. The definition of F implies that v is parallel to the isovalues of g 0 and g 1 . This is ensured by computing its direction vector u from the cross product of the gradient of these coordinates, respectively denoted ∇g 0 and ∇g 1 :

u = v v = k ∇g 0 × ∇g 1 ∇g 0 • ∇g 1 (2)
The factor k in the previous equation denotes the polarity of the displacement. Its sign is the same as the sign of g 0 , which determines in which fault block each point is located. This definition could be problematic for the points laying exactly on the fault surface because their values of g 0 are exactly 0. This is addressed by considering the g 0 value of the neighbors of the problematic points (which implies that the numerical representation of the displaced objects are discontinuous across the fault).

The identification of the fault block is also used when a different behavior is expected on the footwall and the hangingwall of the fault. For example, an uneven repartition of the total relative displacement can be represented by introducing a partition factor γ [START_REF] Georgsen | Fault displacement modelling using 3D vector fields[END_REF] defined as follow:

γ(g 0) = γ h if g 0 ≥ 0 γ f if g 0 < 0 (3) with γ h = 1 -γ f (4)
γ h and γ f are the ratio of the relative displacement affected respectively to the hangingwall and the footwall. One of the fault blocks can also be fixed, setting γ h or γ f to 0, in order to apply the complete displacement on the other fault block.

It is possible to model more complex displacements, as presented by [START_REF] Grasemann | Reverse and normal drag along a fault[END_REF], with this framework, for instance by adding a displacement component orthogonal to g 0 . Also, in certain cases, the displacement direction may not be parallel to the minor or major axis of the fault surface, as required in section 3. Such fault can be handled by creating g 1 and g 2 according to the fault surface axis and by rotating these coordinates around C before the computation of u.

Displacement attenuation profiles

The amplitude of the near-field displacement is supposed to be maximum at a certain point C of the fault surface and it decreases to zero up to a certain distance (Fig. 1b). The profile p(g 0 , g 1 , g 2) of the displacement attenuation evolves differently in each principal direction of F. Three profiles p 0( g 0) , p 1( g 1) and p 2( g 2) are introduced in order to explicitly control the attenuation in the directions g 0 , g 1 and g 2 respectively.

The amplitude of v is computed by multiplying the maximum displacement v max at C by a partition factor γ and an attenuation function ρ: v = v max γ(g 0) ρ(g 0 , g 1 , g 2)

(5)

The attenuation function ρ is computed from the attenuation profiles by applying the following steps:

• The attenuation due to the position along the fault surface is first computed by combining p 1( g 1) and p 2( g 2) so that ρ(0, g 1 , g 2) evolves progressively from p 1( g 1) where g 2 equals 0 to p 2( g 2) where g 1 equals 0. This is obtained by considering the polar coordinates (r, α) in the fault surface, where r denotes the distance between a point (g 1 , g 2 ) and the center C and α the angle between this point and the axis g 1 :

ρ(0, g 1 , g 2) = p 1( r) cos 2 (α) + p 2( r) sin 2 (α) (6) with r = g 2 1 + g 2 2 (7) cos 2 (α) = g 2 1 / g 2 1 + g 2 2 (8) sin 2 (α) = g 2 2 / g 2 1 + g 2 2 (9) Geological Space Regular Fault Space C fault tip non-zero displacement part g 2 g 1 g 0 g 2 g 1 g 0 a.
b.

c. • The effect of the attenuation due to the distance to the fault surface following g 0 is then taken into account by multiplying by p 0( g 0) :

ρ(g 0 , g 1 , g 2) = p 0( g 0) p 1( r) g 2 1 + p 2( r) g 2 2 g 2 1 + g 2 2 (10)
An example of profiles honoring basic geological conditions is depicted in Fig. 4. They can be adapted to specific observations and models of profiles along and across fault [START_REF] Barnett | Displacement geometry in the volume containing a single normal fault[END_REF][START_REF] Walsh | Distributions of cumulative displacement and seismic slip on a single normal fault surface[END_REF][START_REF] Reches | Drag along faults[END_REF][START_REF] Cowie | Fault tip displacement gradients and process zone dimensions[END_REF]. In our implementation, a polynomial spline is used for interactively editing the profiles by modifying the slope at the extremities and introducing control points.

It is possible to obtain more complex profiles by mixing several simple profiles. For example, the profiles presented in [START_REF] Hamblin | Origin of reverse drag on the downthrown side of normal faults[END_REF], where superimposed reverse and normal drag are observed, are modeled by linearly combining two profiles with different fault drag radii. Profiles resulting from segment-linkage can also be implemented by combining several simple profiles with centers placed at different locations.

Time integration of the instantaneous displacement field

As the magnitude and direction of the displacement fields change along the fault, we propose to use a discrete time scheme in order to integrate instantaneous displacement steps instead of directly considering the complete displacement. This approach has numerous interesting properties [START_REF] Theisel | Topological methods for 2d time-dependent vector fields based on stream lines and path lines[END_REF][START_REF] Funck | Vector field based shape deformations[END_REF], with appropriate displacement fields and time steps:

• It avoids non-physical self-intersection of the displaced material.

• It keeps the detailed geometry and continuity of the surrounding structures. The curve on the left presents the attenuation followed by p 1 and p 2 , i.e. along the fault surface. A bell-shaped profile [START_REF] Cowie | Fault tip displacement gradients and process zone dimensions[END_REF] is represented here, but other profiles could be used instead. The curve on the right represents p 0 , with a profile close to the profiles presented in [START_REF] Barnett | Displacement geometry in the volume containing a single normal fault[END_REF]. The value of the derivatives are controlled at the extremities and center of the profiles (arrows) and additional control points can be introduced in order to tune the curves (squares). Note that p 0 is represented as negative for g 0 < 0 (this behavior is managed by the factor k in eq.( 2))

In addition, it is consistent with fault behavior, whereby observed displacements generally result from one or a succession of relatively continuous displacement events.

The time integration process is performed by successively summing up the small displacements computed over small time steps. The duration of the time step τ is determined at the beginning of the process. It is expressed as a percentage of the total displacement duration. It can be refined to better fit the curvature of the fault or to adapt to the displacement intensity variations and to the deformed object discretization.

The fault frame is initialized by computing g 0 , g 1 and g 2 (see section 4.1). For each point of each object to be deformed, the corresponding fault block is computed by affecting the sign of g 0 to the variable k. Then, for each time step τ :

• If the geometry of the fault or its structural parameters have changed, the fault frame is updated.

• For each point X:

-The direction of the displacement vector u is evaluated from eq.( 2).

-The intensity of displacement v is evaluated from eq.( 5).

-The position of the point is updated:

X ← X + τ v u
Another interesting aspect of this approach is that it is possible to accumulate the effect of several tectonic elements, as long as an instantaneous displacement field can be provided. For example, it is possible to accumulate several faults acting simultaneously (Fig. 6). This is also a way to integrate the far-field effect as demonstrated in figure 8 4. Application

Deformation of horizons cut by finite faults

This section presents some examples of the application of our method. A first basic example illustrates the effects of a single fault of limited extent on a volumetric stratigraphic sequence (Fig. 5). Because the fault is planar, the fault frame is obtained by taking the distance to the fault center in the three principal directions of the fault. These directions correspond to the normal of the fault plane for g 0 , the direction of slip g 1 and their orthogonal direction g 2 .

The second example (Fig. 6) presents a more complete application with two non-planar faults simulated by a parametric approach (Cherpeau et al., 2010b[START_REF] Cherpeau | Method for stochastic inverse modeling of fault geometry and connectivity using flow data[END_REF].

In this case, the fault coordinates are computed by a constrained interpolation process [START_REF] Frank | 3D-reconstruction of complex geological interfaces from irregularly distributed and noisy point data[END_REF]Caumon et al., 2012 (in press). g 0 is an output of the parametric simulation of the geometry of the fault. g 1 and g 2 are extrapolated from the fault center with the following constraints:

• The values of g 1 and g 2 are 0 at the fault center C.

• ∇g 1 is orthogonal to the slip direction. • The gradient of the three coordinates are mutually orthogonal, which is ensured by setting:

∇g 1 ⊥ ∇g 0 , ∇g 2 ⊥ ∇g 0 and ∇g 2 ⊥ ∇g 1 .

• The norm of the gradients is set to 1 at C.

The coordinates are finally normalized by dividing g 1 and g 2 by their extreme values observed on the border of the fault and g 0 by its values observed at a specified distance from the fault center. The fault coordinates used for Fig. 6, before and after normalization, are presented in Fig. 7. The near-field displacement alone is generally not sufficient to consistently model the geometry of faulted stratigraphic horizons. Indeed, it only represents additional deformation made possible by the mechanical discontinuity of the fault, but the source of this displacement most often comes from the far-field conditions. Deformations representing only the near-field may seem artificial for this reason and taking the far-field into account is necessary to obtain useful models.

We consider two possibilities to take the far-field into account:

• Data driven approach: which consists of building continuous horizons first, by interpolating the stratigraphic data [START_REF] Caumon | Surface-based 3D modeling of geological structures[END_REF], and then applying the fault operator to add the near-field displacements to the model. With this approach, the far-field is not modeled explicitly but its effect is implicitly accounted for by relying on the data. Therefore, this modeling strategy does not highlight the link between the observable continuous deformation and the actual far-field related to the fault operator. For example, folds appearing before or after the faulting are ignored, possibly leading to inconsistent results. The evolution of the fault can not be represented properly as the integration does not start with an actual geometry of the horizon. Also there is no guaranty that the parameters used to describe the near-field displacements are consistent with the underlying far-field displacements contained in the continuous horizon. Nevertheless, this aproach is a pragmatic modeling choice because it directly honors stratigraphic data.

• Process oriented approach: which starts with an initial horizon geometry, possibly integrating hypothesis on preexisting faults and folds. The far-field displacements are then modeled as an additional vector field, which is superimposed on the near-field displacements during the time integration process (Fig. 8). This approach is more consistent with fault's kinematic than the data driven approach, but, as mentioned earlier, the parameterization of the far-field is delicate, especially concerning the data conditionning. This latter problem should be addressed by an inverse approach as in Allmendinger (1998); [START_REF] Cardozo | Optimized trishear inverse modeling[END_REF].

The second approach makes it possible to take into account the displacement and deformation of the fault surface during the deformation. These displacements change the way the fault affects the surrounding horizons and finally their resulting geometry. Indeed, the evolution of the angle between the fault and the horizons potentially changes the fault drag mode during deformation, from reverse to normal drag or vice versa [START_REF] Grasemann | Sense and non-sense of shear in flanking structures[END_REF].

Large faults and non-planar faults

The treatment of large faults slightly differs from the computation of a short fault displacement field. The profile should be of infinite extension to obtain a constant displacement intensity in each fault block. Only g 0 and g 1 are then needed to control the direction of the displacement. However, with the increase of the displacement intensity, normal drag can become visible [START_REF] Hamblin | Origin of reverse drag on the downthrown side of normal faults[END_REF][START_REF] Reches | Drag along faults[END_REF]. One way to model this consists of introducing an additional profile following g 0 with a limiting effect on the displacement intensity (Fig. 10). This profile is either applied during the deformation, to model possible frictional phenomena, or all at once at the end or at the beginning of the time integration process in order to represent pre-or peri-faulting ductile deformation [START_REF] Reches | Drag along faults[END_REF].

The accentuated curvature of the fault surface produces reverse or normal drag that are not directly explained by the attenuation of the displacement orthogonally to the fault, such as roll-over anticlines observed in the down-thrown block of listric faults [START_REF] Hamblin | Origin of reverse drag on the downthrown side of normal faults[END_REF]. We propose to account for such phenomena by adapting the displacement intensity to the geometry of the fault. The direction of the displacement field is computed as usual with eq.( 2). The intensity of v is locally adapted to ensure that the curvature of the fault is respected and that the displacement path does not diverge

Before normalization

After normalization g0:

-1000 0 1000 g1:

-1000 0 1000 g2:

-1000 0 1000

g0: -1 0 1 g1: -1 0 1 g2: -1 0 1
Figure 7: The fault frame coordinates used for the red fault of Fig. 6. Each picture displays the layers to be deformed, the two faults and the isosurfaces g 0 = 0 and g 1 = 0. On the first row they are respectively painted with g 0 , g 1 and g 2 before their normalization. The second row presents the result of the normalization for each coordinate. The values bigger than 1 or smaller than -1 are painted in red and blue respectively.

from the fault surface during the time integration. A global displacement vector V is associated with each fault block. For example, V can be deduced from regional knowledge about extension rate in a rifting context. v is then integrated until each part of the fault block has been displaced by a distance V in the direction of V :

(x -x 0 ) . V V = V (11)
where x is the current position of a displaced point and x 0 its original position. This process corresponds to vertical or inclined shear hypothesis [START_REF] Yamada | Application of geometric models to inverted listric fault systems in sandbox experiments. paper 1: 2d hanging wall deformation and section restoration[END_REF], depending on the inclination of V .

The algorithm presented in figure 9 is applied to integrate the displacements until this criterion is reached. An example of the application of this algorithm is given in figure 10.

Taking fault displacements into account in stochastic fault network modeling

In the context of geological mapping and natural ressource management, the incomplete observations for building quantitative subsurface geomodels raise significant uncertainties about structural geometries. As a consequence, several interpretations are generally possible. Deterministic approaches, relying on a single geomodel, are then inadequate to properly represent subsurface uncertainties. A better sampling of the possible models can be obtained by stochastic approaches which generate series of equiprobable models honoring available data and generated though prior statistical distributions taken from analog studies.

Fault networks are generally particularly sensitive to uncertainty (number of faults, dimensions, connectivity. . . ) with dramatic impact on rock volume and fluid flow predictions, which justify the use of stochastic approach to fault network modeling [START_REF] Hollund | Havana-a fault modeling tool[END_REF][START_REF] Holden | Stochastic structural modeling[END_REF]250 m Figure 8: Combination of continuous far-field and discontinuous near-field displacements. Stratigraphic horizon are cut and displaced by two normal faults while they are affected by a global vertical shear (black arrows). This far-field is expressed as a 3D vectorial field and integrated over time in the same way as the near-field displacement associated with faults. The near-fields used here are similar to those of Fig. 6. The down-thrown blocks undergo 65 percent of the relative displacement for each fault. Cherpeau et al., 2010a[START_REF] Cherpeau | Method for stochastic inverse modeling of fault geometry and connectivity using flow data[END_REF]. Interpreted seismic sections (Fig. 11a) represent typical data used for fault network modeling. They contain fault sticks, i.e. groups of fault markers localizing the discontinuity introduced by faults in the seismic reflectors. Because of the lack of information between the cross-sections, the number of faults and the way fault surfaces connect the different fault sticks can be highly uncertain.

To address these types of uncertainty, Cherpeau et al. (2010a,b) suggests generating faults sequentially by randomly drawing a set of geometric parameters from conditional probability distributions. These parameters describe the position of the fault center, its strike, dip, sinuosity and extension in strike and dip directions. Each set of parameter values allows us to define a fault surface, whose geometry can be adapted to fit fault sticks located close to the fault surface. Starting from an initially regular surface, the fault is deformed to fit each compatible fault stick. The evaluation of the compatibility of each fault stick is based on its location, on the curvature its integration would produce and on the intensity and orientation of the displacement observed along the stick. This step, called fault stick clustering, is crucial to ensure the likelihood of the output structural model.

Figure 11 presents an application of stochastic fault networks construction based on a simple data set of 2 synthetic cross sections. In this example, two kinds of faults are expected:

• Family 1: a first generation of normal faults with a NW-SE orientation dipping to the NE (yellow and blue sticks).

• Family 2: a second generation of vertical strike slip faults oriented N-S (green sticks).

The data clustering step determines if a given fault trace can be explained by, and thus associated with, the fault being currently modeled. Here, the number of faults to be modeled from family 1 is particularly uncertain. This family could produce a single fault, gathering yellow and blue sticks, or two isolated faults, which affects the connectivity of the fault network.

Among other criteria, the fault trace position is inspected and its consistency with the current geometry of the fault is checked. However, this position may have been assigned by the displacement -get its initial position X 0 -compute its displacement vector v

for each point X:

compute intensity of displacement in the direction of V increment X: X=X+v increment X: X=X+rv d< V no yes increment d b. l = d l > d l < d l = (X + v -X 0 ) • V r= d-l v•V
compute the part of v to apply:

Figure 9: Algorithm of displacement integration for large faults. a: A cross-section view of a listric fault (solid line) with the global displacement constraint (V ) and an example of point (X) affected by the algorithm presented on the right part with its displacement step (v). b: The algorithm for integrating non-planar faults. Each point is progressively displaced until a global displacement target has been reached in the constraint direction. The numerical iterations increment the position of the point with a vector of limited displacement in order to ensure a good integration with the curved part of the fault. This step is repeated for each point until a portion d of the target displacement has been obtained, which represents a kinematic step.

During the numerical iterations of a kinematic step, certain points may have reached this criterion while some others are still being updated. For this reason, the geometry observed before the end of each kinematic step is not representative of a possible natural geometry, but the results of successive kinematic steps represent the evolution of the fault block during the activity of the fault.

implied by the presence of the youngest faults in the network, which is generally not considered. In the example this results in the construction of three faults (Fig. 11, Scenario 1): two faults of family 1 (yellow and blue) are first generated, then a large fault of family 2 is created (green). The yellow and blue sticks are attributed to two separate faults. Indeed, the algorithm of Cherpeau et al. (2010a) does not detect that they may correspond to a single older fault.

One can argue that the way the yellow and the green faults intersect may look inappropriate here. It is possible to adapt the geometry of the green fault or to remove the superfluous part from the yellow fault, but the yellow and the blue faults will remain considered as two distinct faults.

For better kinematic consistency, the displacement field associated with the youngest faults can be considered while simulating the older ones. For this we propose to start with the youngest faults and to go back in time while sequentially modifying the older fault information with the computed displacement fields. Applying this concept to the present example tends to align yellow and blue markers. The probability that they correspond to a single fault thus becomes non-neglegible (Fig. 11, Scenario 2).

Even if the two results are similar in terms of number of fault parts, gathering yellow and blue sticks increases the possible size of the faults of family 1. This also alters the geometry and contact length of the fault networks. Considering the displacement field would then influence the capability of faults to act as drains if the fault are permeable, or would increase the compartmentalization of the fault block if the faults are sealing.

The difference between the two approaches is even more important when considering the implication on the connectivity of the reservoir in the past. Ignoring the effect of younger faults during the clustering step produces a model that implies a rather strong connectivity between the north and south compartments. On the contrary, scenario 2 tends to consider that they were already disconnected after the appearance of the family 1 fault. Such information provides a crucial insight in the understanding of the past circulation of oil or mineralizing fluids. is interpreted here has the heritage of a fault-propagation fold [START_REF] Reches | Drag along faults[END_REF]. It is obtained by adapting the attenuation profile to produce a continuous shear parallel to the fault applied at the beginning of the deformation process. The hanging wall is progressively affected by a displacement vector field following eq. ( 2), until a global displacement of 250m to the West is reached (see algorithm Fig. 9). This process produces a roll-over anticline. It is not imposed explicitly but represents an emergent character implied by the fault geometry. The simulated roll-over anticline and the normal drag (a) are similar to those observed in the natural fault depicted in Hamblin (1965) (b). Note that the fault in this example has been chosen for its well identified normal drag and roll-over anticline. This case may not be representative of most listric faults due to the accentuated difference of dip between the top and bottom parts of the fault.

Conclusions

The method presented in this paper provides a way to parametrize in a 3-dimensional space the deformation associated with faults together with their geometry. This parametrization allows us to apply various models of fault-related displacements in order to deform geological structures affected by normal, reverse or strike-slip faults in a kinematically consistent way. The principal target applications for this parameterization are creating and editing structural models and matching structural parameters to indirect subsurface data, for example history matching [START_REF] Cherpeau | Method for stochastic inverse modeling of fault geometry and connectivity using flow data[END_REF] or potential inversion [START_REF] Guillen | Geological modelling from field data and geological knowledge: Part ii. modelling validation using gravity and magnetic data inversion[END_REF]. Such applications would gain in geological consistency by including more kinematic concepts in the construction of structural models, thereby alleviating the quality control and model editing task.

The key elements of this method lie in the explicit control of the displacements intensity, thanks to attenuation profile p, the time integration approach to the displacement computation, and the curvilinear fault frame used for describing the fault, resulting in a complete 3-dimensional volumetric displacement field. This displacement field represents a basic fault operator suitable for loosely interact-ing faults and can be altered and combined in order to fit more complex fault behavior, for example large listric faults (Fig. 9 and 10). It also offers the possibility of investigating further the time-varying aspect of the fault. For instance, syn-sedimentary faults or active structures could be better represented by progressively increasing the part of the stratigraphic sequence affected by the fault during the time-integration.

The approach presented here represents one more step towards an early integration of structural concepts in geomodeling workflows. This paradigm consists in identifying the tectonic components that impact the model, describing their effect in term of instantaneous displacement field and finally summing them by integrating the displacement field over time. The underlying hypothesis of such an approach is that tectonic features can be considered as independent, which is a rather strong assumption. However, it offers the opportunity to create first-order kinematically consistent quantitative models without relying on mechanical simulation, which generally implies a higher computational cost. Nevertheless, future work will focus on introducing interdependent tectonic features in our modeling approach in order to remove this requirement, for example by including the results of sev- The fault sticks of family 2 are unfaulted by applying the displacement field. b2: After unfaulting the yellow and blue fault sticks appear to be sufficiently aligned to be gathered in one single fault. b3: The fault of family 1 produced this way is then cut by the green fault and displaced. The resulting fault network is more consistent with the kinematics implied by the faults.

eral mechanical studies of fault interaction: faults intersection [START_REF] Maerten | Slip distributions on intersecting normal faults[END_REF], relay [START_REF] Walsh | Fault relays, bends and branch-lines[END_REF][START_REF] Soliva | Displacement and interaction of normal fault segments branched at depth: Implications for fault growth and potential earthquake rupture size[END_REF] and linkage [START_REF] Kim | The relationship between displacement and length of faults[END_REF][START_REF] Walsh | An alternative model for the growth of faults[END_REF].

Figure 2 :

 2 Figure 2: Fault flanking structures. An initially flat horizon is cut by a fault with a normal relative displacement (arrows). This superimposes a ductile deformation (thick line) to the translation of the fault blocks (dashed line). a. Normal drag: bending towards the opposite direction of the block displacement. b. Reverse drag: bending towards the direction of displacement.

Figure 3 :

 3 Figure3: The fault coordinate system. a: A block diagram with a fault cutting a horizon, C denoting the point of maximum displacement. The shaded surface represents the zone with a non-zero displacement delimited by the fault tip (plain line). b and c: The 3D curvilinear fault frame with its 3 coordinates: g 0 (red), g 1 (green), g 2 (blue). It describes an implicit mapping between the actual curvilinear geological space (c) and a regular fault space (b), which simplifies the description of the fault.

Figure 4 :

 4 Figure4: An example of attenuation profiles. The curve on the left presents the attenuation followed by p 1 and p 2 , i.e. along the fault surface. A bell-shaped profile[START_REF] Cowie | Fault tip displacement gradients and process zone dimensions[END_REF] is represented here, but other profiles could be used instead. The curve on the right represents p 0 , with a profile close to the profiles presented in[START_REF] Barnett | Displacement geometry in the volume containing a single normal fault[END_REF]. The value of the derivatives are controlled at the extremities and center of the profiles (arrows) and additional control points can be introduced in order to tune the curves (squares). Note that p 0 is represented as negative for g 0 < 0 (this behavior is managed by the factor k in eq.(2))

Figure 5 :

 5 Figure 5: Single fault example. a: The fault frame used to compute the displacement associated with this fault. b: The deformation resulting from the application of the near-field displacement on initially horizontal layers, red arrows show local displacement vectors.

Figure 6 :

 6 Figure 6: Near-field deformation on two co-genetic normal faults. The far-field has not been applied here to have a clearer view of the near-field modeling (see Fig. 8 for the complete displacement). a: Cross-section showing the stratigraphic horizons affected by the two faults. b: Top view presenting the near-field displacement norm (in meters), i.e. the magnitude of the displacement in 3 dimensions. c: View of the 3D model illustrating the displacement of the horizons affected by the two non-planar faults. Only the layers 4 to 8 are represented in order to make the internal deformation visible. Contour lines representing the elevation are drawn on the top of the layers with an interval of 35 meters.

Figure 10 :

 10 Figure 10: Listric fault displacement field. This figure presents the results of the deformation associated with a listric fault.The numerical model (a) is designed to reproduce typical folds observed in the flanks of natural listric faults (b): a normal drag, present in both blocks, and a reverse drag, in the hanging wall only, forming a roll-over anticline. The normal drag is interpreted here has the heritage of a fault-propagation fold[START_REF] Reches | Drag along faults[END_REF]. It is obtained by adapting the attenuation profile to produce a continuous shear parallel to the fault applied at the beginning of the deformation process. The hanging wall is progressively affected by a displacement vector field following eq. (2), until a global displacement of 250m to the West is reached (see algorithm Fig.9). This process produces a roll-over anticline. It is not imposed explicitly but represents an emergent character implied by the fault geometry. The simulated roll-over anticline and the normal drag (a) are similar to those observed in the natural fault depicted inHamblin (1965) (b). Note that the fault in this example has been chosen for its well identified normal drag and roll-over anticline. This case may not be representative of most listric faults due to the accentuated difference of dip between the top and bottom parts of the fault.

  b2. older faults simulation: taking account of restored faults sticks b1. youngest fault simulation (view from top): -compute associated displacement -restore older fault sticks Scenario 2: fault network simulation + sequential data restoration b3. apply fault displacement field

Figure 11 :

 11 Figure11: An example taking a fault's displacements into account in sequential fault simulation. a: Synthetic data used in this example, two cross-sections (black) with three groups of fault sticks (yellow, green and blue spheres). Scenario 1: the possible effects of the displacement field are not taken into account. Scenario 2: the older fault sticks (yellow and blue) are restored before being interpreted as a fault. b1: A fault of family 2 is modeled (green) together with its displacement field (red vectors). The fault sticks of family 2 are unfaulted by applying the displacement field. b2: After unfaulting the yellow and blue fault sticks appear to be sufficiently aligned to be gathered in one single fault. b3: The fault of family 1 produced this way is then cut by the green fault and displaced. The resulting fault network is more consistent with the kinematics implied by the faults.

Acknowledgments

The authors are grateful to the academic and industrial sponsors of the Gocad Consortium managed by ASGA (Association Scientifique pour la Géologie et ses Applications) for funding this work, and to Paradigm Geophysical for providing the Gocad software and API. We would also like to thank Dr Laurent Aillères for constructive comments which helped improve this paper. This work was performed in the framework of the "Investisse-ments d'avenir" Labex RESSOURCES21 (ANR-10-LABX-21).