Djamel Eddine Khelladi
email: djamel-eddine.khelladi@irisa.fr

Benoit Combemale
email: benoit.combemale@irisa.fr

Mathieu Acher
email: mathieu.acher@irisa.fr

Olivier Barais
email: olivier.barais@irisa.fr

Jean-Marc 2020 Jézéquel
email: jean-marc.jezequel@irisa.fr

Co-Evolving Code with Evolving Metamodels

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

Model-Driven Engineering (MDE) has proven to be effective in the development and maintenance of large scales systems [START_REF] Hutchinson | Model-driven engineering practices in industry[END_REF][START_REF] Hutchinson | Empirical assessment of MDE in industry[END_REF]. Notably, by easing the task of building Software languages (SLs) and their tooling that can play a significant role in all phases of development processes [START_REF] Tolvanen | MetaEdit+: defining and using integrated domain-specific modeling languages[END_REF].

When comes the time to specify a SL, the metamodel is pivotal and a central artifact [START_REF] Hutchinson | Empirical assessment of MDE in industry[END_REF]. It introduces the domain concepts, which are later used as a cornerstone to build the expected tooling (e.g., editor, checker, compiler, data access layers, etc.). Regular Expressions, CSS, GraphViz are examples of SLs that use an explicit initial specification to later drive the development of the corresponding tooling [START_REF] Fowler | Domain-specific languages[END_REF]. However, it is often overlooked that metamodels are inputs for complex code generators that leverage on the abstractions used in metamodels (e.g., bidirectional references, containments, . . . ) to generate complex and industrial code that may include advanced architecture and different ways to manipulate instances.

As a prominent example, Eclipse Modeling Framework (EMF) [START_REF] Steinberg | EMF: eclipse modeling framework[END_REF] supports the generation of Java code consisting of a core API (both generic and domain-specific), adapters, serialization facilities, and an editor, all from metamodels. This code can be used to load and manipulate the model instances. More importantly, it is often enriched by developers to offer additional functionalities and on top of which the SLs' toolings are further developed.

Today, many Eclipse projects 1 leverage on the code generated by EMF to support developers in various tasks (e.g., developing advanced editors, specific model transformations, etc.). As a result, many languages have been developed using EMF and its ecosystem, such as BPMN [START_REF] Omg | Object Management Group. Business Process Model And Notation (BPMN)[END_REF] or UML [START_REF] Omg | Object Management Group. Unified Modeling Language (UML)[END_REF]. In all these examples, the initial core API generated with EMF have been complemented, either manually or through other tools from the EMF ecosystem, to provide advanced functionalities and tooling such as simulation, debuging, etc. A metamodel and its corresponding implementation are hence a corner stone when building a SL and its tooling.

One of the foremost challenges to deal with in MDE is the evolution of metamodels and its impact on artifacts that use metamodels as a foundation. As a result of a metamodel evolution, the core API is re-generated and may impact the additional code implemented by developers. Hence, it may need to be co-evolved accordingly. Knowing that other projects/SLs can depend on the impacted code, it is essential to efficiently co-evolve it whenever the metamodel evolves. In the last decade, the co-evolution challenge has been extensively addressed in MDE. However, the literature so far has focused on the co-evolution between metamodel and models [START_REF] Cicchetti | Automating co-evolution in model-driven engineering[END_REF][START_REF] Garcés | Managing model adaptation by precise detection of metamodel changes[END_REF][START_REF] Herrmannsdoerfer | COPE-automating coupled evolution of metamodels and models[END_REF][START_REF] Kessentini | Automated metamodel/model co-evolution: A search-based approach[END_REF][START_REF] Kessentini | Integrating the Designer in-the-loop for Metamodel/Model Co-Evolution via Interactive Computational Search[END_REF][START_REF] Wachsmuth | Metamodel adaptation and model coadaptation[END_REF], constraints [START_REF] Batot | Heuristic-Based Recommendation for Metamodel-OCL Coevolution[END_REF][START_REF] Correa | Refactoring object constraint language specifications[END_REF][START_REF] Djamel | A semi-automatic maintenance and co-evolution of OCL constraints with (meta) model evolution[END_REF][START_REF] Kusel | Systematic Co-Evolution of OCL Expressions[END_REF], and transformations [START_REF] Garcés | Adapting transformations to metamodel changes via external transformation composition[END_REF][START_REF] García | Model transformation co-evolution: A semi-automatic approach[END_REF][START_REF] Kessentini | Automated Co-evolution of Metamodels and Transformation Rules: A Search-Based Approach[END_REF][START_REF] Djamel | Change Propagation-based and Composition-based Co-evolution of Transformations with Evolving Metamodels[END_REF][START_REF] Kusel | Consistent co-evolution of models and transformations[END_REF]. To the best of our knowledge, there exist no work that aims at supporting the co-evolution between metamodels and code.

This paper contributes in filling this gap. We propose a semi-automatic approach to co-evolution of code when the metamodels evolve. The co-evolution propagates the metamodel changes to the code. Change propagation showed to be efficient in many domains, such as in the context of artifact recommendation or co-evolution of models or transformations [START_REF] Čubranić | Hipikat: Recommending pertinent software development artifacts[END_REF][START_REF] Demuth | Co-evolution of metamodels and models through consistent change propagation[END_REF][START_REF] Djamel | Change Propagation-based and Composition-based Co-evolution of Transformations with Evolving Metamodels[END_REF]. This paper leverages on the knowledge provided by the metamodel changes during evolution which are propagated by means of code co-evolution. For example, knowing that a property p in the metamodel is moved from one class to another through a reference ref, it is likely that the navigation path to p will also be updated in the code by either extending it with ref or reducing it. The current approach first runs an impact analysis to identify all impacts at the code level.

Our approach leverages on the abstraction level of metamodels where a given metamodel element has often a complex implementation in the code with various usages of such a code. For example, for a single metamodel class, an interface and a class are generated in the code, but also a method creation in the factory class, etc. Our impact analysis is designed to identify all impacted usages of a given changed metamodel element. Then, we propose a catalog of resolutions that is used as a basis for our code co-evolution in response to metamodel evolution. Developers can select among our supported alternative resolutions.

Our approach with the proposed resolutions are evaluated on three Eclipse EMF-based SL implementations, namely OCL [START_REF] Mdt | Model Development Tools. Object Constraints Language (OCL)[END_REF], Modisco [START_REF] Mdt | Model Development Tools[END_REF], and Papyrus [START_REF] Mdt | Model Development Tools[END_REF], which have been developed for more than 10 years and have been evolved several times. We collected projects consisting of both original and evolved versions of metamodels and code. Thus, we evaluated to what extent can our approach correctly co-evolve the projects' code in response to the metamodels evolutions. Within our case studies, 477 changes were considered on five metamodels, which caused 976 impacts on 18 projects. When comparing our performed co-evolution to the expected co-evolution from the original to the evolved versions of the projects, our approach showed to support useful alternative co-evolutions that developers have manually performed. Our approach reached a weighted average (by the number of metamodels changes) of 87.4% and 88.9% respectively of precision and recall of co-evolution. Our approach also showed to be useful by maintaining imapacted code while it was unnecessary deleted. This paper is structured as follows. Section 2 introduces a background and a motivating example. Section 3 presents our overall approach. Section 4 details our evaluation results Section 5 threats to validity. Section 6 discusses related work and how our work distinguishes from code migration approaches. Finally, Section 7 concludes the paper and reflects on future works.

BACKGROUND AND MOTIVATING EXAMPLE

This section first gives a background on the role of metamodels when building SLs and their toolings and discusses the scenarios of co-evolution that arises between metamodels and code. After that, it introduces an example taken from the OCL case study in Eclipse of metamodel and code co-evolution to motivate our work.

Background

Figure 1 depicts an SL structure and usage as in practice, for instance, as in the Eclipse platform. Metamodels are cornerstones when building a SL and its toolings [START_REF] Richard | Eclipse modeling project: a domainspecific language (DSL) toolkit[END_REF][START_REF] Van Deursen | Domainspecific languages: An annotated bibliography[END_REF]. Metamodels define the aspects of a business domain, i.e. the main concepts, their properties, and relationships between them [START_REF] Cabot | Object constraint language (OCL): a definitive guide[END_REF].

Once the metamodels are carefully specified, working piece of code is generated consisting of a core API [START_REF] Steinberg | EMF: eclipse modeling framework[END_REF], which is often enriched with additional code to offer more advanced functionalities. For instance, methods defined in classes in the metamodel, are generated without their bodies, which developers must implement. This part consists of enriching the generated core API. Developers also integrate additional classes to implement advanced functionalities (language services, language tooling, . . . ). This code is built on top of the generated core API.

Take for example the Regular Expression SL. From its metamodel the generated core API allows to parse a given regular expression and to manipulate it problematically (navigate in the AST, modifies it, etc.). However, no advanced functionalities are provided, such as the interpretation and evaluation of the given regular expression. This is the tooling that developers will add on top of the core API. Both depend and uses the generated core API from the metamodel, as shown in Figure 1.

When the metamodel evolves, the same part of the generated code can be re-generated again. As a consequence, the code manually integrated by developers may be impacted and may need to be co-evolved accordingly as well. The next section illustrates these challenges on a real-world use case.

Motivating Example

To motivate our work, this section introduces an example of metamodel and code co-evolution taken from the implementation of the OMG standard Object Constraint Language (OCL) [START_REF] Mdt | Model Development Tools. Object Constraints Language (OCL)[END_REF][START_REF] Omg | Object Management Group. Object Constraints Language (OCL)[END_REF] in the Eclipse platform, which has evolved several times in the last years. OCL is a declarative language based on first-order logic, to define constraints on models, i.e., conditions that must hold on the models. Part of the OCL implementation is the OCL Pivot which provides the unified merged OCL functionality for the Ecore or UML metamodel instances. The OCL Pivot is implemented as an SL with a metamodel consisting of 111 classes in the version of 3.2.2.

Figure 2 depicts a small excerpt of the Pivot metamodel (6 out of 111 metaclasses). It illustrates the domain concepts Property, MultiplicityElement, Model, NamedElement, Namespace and ClassifierType. From all of these metaclasses, a first code is generated (i.e., core API), consisting of interfaces, classes, a factory, a package, etc. Listing 1 shows a snippet of the generated interfaces and classes for the metamodel excerpt in Figure 2 For some impacts, the co-evolution can be unique. For example, the renamed method validateCompatibleInitialiser to validateCompatibleDefaultExpression or the renamed class ClassifierType to MetaClass are propagated by developers to the code level as well and renamed in the Java class PivotValidator (lines 2 and 6 in Listing 4). However, some other impacts can be co-evolved in many different ways. Indeed, in the class Pivot2CSConversion, due to the deleted properties lower, upper, isOrdered, isUnique in the class MultiplicityElement, two different co-evolutions were applied. The navigation paths to the methods getLower() and getUpper() were deleted without deleting the whole variable declaration statement (lines 7 and 8 in Listing 3). Whereas, the whole statement of the method call refreshQualifiers is deleted rather than deleting the navigation paths to the methods getIsOrdered() and getIsUnique() (lines 12 and 13 in Listing 3). In the class PivotValidator the import and the whole method is deleted due to the delete of MultiplicityElement (lines 3 and 12 in Listing 2). Here we already see that a delete change in the metamodel can be co-evolved in different manners to cover different possible co-evolutions to different parts of the code.

Similar situation can even occur for the same impacting change. Take the pushed property ownedRule from NamedElement to Namespace. In the method refreshClassifier, the call to getOwnedRule() is replaced by getOwnedInvariant() (line 17 in Listing 3). Whereas, for the same impacted call to getOwnedRule() in the method refreshTypedElement, its whole statement of the method call visitDeclarations is deleted (line 22 in Listing 3). Other possible co-evolution could also be to introduce a cast or a cast with an if test of the type. This emphasizes the fact that developers should be involved in the co-evolution process to decide which alternatives should be applied.

Listing 1: Excerpt of the generated code in the SL org.eclipse.ocl.examples.pivot.

1 // M u l t i p l i c i t y E l e m e n t Interface 2 public i n t e r f a c e M u l t i p l i c i t y E l e m e n t extends Element

APPROACH

This section presents our approach of code co-evolution with evolving metamodels. Figure 3 shows our overall co-evolution process of code. It first identifies the metamodel changes and then runs an impact analysis on the code 1 . Here, we bridge the gap between the two abstraction levels of metamodels and code. Then a change propagation-based coevolution is applied 2 . The developer is involved throughout the co-evolution to decide among alternatives, similarly as developers would choose among quick fixes in today's IDEs to fix code errors. Finally, the chosen resolutions are applied to co-evolve the code 3 .

Change Detection and 1 : 𝑛 Impact Analysis

Before we start co-evolving code, we run an impact analysis to detect the code parts that must be co-evolved. The starting point of our impact analysis is the changes applied during the metamodel evolution. Numerous approaches exists in the literature, such as [START_REF] Cicchetti | Managing dependent changes in coupled evolution[END_REF][START_REF] Djamel | Detecting complex changes during metamodel evolution[END_REF][START_REF] Djamel | Detecting complex changes and refactorings during (meta) model evolution[END_REF][START_REF] Langer | A posteriori operation detection in evolving software models[END_REF][START_REF] Sander D Vermolen | Reconstructing complex metamodel evolution[END_REF][START_REF] James R Williams | Searching for model migration strategies[END_REF]. To not restrict our approach to a particular detection approach, we define a Change interface (see Figure 3) that serves to bridging existing change detection approaches with our internal representation of metamodel changes. Thus, we can replace a given detection approach with another one.

The change interface considers both atomic and complex changes [START_REF] Herrmannsdoerfer | An Extensive Catalog of Operators for the Coupled Evolution of Metamodels and Models[END_REF]. Atomic changes are additions, removals, and updates of a metamodel element. Complex changes consist in a sequence of atomic changes combined together. For example, push property is a complex change where a property is moved from a superclass to its subclasses. This is composed of two atomic changes: delete a property and add a property.

With the detected metamodel changes, we run an impact analysis [START_REF] Iovino | On the Impact Significance of Metamodel Evolution in MDE[END_REF] on the code. To this end, we parse the code to access the Abstract Syntax Tree (AST). We then establish a mapping between the metamodel elements 𝑒 1 , ..., 𝑒𝑛 and the code AST nodes corresponding to 𝑒 1 , ..., 𝑒𝑛, as illustrated in Algorithm 1.

In contrast to the co-evolution of metamodel and models/transformations/constraints that considers a 1 : 1 relationship between the metamodel elements and their use in the impacted artifacts. Our impact analysis considers a 1 : 𝑛 relationship where a metamodel element is used in 𝑛 different forms at the code level. Indeed, take for example the class Property in the OCL Pivot metamodel. At the code level, an interface Property and a class PropertyImpl are generated. But also in the Algorithm 1 Code impact analysis end for 27: end function classes PivotFactoryImpl and PivotPackageImpl, respectively the methods createProperty() and getProperty() are generated. All this 𝑛 code correspondences of the class Property are identified in Algorithm 1 (Lines 5-11). Other metamodel elements also have 𝑛 code correspondences, such as attributes. For example, for the attribute implicit in the class Property, a getter and a setter getImplicit()/setImplicit() are generated, but also the method getProperty_Implicit() is generated in the class PivotPackageImpl (handled in Algorithm 1 (Lines 13-20)).

Therefore, one change in the metamodel may impact its 𝑛 different usages. Our 1 : 𝑛 impact analysis leverages on the power of abstraction of the metamodels [START_REF] Djamel | On the Power of Abstraction: a Model-Driven Co-evolution Approach of Software Code[END_REF] to identify those 𝑛 impacts on the code to co-evolve. We leverage on the knowledge of how the core API code is generated to cover all 𝑛 different usages of a given metamodel element. Finally, the impact analysis consists in accessing for each impacting metamodel change on an element 𝑒 𝑖 , the set of impacted code parts and their impacted AST nodes.

Code Co-evolution

Algorithm 2 illustrates the overall co-evolution process. It first retrieves the impacted code parts by a given metamodel change from the mapping table computed by our impact for all c ∈ changes do 3:

List<CE> codeImpacts ← MT.get(c.element)

4:
for all cImpact ∈ codeImpacts do 5:

resolutions ← RC.getResolutions(cImpact, c) end for 16: end function analysis (line 3). Then for each impact, it proposes the appropriate resolutions that can propagate the impacting change (line 5). A developer can choose among the alternative resolutions which one fits her needs. This acts as a user acceptance of the resolution to be applied (line 7). Finally, the chosen resolution is applied on the impacted code part (line 13).

Table 1 depicts our resolutions that propagate to the code metamodel changes that are known to be impacting [START_REF] Iovino | On the Impact Significance of Metamodel Evolution in MDE[END_REF]. The resolutions are inspired from the co-evolution of other artifacts, such as models, constraints, or transformations [2, 5, 6, 11-13, 15, 18, 23-25, 27, 32-34, 49], where they showed to be efficient and useful. Nonetheless, they are further adapted to fit the code co-evolution.

For example, resolutions 𝐶𝑅15, 𝐶𝑅16, 𝐶𝑅17 aims to coevolve the different code impacts of changing a property type in the metamodel. To react to delete changes in the metamodels, we propose to delete different granularities of the impacted code, such as deleting the direct element, its direct parent expression (whole call path), or its instruction. Here we make sure that no error is introduced by further propagating the resolutions when necessary (e.g., to update method calls due to delete of a parameter) and by adding a cast (e.g., to the expression using a deleted property).

Note that, similarly as for code repair approaches [START_REF] Monperrus | Automatic software repair: a bibliography[END_REF], we do not claim completeness of our catalog of resolutions. It rather represents meaningful propagations of the metamodel changes at the code level. Enriching our catalog with other resolutions is left for future work. However, as we will show in section 4, the catalog of resolutions in Table 1 is sufficient and useful in practice, and allows to correctly co-evolve impacted code w.r.t. the developers intent.

Prototype Implementation

Our approach's Java implementation handles Ecore/EMF metamodels and Java code. It interfaces with our previous work of change detection [START_REF] Djamel | Ad-room: a tool for automatic detection of refactorings in object-oriented models[END_REF] to then run an impact analysis on the code. For each impacted part, we propose resolutions. The resolutions are implemented in Java. Both the impact analysis and the resolutions work on the JDT API2 to parse and to manipulate the code, application of resolutions.

EVALUATION

This section presents the evaluation results of our co-evolution approach. First, we present the dataset and evaluation process. Then, we present the research questions and discuss the obtained results. Time performance of the co-evolution is measured as well.

Data Set

This section presents the used data set in our evaluation. Our data set covers three case studies covering three different language implementations in Eclipse, namely OCL [START_REF] Mdt | Model Development Tools. Object Constraints Language (OCL)[END_REF], Modisco [START_REF] Mdt | Model Development Tools[END_REF], and Papyrus [START_REF] Mdt | Model Development Tools[END_REF].

These languages have been developed for more than 10 years and have been evolved several times. OCL is a standard language defined by the Object Management Group (OMG) to specify First-order logic constraints. Modisco is an Academic initiative to support development of model-driven tools, verification, and transformation of existing software systems. Papyrus is an industrial project led by CEA3 to support model-based simulation, formal testing, safety analysis, etc. Thus, our data set covers standard, academic, and industrial languages. We considered original and evolved versions of the case studies, as shown in Table 2.

The three case studies cover code co-evolution in response to five evolved metamodels respectively, two in the OCL language, one in the Modisco Language, and two in the Papyrus language, as shown in Table 2. The total of applied metamodel changes was 477 in the five metamodels.

For those three case studies, we collected 15 Java projects that were impacted by those five evolving metamodels (i.e., projects that are dependent on the metamodels' generated core API). We collected the original and evolved Java code of those projects. Three projects were impacted by two metamodels. Thus, in total we considered 18 projects' co-evolutions due to the five evolved metamodels.

All evolutions in our case studies of the metamodels and code have been performed manually by developers between the different releases that we collected. Thus, we considered the manual code co-evolution as the reference (i.e., ground truth) to which we compare our performed co-evolution.

Table 2 gives details about the selected case studies in particular about their metamodels and the applied changes during evolution. Table 3 further gives details on the size of the projects and code of the original versions that we co-evolve.

Evaluation Process

We evaluate our change propagation-based co-evolution of code by measuring the correctness of our co-evolution. We compare for the same set of code how it was manually coevolved by developers against our proposed co-evolution. In this experiment, we measure the correctness of our approach by using the two metrics precision and recall that vary from 0 to 1, i.e., 0% to 100%. They are defined as follows:

Table 3: Details of the projects and their impacts caused by the metamodels evolution.

Evolved metamodels

Projects to co-evolve in response to the evolved metamodels

𝑁 𝑜 of packages 𝑁 𝑜 of classes 𝑁 𝑜 of LOC 𝑁 𝑜 of
𝑝𝑟𝑒𝑐𝑖𝑠𝑖𝑜𝑛 = 𝑃 𝑟𝑜𝑝𝑜𝑠𝑒𝑑𝑅𝑒𝑠𝑜𝑙𝑢𝑡𝑖𝑜𝑛𝑠 ∩ 𝐸𝑥𝑝𝑒𝑐𝑡𝑒𝑑𝑅𝑒𝑠𝑜𝑙𝑢𝑡𝑖𝑜𝑛𝑠 𝑃 𝑟𝑜𝑝𝑜𝑠𝑒𝑑𝑅𝑒𝑠𝑜𝑙𝑢𝑡𝑖𝑜𝑛𝑠 𝑟𝑒𝑐𝑎𝑙𝑙 = 𝑃 𝑟𝑜𝑝𝑜𝑠𝑒𝑑𝑅𝑒𝑠𝑜𝑙𝑢𝑡𝑖𝑜𝑛𝑠 ∩ 𝐸𝑥𝑝𝑒𝑐𝑡𝑒𝑑𝑅𝑒𝑠𝑜𝑙𝑢𝑡𝑖𝑜𝑛𝑠 𝐸𝑥𝑝𝑒𝑐𝑡𝑒𝑑𝑅𝑒𝑠𝑜𝑙𝑢𝑡𝑖𝑜𝑛𝑠
The 𝑃 𝑟𝑜𝑝𝑜𝑠𝑒𝑑𝑅𝑒𝑠𝑜𝑙𝑢𝑡𝑖𝑜𝑛𝑠 are the resolutions applied by our approach (from Table 1) and the 𝐸𝑥𝑝𝑒𝑐𝑡𝑒𝑑𝑅𝑒𝑠𝑜𝑙𝑢𝑡𝑖𝑜𝑛𝑠 are the actual manually performed resolutions by developers.

Note that we also studied to confirm the metamodel changes from version 1 to version 2. Based on those changes we run an impact analysis to identify the impacted parts of the code that will be co-evolved.

Research Questions

This section defines two research questions (RQs) to assess in particular the applicability, correctness, and the usefulness of our work. The research questions are as follows:

RQ1: To what extent and how fast can our co-evolution approach handle the impacted parts of the code? If we cannot propose a co-evolution for all impacted parts due to evolving metamodels, we can conclude that our change propagationbased co-evolution of code is insufficient. This assesses the overall applicability of our co-evolution approach.

RQ2: To what extent does our approach correctly co-evolve the impacted code? This measures its correctness level with the precision and recall metrics, i.e., to which extent it proposes the expected resolutions from developers. It also allows us to assess its usefulness. For alternative resolutions, we select the one as close as possible to the expected resolution to be able to compare them afterward.

RQ3: Can our approach still propose relevant co-evolution when it fails to propose the expected one by developers? This aims to investigate, when our co-evolution diverges from the expected one, whether our proposed resolutions are still useful to developers.

Results

We now discuss the results w.r.t. our research questions.

RQ1.

With our impact analysis, we were able to find the impacted code parts by the 477 metamodel changes for which resolutions were proposed by our co-evolution approach. One to five resolutions were proposed for each impact in the code. To co-evolve all impacted parts, a total of 983 resolutions were applied during the change propagation. This shows the applicability of our co-evolution approach that was able to handle all different impacts in the code caused by the metamodel evolution changes. Table 4 gives the list of the applied resolutions for each co-evolved project during co-evolution. In total, the 972 applied resolutions were covered by 11 resolution types from our catalog in Table 1, namely 𝐶𝑅1, 𝐶𝑅2, 𝐶𝑅3, 𝐶𝑅5, 𝐶𝑅6, 𝐶𝑅6, 𝐶𝑅9, 𝐶𝑅11, 𝐶𝑅12, 𝐶𝑅15, 𝐶𝑅16, 𝐶𝑅17. These resolutions allowed us to cover different possible co-evolution for the different impacting metamodel changes. Note that all occurrences of 𝐶𝑅1 were due to deleted classes (DC) in the metamodels that were used in the code by subclasses (i.e., extends DC). Figure 4 further shows the frequency of each applied resolution. As some impacting change types did not occur in our case studies, their resolutions were never applied, such as 𝐶𝑅7, 𝐶𝑅8, 𝐶𝑅14.

The parsing of the project's code took from couple of seconds to a minute for the largest project. For each impacting metamodel change, the detection of the impacted code parts took from milliseconds to few seconds, and the application of each resolution was a matter of less than half a second. The evaluation was run on a Windows 7 PC with a Core i7 3.4GHz and 16GB RAM.

RQ2.

The metamodel changes present in our case studies (see Table 2 columns 5 and6) covered different types of changes and evolution scenarios. In particular the impacting changes for which we proposed resolutions in Table 1. To assess to what extent can our approach correctly co-evolve the impacted code w.r.t. the manual co-evolution of the code that developers went through, we measured precision and recall.

Figure 5 depict the reached precision and recall in our 18 scenarios of metamodel evolution and code co-evolution. The projects [P9] and [P16] were not impacted by any metamodel evolution with no performed co-evolution. Thus, we did not measure precision and recall for them. For the rest of the projects, both the measured precision and recall varied from 48% to 100%. Our code co-evolution was able to reach respectively the average of 94,5% and 95,5% of precision and recall. However, as the size of the projects and the metamodels' evolution varied from small to large, we computed a weighted average to better characterize precision and recall. We computed the weighted average by the number of the metamodel changes for each project, resulting in a weighted average of 87.4% and 88.9% of precision and recall. This means that our applied resolutions propagating the impacting metamodel changes covered the expected resolutions by the developers in 88.9% and were correct in 87.4%. This shows that the catalog of resolutions in Table 1 showed to be both sufficient and useful in co-evolving code in our case studies while meeting the developer's co-evolution needs.

RQ3.

We further looked into the cases where our proposed resolutions mismatched the expected resolutions, i.e., the cause of lowering precision and recall. We observed that for all those impacts in the code caused by non-delete changes in the metamodels, the expected resolutions applied in the co-evolution by the developers were to delete the impacted code. Rather than to delete them, our co-evolution approach was able to successfully maintain them.

For example, in the project [P12], we observed that 28 impacted parts of the code were due to move and rename changes in the metamodel. Listing 6 shows an excerpt of two impacted part of the code from the project [P12] due to: 1) rename discoveryError to discoveryErrors, and 2) Our co-evolution approach was able to maintain impacted parts of the code, whereas, the developers' manual co-evolution consisted in deleting them as the whole project was deleted in the new version. The applied resolutions were replacing the call to discoveryDate to a call from its new class and renaming discoveryError, as shown in Listing 8.

THREATS TO VALIDITY

This section discusses threats to validity [START_REF] Wohlin | Experimentation in software engineering[END_REF].

Internal Validity

To measure precision and recall, we had to analyze the developers manual co-evolution to identify the expected resolutions. To reduce the risk of misidentifying an expected resolution, for each impacted part of the code, we investigated the whole co-evolved class. In case we did not find it, we further searched in other classes in case the original impacted part of the code was moved into another class. Thus, we aimed at reducing the risk of missing any correspondence between the original and evolved impacted parts of the code. Moreover, as our co-evolution relies on the quality of detected metamodel changes. We analyzed each detected change and checked whether it occurred between the original and evolved metamodels to alleviate the risk of relying on an incorrect metamodel change. Moreover, in the evaluation, when alternative resolutions were proposed during co-evolution, we selected the one as close as possible to the expected resolution to be able to compare them afterward. Our goal in the evaluation was to assess (RQ1) whether we can co-evolve the impacted code and (RQ2) usefulness of our approach in terms of how far can it allow a developer to reach the same manual co-evolution of our case studies. Assessing how different developers would co-evolve the same impacts, possibly diverging with different resolutions, was not the goal of our evaluation. It is left for a future user study where we would also investigate what influences developers to select a given resolution over another one. Still, We observed a start of an answer in (RQ3) where alternative resolutions were actually relevant.

External Validity

We implemented and evaluated our approach for EMF/Ecore metamodels and Java code. Other languages, such as C# or C++, use different syntax but conceptually use the same constructions as in Java. Although we think that the coevolution would be applicable for other languages, we cannot generalize our results. Further experimentation on other languages is necessary. However, the only requirement to apply our approach to other languages is to have access to the ASTs of the parsed code and to adapt our resolutions to the new ASTs' structure. Furthermore, our evaluation was performed on SLs and their toolings in the Eclipse platform. Thus, we cannot generalize our findings to all SLs or DSLs neither domain models for instance. However, our approach could be used to coevolve Java code built on a generated API from a metamodel or a domain model. This is the case for example for tools such as JHipster that generates a complete and modern Web app or microservice architecture based on a domain model and a set of architectural choices. Further evaluations remain necessary here.

Conclusion Validity

Our evaluation gave promising results, showing that our code co-evolution is fast and useful with a weighted average precision and recall respectively of 87.4% and 88.9%. The evaluation results also showed the usefulness of our approach and the proposed resolutions in our catalog in Table 1. Even though we evaluated it on 18 scenarios of metamodel evolution and code co-evolution. To have more insights and statistical evidence, further evaluation is needed on more case studies.

RELATED WORK

The main idea of change propagation was investigated for different purposes, such as by Cubranic et al. [START_REF] Čubranić | Hipikat: Recommending pertinent software development artifacts[END_REF] who used it to recommend relevant software development artifacts, or Demuth el al. [START_REF] Demuth | Co-evolution of metamodels and models through consistent change propagation[END_REF] who used it for models co-evolution. In our work we use change propagation to co-evolve code.

In this section, we present the main related work w.r.t. co-evolution. Many approaches proposed to co-evolve models [START_REF] Cicchetti | Automating co-evolution in model-driven engineering[END_REF][START_REF] Garcés | Managing model adaptation by precise detection of metamodel changes[END_REF][START_REF] Hebig | Approaches to co-evolution of metamodels and models: A survey[END_REF][START_REF] Herrmannsdoerfer | COPE-automating coupled evolution of metamodels and models[END_REF][START_REF] Kessentini | Automated metamodel/model co-evolution: A search-based approach[END_REF][START_REF] Kessentini | Integrating the Designer in-the-loop for Metamodel/Model Co-Evolution via Interactive Computational Search[END_REF][START_REF] Richard F Paige | Evolving models in model-driven engineering: State-of-the-art and future challenges[END_REF][START_REF] Wachsmuth | Metamodel adaptation and model coadaptation[END_REF], constraints [START_REF] Batot | Heuristic-Based Recommendation for Metamodel-OCL Coevolution[END_REF][START_REF] Correa | Refactoring object constraint language specifications[END_REF][START_REF] Djamel | A semi-automatic maintenance and co-evolution of OCL constraints with (meta) model evolution[END_REF][START_REF] Djamel | Metamodel and constraints co-evolution: A semi automatic maintenance of ocl constraints[END_REF][START_REF] Kusel | Systematic Co-Evolution of OCL Expressions[END_REF], and model transformations [START_REF] Garcés | Adapting transformations to metamodel changes via external transformation composition[END_REF][START_REF] García | Model transformation co-evolution: A semi-automatic approach[END_REF][START_REF] Kessentini | Automated Co-evolution of Metamodels and Transformation Rules: A Search-Based Approach[END_REF][START_REF] Djamel | Change Propagation-based and Composition-based Co-evolution of Transformations with Evolving Metamodels[END_REF][START_REF] Kusel | Consistent co-evolution of models and transformations[END_REF]. Our work first distinguishes by co-evolving code. It further distinguishes with related work by performing a 1 : 𝑛 impact analysis. Whereas, a 1 : 1 impact analysis is performed in the approaches of models/constraints/transformations co-evolution [2, 5, 6, 11-13, 18, 23-25, 27, 32-34, 49].

In this paper, we focus so far the co-evolution on the direction metamodel to code, which is not trivial. Yu et al. [START_REF] Yu | Maintaining invariant traceability through bidirectional transformations[END_REF] proposes to co-evolve the metamodels and the generated API in both directions. However, they do not co-evolve the additional code on top of it, which our approach does.

Existing approaches of code migration are related to our work. We focus on the main existing approaches. Henkel et al. [START_REF] Henkel | CatchUp! Capturing and replaying refactorings to support API evolution[END_REF] proposed an approach that captures refactoring actions and replay them on the code to migrate. However, they support only the changes renames, moves, and type changes.

Nguyen et al. [START_REF] Hoan | A graph-based approach to API usage adaptation[END_REF] also proposed an approach that guides developers in adapting code by learning adaptation patterns from previously migrated code. Similarly, Dagenais et al. [START_REF] Dagenais | Recommending adaptive changes for framework evolution[END_REF] also uses a recommendation mechanism of code changes by mining them from previously migrated code. Anderson et al. [START_REF] Andersen | Generic patch inference[END_REF] proposed to migrate drivers in response to evolutions in Linux internal libraries. It identifies common changes made in a set of files to extract a generic patch that can be reused on other code parts.

Our current work distinguishes from these code migration approaches [START_REF] Andersen | Generic patch inference[END_REF][START_REF] Dagenais | Recommending adaptive changes for framework evolution[END_REF][START_REF] Henkel | CatchUp! Capturing and replaying refactorings to support API evolution[END_REF][START_REF] Hoan | A graph-based approach to API usage adaptation[END_REF] by considering and reasoning on the changes at the metamodel level and not at the code level. Thus, our approach treats way less changes to identify the impacted parts in the code than code migration approaches [START_REF] Andersen | Generic patch inference[END_REF][START_REF] Dagenais | Recommending adaptive changes for framework evolution[END_REF][START_REF] Henkel | CatchUp! Capturing and replaying refactorings to support API evolution[END_REF][START_REF] Hoan | A graph-based approach to API usage adaptation[END_REF]. This is possible thanks to the abstraction offered by the metamodels. Whereas our work considers 𝑜𝑛𝑒 change of a given element, 𝑛 changes must be considered in order to fully migrate the code. To the best of our knowledge, our work is the first attempt to address the challenge of code co-evolution with evolving metamodels in SLs.

CONCLUSION

This paper proposed a code co-evolution approach when metamodels evolve. It supports a change propagation-based co-evolution by leverages on the metamodel changes and propagates them on the code. Our code co-evolution was evaluated on 15 projects and five metamodels from three Eclipse language implementations and it was stress tested on 18 scenarios of metamodel evolution and code co-evolution. It showed to be efficient and useful in co-evolving code with a weighted average precision and recall respectively of 87.4% and 88.9%.

As future work, we first plan to enrich our catalogue of resolutions. Although, we envision and present our resolutions similarly as code quick fixes or repair for developers to choose from, we plan to explore ranking heuristics to ease the task of selecting resolutions for developers. Building on that, we can then propose a full automation based on the highest ranked resolutions to compute co-evolutions plans. Finally, we will evaluate our approach on more case studies other than Eclipse language implementations.

Figure 1 :

 1 Figure 1: An SL structure and evolution in the Eclipse platform.

  for the OCL Pivot SL. It is further enriched by the developers with a validation functionality in the new class PivotValidator, as shown in Listing 2. Another SL, namely OCL Base, also leverages on the OCL Pivot and its generated interfaces and classes in a conversion functionality from OCL Pivot to OCL Base, as shown in Listing 3. In version 3.4.4, the OCL Pivot metamodel evolved significantly with various changes, among which, developers: 1) deleted the metaclasses MultiplicityElement and Model along with their properties, 2) renamed the method validateCompatibleInitialiser in the class Property to val-idateCompatibleDefaultExpression, 3) renamed the class ClassifierType to MetaClass, and 4) pushed the property ownedRule from NamedElement to Namespace. Naturally, the Code of Listing 1 was re-generated from the evolved version of the metamodel, which impacted the existing code depicted in Listings 2 and 3. The different caused impacts are underlined in red and their co-evolution is underlined in green in Listings 4 and 5 corresponding to version 3.4.4.

Figure 2 :

 2 Figure 2: Excerpt of the Pivot metamodel.

{ ... } 3 /

 3 / Property Interface 4 public i n t e r f a c e Property extends Feature , ... { 5 ...

6 boolean

 6 v a l i d a t e C o m p a t i b l e I n i t i a l i s e r ( D ia gn os t ic Ch ai n diagnostics , Map < Object , Object > context ) ; 7 ...

8 } 9 /}Listing 2 :Listing 3 :} p r o t e c t e dListing 4 :

 89234 / PropertyImpl Class 10 public class PropertyImpl i m p l e m e n t s Property { 11 public boolean v a l i d a t e C o m p a t i b l e I n i t i a l i s e r ( Diagnostic Ch ai n diagnostics , Map < Object , Object > context ) { try {...} catch ( I n v a l i d V a l u e E x c e p t i o n e ) {...} } Excerpt of the additional internal code in the SL org.eclipse.ocl.examples.pivot.// Provides validation support import org.eclipse.ocl.examples.pivot.ClassifierType; import org.eclipse.ocl.examples.pivot.MultiplicityElement; public class Pivo tValidat or extends O bj ec tV a li da to r { ... public boolean v a l i d a t e P r o p e r t y _ v a l i d a t e C o m p a t i b l e I n i t i a l i s e r (...) { return property .validateCompatibleInitialiser( dia , con ) ; } public boolean v a l i d a t e C l a s s i f i e r T y p e (Classifier Type classifierType , ...) { ... } public boolean v a l i d a t e M u l t i p l i c i t y E l e m e n t (Multi plicityElement multiplicityElement, ...) { ... } } Excerpt of the additional external code in the SL org.eclipse.ocl.examples.Base. // Supports conversion of OCL Pivot to OCL Base public class P iv o t 2 C S C o n v e r s i o n extends ... { ... public TypedElement CS r e f r e s h T y p e d M u l t i p l i c i t y E l e m e n t (...) { ... if ( csTypeRef != null ) { int lower = object .getLower(). intValue () ; int upper = object .getUpper(). intValue () ; a l i f i e r s (... , object .getIsOrdered() ? Boolean . TRUE : null ) ; re fre shQu a l i f i e r s (... , object .getIsUnique() ? null : Boolean . FALSE ) ; return csElement ; ClassifierCS r e f r e s h C l a s s i f i e r (...) { ... refreshList ( csElement . g e t O w n e d C o n s t r a i n t () , vi sit D e c l a r a t i o n s ( ConstraintCS . class , object .getOwnedRule(), null ) ) ; ... } public TypedElement CS r e f r e s h T y p e d E l e m e n t (...) { ... refreshList ( csElement . g e t O w n e d C o n s t r a i n t () , vi sit D e c l a r a t i o n s ( ConstraintCS . class , object .getOwnedRule(), null ) ) ; return csElement ; } } Excerpt of the additional internal code in the SL org.eclipse.ocl.examples.pivot.

19 } 20 21 public

 192021 Type dElement CS r e f r e s h T y p e d E l e m e nt (...) { 22

Figure 3 :

 3 Figure 3: Overall approach [𝐼𝑛: an impact in the code, 𝑅𝑛: a resolution].

Algorithm 2

 2 Code co-evolution 1: function Code_Co-Evolution(Set<changes> changes, Map<MME, List<CE>> MT, ResolutionCatalog RC) 2:

Figure 4 :

 4 Figure 4: Frequency of the applied resolutions.

Figure 5 :

 5 Figure 5: Precision measurement for our code co-evolution.

  for ( Constraint constraint : (Namespace namedElement).getOwnedRule()) { 10

  1 // Provides validation support 2 import org.eclipse.ocl.examples.pivot..MetaClass; 3 public class Pivot Validat or extends O bj ec tValidator { 4

		...
	5	public boolean
		v a l i d a t e P r o p e r t y _ v a l i d a t e C o m p a t i b l e I n i t i a l i s e r
		(...) {
	6	return property .validateCompatibleDefaultExpres
		sion( dia , con ) ;
	7	}
	8
	9	public boolean v a l i d a t e C l a s s i f i e r T y p e (MetaCLass
		metacLass , ...) { ... }
	10 int lower;
	8	int upper;
	9	...
	10	}
	11	...
	12	return csElement ;
	13

} Listing 5: Excerpt of the additional external code in the SL org.eclipse.ocl.examples.Base. 1 // Supports conversion of OCL Pivot to OCL Base 2 public class P i v o t 2 C S C o n v e r s i o n extends A b s t r a c t C o n v e r s i o n i m p l e m e n t s Pivot Constants { 3 ...

4 public Type dElement CS r e f r e s h T y p e d M u l t i p l i c i t y E l e m e n t (...) { 5 ... 6 if ( csTypeRef != null ) { 7 } 14 15 p r o t e c t e d ClassifierCS r e f r e s h C l a s s i f i e r (...) { 16 ...

17

refreshList ( csElement . g e t O w n e d C o n s t r ai nt () , v i s i t D e c l a r a t i o n s ( ConstraintCS . class , object .getOwnedInvariant(), null ) ) ; 18 ...

  if (ce 𝑖 .name = mme 𝑖 .name ∨ (ce 𝑖 .name = "create"+mme 𝑖 .name ∧ ce 𝑖 .class.isFactory()) . . . ) then

	4:	List<CE> cList;
	7:
	10:
	11:	cList.add(ce 𝑖 )
	12:	end if
	13:	else
	14:
	19:	end if
	20:
	22:	end if
	23:	end if
	24:
	25:	end for
	26:

1: function ConstructMappingTable(Set<MME> mmElements, Set<CE> cElements)

2: Map<MME, List<CE>> mTable; 3: for all mme 𝑖 ∈ mmElements do 5: for all ce 𝑖 ∈ cElements do 6: if mme 𝑖 .isTypeOf(Class) then 8: ce 𝑖 .name = mme 𝑖 .name+"Impl" ∨ 9: if mme 𝑖 .isTypeOf(Attribute) then 15: if ( (ce 𝑖 .name=mme 𝑖 .name ∨ 16: ce 𝑖 .name="get"+mme 𝑖 .name ∨ 17: ce 𝑖 .name='set'+mme 𝑖 .name) ∧ 18: ce 𝑖 .class=mme 𝑖 .class) . . . ) then cList.add(ce 𝑖 ) 21: else . . . //Operations, References, Ennum, etc. mTable.put(mme 𝑖 , cList)

Table 1 :

 1 Catalog of resolutions that propagates the metamodel changes. ▷𝐶𝑅7 Introduce a for loop statement to iterate on a collection (e.g., lng.p.value → for(v in lng.p) {v.val}) ▷𝐶𝑅8 Retrieve the first value of a collection (e.g., value = lng.p → value = lng.p.toArray()[0]) ◇ Move property 𝑝 𝑖 from class S to T through ref ◇ Extract class of properties 𝑝 1 , ..., 𝑝𝑛 from S to T through ref ▷𝐶𝑅9 Extend navigation path of 𝑝 𝑖 (e.g., lng.𝑝 𝑖 → lng.ref.𝑝 𝑖 ) ▷𝐶𝑅10 Reduce navigation path of 𝑝 𝑖 (e.g., lng.ref.𝑝 𝑖 → lng.𝑝 𝑖 ) ▷𝐶𝑅11 Extend navigation path of 𝑝 𝑖 and add a for loop (e.g., lng.𝑝 𝑖 → for(v in lng.ref) {v.𝑝 𝑖 }) ▷𝐶𝑅5 Replace 𝑝 𝑖 with another property 𝑝 ′

	Impacting Metamodel Changes	Proposed Code Resolutions
		▷𝐶𝑅1 Remove the direct use of p (e.g., label = s.name + s.m1().p.m2() → label = s.name + ( (Type_Of_P) s.m1() ).m2())
		▷𝐶𝑅2 Remove the statement using p (i.e., if, loop, assignment, etc.)
	◇ Delete property p	▷𝐶𝑅3 Remove the whole call path of p (e.g., label = s.name + s.m1().m2().p → label = s.name)
		▷𝐶𝑅4 Replace the whole call path of p with a default value (e.g., id = s.id + s.m1().m2().p → id = s.id + 0)
		▷𝐶𝑅5 Replace p with another property p' (can be given as input by the developer)
		▷𝐶𝑅1 Remove the direct use of the type c (e.g., extending/implementing c, in method argument/returned
		type and not the whole method declaration. Calls to the updated methods are subsequently updated)
	◇ Delete class C	▷𝐶𝑅2 Remove the statements using the type C (e.g., import, variable declaration, method argument/returned type,
		method declaration, type instantiation, etc. Calls to the deleted variables and methods are subsequently removed)
		▷𝐶𝑅5 Replace C with another type C'
	◇ Rename element e	▷𝐶𝑅6 Rename e in the code
	◇ Generalize property p
	multiplicity from a single
	value to multiple values

𝑖 ◇ Push property p from class Sup to Sub 1 ,...,Sub𝑛 ▷𝐶𝑅12 Introduce a type test with an If statement (e.g., t.name = s.p.name → 𝑖𝑓 𝑠.𝑝.𝑖𝑠𝑡𝑦𝑝𝑒𝑜𝑓 𝑆𝑢𝑏 1 ) {t.name = (Sub 1 s).p.name} ... 𝑒𝑙𝑠𝑒 𝑖𝑓 𝑠.𝑝.𝑖𝑠𝑡𝑦𝑝𝑒𝑜𝑓 𝑆𝑢𝑏𝑛 {t.name = (Sub𝑛 s).p.name}) ▷𝐶𝑅13 Cast p to one specific sub class 𝑆𝑢𝑏 𝑖 (e.g., t.name = s.p.name → t.name = ((𝑆𝑢𝑏 𝑖 )s).p.name) ▷𝐶𝑅5 Replace p with another property p' ◇ Inline class S to T with properties 𝑝 1 , ..., 𝑝𝑛 ▷𝐶𝑅14 Change the class type from S to T (e.g., List<S> l = ...; → List<T> l = ...; ) ▷𝐶𝑅10 Reduce navigation path of 𝑝 𝑖 (e.g., lng.ref.𝑝 𝑖 → lng.𝑝 𝑖 )

◇ Change property p type from S to T ▷𝐶𝑅15 Change variable declaration type initialized with p from S to T (e.g., S var = s.p; → T var = s.p;) ▷𝐶𝑅16 Add a cast of p when used as argument of a method (e.g., mCall(p) → mCall((ExpectedType) p)) ▷𝐶𝑅17 Cast the value for the setter of p to T (e.g., s.setP(v); → s.setP((T) v);)

Table 2 :

 2 Details of the metamodels and their evolutions.

	Case study	Evolved metamodels	Versions	Atomic changes in the metamodel	Complex changes in the metamodel
				Deletes: 2 classes, 16 properties, 6 super types
		Pivot.ecore in project ocl.examples.pivot	3.2.2 to 3.4.4	Renames: 1 class, 5 properties Property changes: 4 types; 2 multiplicities Adds: 25 classes, 121 properties, 36 super types	1 pull property 2 push properties Total = 3
	OCL			Total = 218
				Deletes: 4 classes, 11 properties, 17 super types	5 moves properties
		Base.ecore in project ocl.examples.xtext.base	3.1.0 to 3.4.4	Renames: 1 property Property changes: 2 type Adds: 10 classes, 38 properties, 19 super types	1 push property 2 extract class 2 extract super class
				Total = 102	Total = 10
	Modisco	Benchmark.ecore in project modisco.infra.discovery.benchmark	0.9.0 to 0.13.0	Deletes: 6 classes, 19 properties, 5 super types Renames: 5 properties Adds: 7 classes, 24 properties, 4 super types Total = 70 changes	4 moves property 6 pull property 1 extract class 1 extract super class Total = 12
				Deletes: 10 properties, 2 super types	2 pull property
		ExtendedTypes.ecore in project	0.9.0 to	Renames: 3 classes, 2 properies	1 push property
		papyrus.infra.extendedtypes	1.1.0	Adds: 8 classes, 9 properties, 8 super types	1 extract super class
	Papyrus			Total = 42	Total = 3
		Configuration.ecore in project papyrus.infra.queries.core.configuration	0.9.0 to 1.1.0	Deletes: 6 classes, 7 properties, 4 super types Total = 17	none

  𝑃 9 modisco.infra.discovery.benchmark 𝑃 10 gmt.modisco.java.discoverer.benchmark 𝑃 11 modisco.java.discoverer.benchmark 𝑃 12 modisco.java.discoverer.benchmark.javaBenchmark

						Impacted	𝑁 𝑜 of total
						classes	impacts
	OCL	𝑃 1 ocl.examples.pivot	22	439	74002	47
	Pivot.ecore	𝑃 2 ocl.examples.xtext.base	12	181	17599	9
		𝑃 3 ocl.examples.xtext.base	12	181	17599	21
		𝑃 4 ocl.examples.xtext.completeocl	14	83	33807	9
	OCL	𝑃 5 ocl.examples.xtext.essentialocl	17	135	30920	15
	Base.ecore	𝑃 6 ocl.examples.xtext.oclinecore	15	46	45862	3
		𝑃 7 ocl.examples.xtext.oclinecore.ui	10	25	55648	3
		𝑃 8 ocl.examples.xtext.oclstdlib	14	58	38639	1
			3	28	2333	0	0
	Modisco		8	21	1947	3	13
	Benchmark.ecore		10	28	2794	5	24
			3	16	1654	7	58
	Papyrus ExtendedTypes.ecore		7 5 5	37 12 15	2057 374 725	8 5 5	47 14 14
	Papyrus Configuration.ecore		4 7 5	23 37 12	1045 2057 374	0 4 4	0 12 18

𝑃 13 papyrus.infra.extendedtypes 𝑃 14 papyrus.infra.extendedtypes.emf 𝑃 15 papyrus.uml.tools.extendedtypes 𝑃 16 papyrus.infra.queries.core.configuration 𝑃 17 papyrus.infra.extendedtypes 𝑃 18 papyrus.infra.extendedtypes.emf

Table 4 :

 4 Number of applied resolutions in our code coevolution for each project and per evolved metamodel.

	Evolved	Co-evolved	𝑁 𝑜 of applied
	metamodels	projects	resolutions
			𝐶𝑅1 : 3, 𝐶𝑅2 : 57, 𝐶𝑅3 : 298,
		𝑃 1	𝐶𝑅5 : 19, 𝐶𝑅6 : 108, 𝐶𝑅12 : 18,
	OCL		𝐶𝑅15 : 15, 𝐶𝑅16 : 3, 𝐶𝑅17 : 11
	Pivot.ecore		𝐶𝑅2 : 10, 𝐶𝑅3 : 8, 𝐶𝑅5 : 2,
		𝑃 2	𝐶𝑅6 : 5, 𝐶𝑅12 : 2, 𝐶𝑅15 : 1
			𝐶𝑅16 : 1, 𝐶𝑅17 : 1
			𝐶𝑅1 : 3, 𝐶𝑅2 : 91, 𝐶𝑅3 : 5,
		𝑃 3	𝐶𝑅5 : 18, 𝐶𝑅6 : 1, 𝐶𝑅9 : 10, 𝐶𝑅11 : 1, 𝐶𝑅12 : 1, 𝐶𝑅15 : 7
			𝐶𝑅16 : 3, 𝐶𝑅17 : 1
	OCL	𝑃 4	𝐶𝑅2 : 10, 𝐶𝑅5 : 6
	Base.ecore	𝑃 5	𝐶𝑅1 : 1, 𝐶𝑅2 : 25, 𝐶𝑅5 : 17, 𝐶𝑅9 : 1, 𝐶𝑅11 : 1, 𝐶𝑅12 : 1
		𝑃 6	𝐶𝑅2 : 9, 𝐶𝑅3 : 1
		𝑃 7	𝐶𝑅1 : 1, 𝐶𝑅2 : 5
		𝑃 8	𝐶𝑅5 : 2
		𝑃 9	𝑛𝑎
	Modisco	𝑃 10	𝐶𝑅2 : 13
	Benchmark. ecore	𝑃 11 𝑃 12	𝐶𝑅2 : 24 𝐶𝑅1 : 2, 𝐶𝑅2 : 28, 𝐶𝑅5 : 20, 𝐶𝑅6 : 8
	Papyrus	𝑃 13	𝐶𝑅2 : 1, 𝐶𝑅6 : 47
	ExtendedTypes.	𝑃 14	𝐶𝑅6 : 16
	ecore	𝑃 15	𝐶𝑅6 : 14
	Papyrus	𝑃 16	𝑛𝑎
	Configuration.	𝑃 17	𝐶𝑅2 : 12
	ecore	𝑃 18	𝐶𝑅2 : 18

Eclipse Java development tools (JDT): https://www.eclipse.org/jdt/core/

http://www-list.cea.fr/en/

ACKNOWLEDGMENTS

The research leading to these results has received funding from the CNRS PEPS, and from the AIS Rennes Métropole under grant no. 190270.

Similarly, in the project [P1], we observed the same cases of 41 deleted code impacted by the metamodel changes rename, push property, and type change. Listing 7 shows an excerpt of two impacted part of the code from the project [P1] due to: 1) push property ownedRule and 2) change type of the property specification. The developer manual co-evolutions consisted in deleting the statements where they were used. Instead, we maintained them by changing the type of the variable declaration and by introducing an if statement with a cast, as shown in Listing 9.

The above observations show that developers unnecessarily deleted many impacted parts in the code, while our approach instead successfully maintain them during co-evolution. This is the cause of lowering precision and recall in our case studies, e.g., 48% in the project [P12]. From our point of view, it is surprising to delete code statements, whereas it would have been possible to maintain them (e.g., with rename, move, or push resolutions). This obviously hints on the lack of co-evolution support. Otherwise, they would have been easily maintained in the new version. Another explanation is that the code functionality (or its requirement) became unnecessary in the new version, which explains its deletion.

Whereas these cases lower precision and recall, we do not consider it as wrong co-evolution. It could actually be preferred over simply deleting the impacted code. Indeed, it is more relevant to propagate a rename of a method than to delete its call elsewhere in the code, which our approach was able to achieve for those cases. Nonetheless, it could be possible that for some case simply deleting the impacted code (by non-delete metamodel changes) is desired. In future work, we could propose delete resolutions as an alternative for any impacted code. That is why it is essential for developers to guide the co-evolution process and decide which resolutions best fit their needs.


Nicolas Cherpeau
email: nicolas.cherpeau@gocad.org

Guillaume Caumon

Jef Caers

Bruno Lévy

Method for stochastic inverse modeling of fault geometry and connectivity using flow data *

Keywords: structural modeling, uncertainty, topology, inverse modeling

This paper focuses on fault-related uncertainties in the subsurface, which can significantly affect the numerical simulation of physical processes. Our goal is to use dynamic data and process-based simulation to update structural uncertainty in a Bayesian inverse approach. We propose a stochastic fault model where the number and features of faults are made variable. In particular, this model samples uncertainties about connectivity between the faults. The stochastic 3D fault model is integrated within a stochastic inversion scheme in order to reduce uncertainties about fault characteristics and fault zone layout, by minimizing the mismatch between observed and simulated data. The stochastic fault model uses a priori information such as fault orientation, location, size and sinuosity, to sample both geometrical and topological uncertainties with realistic fault descriptions. Each fault object is parameterized by the random vector used to simulate fault features. Then, during inversion, the random vector of the current model is stochastically perturbed, producing a new parameter vector used as input by the stochastic fault model to produce a new model.

Even if the topology varies from one model to another, the algorithm produces correlated models so that their flow responses evolve quite smoothly.

The methodology is applicable in general and illustrated on a synthetic two-phase flow example. A first set of models is generated to sample the prior uncertainty space. Then, models minimizing reference water-saturation data misfit are used as seeds to generate continuous Markov chains of models with discrete states. Posterior models reduce uncertainties about fault position, while the topology varies from one model to another. A second example highlights the interest of the parameterization when interpreted data is available, by perturbing geological scenarios and falsifying those that do not match two-phase flow observations.

Introduction

In subsurface studies, uncertainties about faults and fault network primarily stem from the lack or poor quality of observation data. Due to the inherent uncertainties and possible human bias during interpretation of geological data [START_REF] Bond | What do you think this is? "conceptual uncertainty" in geoscience interpretation[END_REF], model predictions seldom match physical observations such as reservoir pressure, seismograph measurements or gravity anomalies. One may want to use such data to infer some characteristics of the model, which is known as the inverse problem.

In this paper, we propose a general method to reduce structural uncertainty about fault networks. The method generates models with varying number of faults, fault geometry and fault connectivity, i.e. how faults truncate each others. The method accounts for large to small scale faults (km to decameter scale), for which uncertainties may be large due to the lack of data (incomplete sampling or sensing), imaging ambiguities (e.g. in potential field and seismic methods), limited resolution and positioning errors (borehole trajectory, time-to-depth conversion, etc.). Therefore, in the context of hydrocarbon reservoir management, our method is applicable at several stages :

• In exploration of "green fields", in which a limited amount of data is generally available (a set of seismic lines, a few borehole data). In this case, we focus on large-scale, kilometric faults or larger.

• In "brown field" management, in which 3D seismic data is generally available but only provides approximate view on structures because of resolution, acquisition or processing issues. In this case, the method is applicable to aid the interpretation of large to mid-scale faults in poorly imaged areas or to simulate sub-seismic faults following geological priors.

In such uncertain situations, the method mainly relies on prior information about faults, which may come from the interpretation of available data, structural geological concepts, regional context or reservoir analogs.

We apply this method to the inversion of reservoir or hydrogeologic flow data. Indeed, faults impact fluid flow in different ways:

1. Individually, faults have a direct impact on the permeability field due to fracturing, alteration, cementation or shale smearing nearby the fault. Uncertainties about fault rock hydraulic properties in flow models are generally high and may have a significant impact on fluid flow [START_REF] Manzocchi | Fault transmissibility multipliers for flow simulation models[END_REF].

2. At the reservoir scale, fault zone complexity and spatial layout directly impact fluid flow by creating connected flow conduits or flow barriers isolating fault blocks, depending on the nature of the fault rock. In this paper, we focus on uncertainties about fault layout and geometry including fault connectivity uncertainties.

3. The displacement of geological layers along fault surfaces perturbs both the horizontal and vertical connectivity of transmissive layers.

3D structural models that represent geological layers and structures such as faults are of prime importance for accurately representing fault displacement, trends, anisotropy of permeability and porosity heterogeneities. The rock properties greatly depend on the nature of the rocks, thus are often estimated in a coordinate sytem representing the depositional space [START_REF] Mallet | Space-time mathematical framework for sedimentary geology[END_REF][START_REF] Gringarten | New Geologic Grids for Robust Geostatistical Modeling of Hydrocarbon Reservoir[END_REF].

The proposed methodology considers geological structures as a reduced set of parameters, thus opens the ways to the management of structural uncertainty in classical stochastic inversion methods like MCMC (figure 1).

A simple non-conditional case study (section 5.1) shows the applicability of the proposed parameterization to reduce fault location uncertainty using a rejection sampler. Indeed, the goal is to sample the posterior distribution and obtain very different realistic models, hence rejection sampling is more appropriate than optimization techniques.

In sparse data situations, several fault arrays may be interpreted, all consistent with available data. In this context, we propose to perturb geological scenarios while maintaining consistency with interpreted data, in order to falsify scenarios that do not match two-phase flow observations (section 5.2).

Related Work and Challenges of Structural Parameters Inversion

The inverse problem is known as history matching or model calibration in the case observed data is the reservoir/aquifer dynamic flow response. The idea behind inverse modeling is obvious: if a model does not reproduce the observations, it is unlikely that its predictions are correct. Moreover, if a model matches the observations, it has a higher chance to make relevant predictions. However, only stochastic approaches provide a relevant framework for risk assessment, by considering all possible scenarios and not just one deterministic calibrated model.

Most history matching techniques focus on permeability and porosity heterogeneities, either in the field of hydrogeology [START_REF] Carrera | Inverse problem in hydrogeology[END_REF] 3D structural modeling is based on the interpretation of subsurface data, which may entail uncertainties due to model approximations, sparse/low resolution data and possible interpretive bias. Once geological structures have been modeled, petrophysical heterogeneities are computed using geostatistical methods conditionning to subsurface data. Fluid-flow simulation then enables to compute flow predictions that can be compared to observations. The proposed parameterization opens the way to inversion strategies by updating not only petrophysical properties, as proposed by classical inverse methods, but also 3D geological structures.

de [START_REF] De Marsily | Dealing with spatial heterogeneity[END_REF][START_REF] Zimmerman | A comparison of seven geostatistically based inverse approaches to estimate transmissivities for modeling advective transport by groundwater flow[END_REF], or for the hydrocarbon industry [START_REF] Oliver | Recent progress on reservoir history matching: a review[END_REF]. Other techniques focus on connectivity features [START_REF] Western | Toward capturing hydrologically significant connectivity in spatial patterns[END_REF] or sedimentary objects such as channels [START_REF] Ronayne | Identifying discrete geologic structures that produce anomalous hydraulic response: An inverse modeling approach[END_REF]) that may have a great impact on fluid flows by creating high-conductivity preferred flow paths. However, uncertainties about geological structures are neglected in most cases. [START_REF] Jenni | History matching of stochastic models of field-scale fractures: Methodology and case study[END_REF] propose to handle uncertainties about fractures using gradual deformations to update field scale planar fractures while honoring a fracture density map.

In this work, we propose an approach similar to [START_REF] Jenni | History matching of stochastic models of field-scale fractures: Methodology and case study[END_REF], but use a 3D flexible and more realistic representation of fault objects in terms of geometry and topology. [START_REF] Suzuki | Dynamic data integration for structural modeling: model screening approach using a distance-based model parameterization[END_REF] focus on larger structural uncertainties. They select models that match reflection seismic data in a large set of prior structural models using a stochastic search method and the Hausdorff distance as metric to avoid running the forward problem for each ensemble member. The prior set of models is generated by perturbing the geometry of different seismic interpretations representing different geological scenarios. As [START_REF] Suzuki | Dynamic data integration for structural modeling: model screening approach using a distance-based model parameterization[END_REF], we also focus on large structural uncertainties but our parameterization is flexible enough to automatically generate models with varying number of faults and fault connections.

Other uncertainties about reservoir boundaries do not only affect reservoir volume but also control the spatial layout of geological layers and thus possibly affect fluid flows. In this context, [START_REF] Seiler | Structural Surface Uncertainty Modeling and Updating Using the Ensemble Kalman Filter[END_REF] use the Ensemble Kalman Filter (EnKF) to update top and bottom stratigraphic horizons of a reservoir. They use a corner-point grid to represent the most likely initial state and deform it to match simulated top and bottom horizons. As the grid is deformed, their algorithm keeps the same number of cells over realizations and thus a state vector with same dimensions as required by the EnKF.

However, none of these methods is flexible enough to sample posterior, i.e. after data inversion, realistic 3D models of geological structures. Geological surfaces (faults, fractures or stratigraphic horizons) are either planar or their geometry is limited by the gridding method, e.g. pillar-based grids.

Our fault representation does not suffer such limitations, for implicit modeling [START_REF] Calcagno | Geological modelling from field data and geological knowledge: Part i. modelling method coupling 3d potential-field interpolation and geological rules[END_REF][START_REF] Frank | 3D-reconstruction of complex geological interfaces from irregularly distributed and noisy point data[END_REF][START_REF] Guillen | Geological modelling from field data and geological knowledge: Part ii. modelling validation using gravity and magnetic data inversion[END_REF] is used to represent faults as isopotentials of scalar fields, namely implicit surfaces.

Inverting structural parameters using dynamic data raises both conceptual and technical issues:

1. One needs a robust method to sample structural uncertainties. While it is a common practice to generate several geostatistical realizations of facies or rock properties, few methods propose to sample uncertainties about geological structures, either geometrical [START_REF] Abrahamsen | Bayesian kriging for seismic depth conversion of a multilayer reservoir[END_REF][START_REF] Caumon | Elements for stochastic structural perturbation of stratigraphic models[END_REF][START_REF] Lecour | Modeling of stochastic faults and fault networks in a structural uncertainty study[END_REF][START_REF] Mallet | Solid earth modeling and geometric uncertainties[END_REF][START_REF] Thore | Structural uncertainties: determination, management and applications[END_REF][START_REF] Wellmann | Uncertainties have a meaning: Information entropy as a quality measure for 3-D geological models[END_REF] or topological (Cherpeau et al., 2010a;[START_REF] Holden | Stochastic structural modeling[END_REF]. The algorithm used in this work is adapted from Cherpeau et al. (2010a) in order to meet the needs of inverse modeling. It generates stochastic fault networks from prior information about fault orientation, size, sinuosity and location. 2. It is not an easy task to decribe a structural model as a set of relevant parameters. In this paper, we focus on fault network parameters, considering both individual fault parameters such as orientation, size and location and also the number of faults in the model. Connections between faults are not explicitly parameterized but result from the fault network generation process (section 3).

3. From a technical standpoint, modeling realistic geological structures requires a flexible 3D representation while modeling physical phenomena requires meshes for discretization of partial differential equations describing the physical phenomenon of interest [START_REF] Caumon | High resolution geostatistics on coarse unstructured flow grids[END_REF].

Consequently, using dynamic data to constrain geological structures calls for either back and forth and tedious conversions between the two representations, or a simplified geological representation that only enables to update the reservoir model, as proposed in the seminal work of [START_REF] Holden | Stochastic structural modeling[END_REF]. Our implicit fault surfaces can be discretized into a control volumes connectivity graph that can be fed into a flow simulator [START_REF] Cherpeau | Assessing the impact of fault connectivity uncertainty in reservoir studies using explicit discretization[END_REF]. Hence, the geological representation does not suffer from any simplifications while the flow response can be computed from the connectivity graph.

Parameterization of Fault Networks

Overview of the Stochastic Fault Model

The fault generation process is a modified version of the method proposed by Cherpeau et al. (2010a,b) to sample fault-related uncertainties in sparse data situations. The algorithm sequentially simulates faults in the domain of interest. Faults are grouped into fault families representing different faulting events in terms of timing and regional stress, resulting in different input parameters about fault orientation, size, localization and sinuosity. Moreover, fault families are simulated in a given chronological order reflecting the timing of faulting events. Consequently, new simulated faults are branching faults in case of intersections with previously simulated faults. When the timing is uncertain, fault families are considered cogenetic so that truncation is randomly simulated.

Each fault object is simulated through a number of substeps and the global simulation process is as follows:

For each fault family S, from the oldest to the youngest:

For each fault F in S:

1. Draw dip and strike from prior distributions.

2. Draw fault center coordinates from the fault center probability distribution function f built from input parameters (section 3.1.1).

3. Generate a fault geometry from orientation parameters, fault center and input sinuosity parameters.

4. Draw an extension along the strike direction and along the dip direction (section 3.1.2).

Look for intersections with previously simulated faults.

The simulation can also be conditional to interpreted fault traces. In this case, extra simulation steps assign data to faults and randomly cluster data points using geometrical concepts (Cherpeau et al., 2010b). The algorithm assigns data points to faults in different orders and generates different data combinations, hence produces different models that all honor hard data. The following sections describe some substeps of the fault object simulation.

Fault center coordinates

The fault nucleation step is a key step during the simulation process for it determines the neighborhood of a fault thus possible intersections with neighbor faults. In some geological contexts, major faults tend to repulse each other, meaning that new faults occur far from existing ones, whereas secondary faults tend to occur around major faults. If one wants to mimic such behaviors, the fault center probability distribution function f may be non-uniform in the domain of interest and may depend on neighbor faults.

Consequently, the fault center is drawn from f (u, v, w) defined in a Cartesian grid, by successive drawing of

u ∈ [1; n u ], v ∈ [1; n v ] and w ∈ [1; n w ]
coordinates, with n u , n v , n w the number of voxels along the three orthogonal axis a u , a v and a w respectively. The probability p(U = a) that the u coordinate is equal to a is:

p(U = a) = ∫ v ∫ w f (a, v, w)dvdw = nv ∑ j=1 nw ∑ k=1 f (a, j, k) (1) Once u coordinate has been drawn, p(V = b|U = a) with b ∈ [1; n v ] is computed according to: p(V = b|U = a) = ∫ w kf (a, b, w)dw = nw ∑ k=1 kf (a, b, k) with k a normalization constant, k = 1 ∫ v ∫ w f (a, v, w)dvdw = 1 nv ∑ j=1 nw ∑ k=1 f (a, j, k) (2)
Then, w can be drawn in the same manner, i.e. by computing a discrete pdf along the axis a w given u = a and v = b coordinates. Continuous values of these discrete probability density functions are obtained by linear interpolation between discrete values, that can then be back-transformed into x, y, z coordinates.

Fault extension

In order to model laterally-terminating and synsedimentary faults, the algorithm restricts the existence of a fault by a 3D ellipsoid E(x, v 1 , v 2 , v 3 ) with x the fault center coordinates, v 1 , v 2 , v 3 three orthogonal vectors. The vectors v 1 and v 2 are respectively computed from extension values randomly drawn from length distribution along the strike direction and height distribution along the dip direction. The third vector v 3 is fully defined by

v 3 = v 1 × v 2 .
This elliptical fault shape is not a restriction of our algorithm and is consistent with the litterature. Indeed, fault's tipline is considered to be circular to elliptical for an isolated fault [START_REF] Barnett | Displacement geometry in the volume containing a single normal fault[END_REF][START_REF] Walsh | Formation of segmented normal faults: a 3-D perspective[END_REF] but may evolve to elongated ellipses if the vertical fault height is restricted by the thickness of a geological layer (Benedicto et al., 2003), which can be approximated by choosing ||v 1 || >> ||v 2 ||.

Fault number parameter

In case of large fault-related uncertainties, not only is the fault connectivity uncertain but also the number of faults occuring in the studied area. The algorithm associates, for each fault family, a theorical number of faults n with a probability of fault existence p e ∈ [0, 1], so that the output number of faults m varies in [0, n] given the following probability density function

Q(m): Q(m) = ( n m ) p m e × (1 -p e ) n-m (3)
The following section presents the parameterization of our stochastic fault model to fit the needs of an inversion scheme.

Model Parameterization

The

Inversion Methodology

The parameterization described in section 3.2 opens ways to inversion schemes in which geological structures are part of the inversion loop, as shown in fig-

ure 1. In this section, we present how the parameterization can be used to sample fault networks using two-phase flow response as a misfit to minimize. The inversion scheme relies on a version of the Metropolis algorithm [START_REF] Metropolis | Equation of State Calculations by Fast Computing Machines[END_REF] proposed by [START_REF] Mosegaard | Monte carlo sampling of solutions to inverse problems[END_REF]. The problems can be found in [START_REF] Mosegaard | Monte Carlo analysis of inverse problems[END_REF].

Random Walk

Inverse methods search for solutions that minimize the misfit between predictions and observations and thus need random walks or perturbation processes to explore the parameter space. In this space, the neighborhood of point i is defined as all the points j that can be reached from i in a single move. We define the neighborhood N i of a model m i (r i 1 , ..., r i 8n ) in 8n-D search space as the set of all models m j (r j 1 , ..., r j 8n ) that can be obtained through the following steps:

For each random number r i k :

1

m i to model m j : p(m i |m j ) = p(m j |m i ).
For a point i in parameter space, the perturbation amplitude p defines the size of the neighborhood around i, thus controls the speed of exploration of the parameter space. Moreover, p reflects the trade-off between exploitation and exploration of the parameter space. Indeed, if p is small, it means that the algorithm searches for models in a close neighborhood and thus tries to exploit current position in parameter space. If p is large, the neighborhood is large and thus the algorithms tends to explore the parameter space. However, high-dimensional spaces tend to be empty spaces and the probability of moving from a high-misfit point i to a low-misfit point j may be high if the move, i.e. the perturbation amplitude p, is too large.

Parameter variation speed

Another way of generating a new model could be to directly perturb the fault parameters of the current model, i.e. to bypass the statistical distributions. When perturbing random numbers, all numbers evolve at the same "speed" in the range [0, 1] depending on the input perturbation amplitude p. Then, each fault parameter is drawn from the corresponding distribution, meaning that fault parameters evolve at a speed proportional to the spread of their distribution. Indeed, the narrower the distribution, the slower the parameter variation speed. Consequently, the method accounts for the prior information and uncertainty. This could be hardly achieved by perturbing directly the fault parameters.

Structural parameters evolution

The parameterization and perturbation process described above ensures that fault parameters evolve smoothly. However, the fault sinuosity is described by a single random number used as a seed for generating a correlated random noise using a geostatistical algorithm, e.g. Sequential Gaussian Simulation.

The method is similar to [START_REF] Caumon | Elements for stochastic structural perturbation of stratigraphic models[END_REF], the generated random field is added to the scalar field representing the fault surface. The input variogram ranges needed for the geostatistical method is computed from sinuosity parameters. In contrast with gradual deformation methods, the fault sinuosity does not evolve smoothly since perturbing the seed entails the generation of a new perturbation field.

Moreover, the fault connectivity is not parameterized as it results from the stochastic fault simulation algorithm. However, if two faults are no longer connected, the fault positions evolve smoothly so that the flow has a small volume of matrix to pass trough to connect the two disconnected faults. Consequently, even if some fault connections are lost from one model to another, the corresponding flow responses are still correlated.

Fault number evolution

Adding or removing a fault from one model to another may entail a large perturbation of the model predictions, e.g. if the new or removed fault is localized in a flow-influencing zone, the flow response change may be too large to consider correlated models. In order to get smooth transitions, the existence of a fault for a given fault family is coupled with its length distribution so that the probability of a zero-length fault is 1-p e with p e the probability of existence. Consequently, the length distribution is modified as illustrated in figure 2. The ϵ parameter should be chosen carefully in function of the perturbation amplitude p, for it controls the fault growth/contraction rate. Indeed, if p is too large as compared to ϵ, the fault may get its normal size or disappear during one iteration.

Acceptance Rule

The acceptance rule modifies the random walk (presented in section 4.1) by introducing a probability p ij to move from point i to point j, so that the posterior probability density function is sampled. We use the Metropolis observations w obs t :

S(m) = l ∑ t=1 (w t (m) -w obs t ) 2 l (4)
The likelihood function is of the form L(m) = exp(-S(m) s 2 ), with s 2 the total noise variance. Then, the probability p ij of moving from point i to point j is:

p ij = 1 if L(m j ) ≥ L(m i ) = L(m j) L(m i ) otherwise = exp(- S(m j ) -S(m i ) s 2 ) (5)
Note that the misfit function could be any function of interest for the study at hand. We here use it for the two examples presented in section 5.

Moreover, the methodology could be used in other scientific areas. It could for instance be used in geophysical inversions to better constrain subsurface structure. In this case, a wave propagation method would replace the flow simulation during the forward modeling step in the inverse methodology.

Applications

Updating Prior Structural Uncertainty

This case study aims at understanding the behavior of the method for a simple case and assessing the amount of structural information introduced by a simple history matching problem. Fault occurence and location are not constrained by any interpreted data.

Input Data

The studied area of size 1400x1400x200m has one injector and one producer wells on both sides of a fault zone composed of two fault families:

1. The first fault family corresponds to major faults aligned along the main flow path:

• strike: uniform, [110 -130] • • dip: Gaussian, mean = 70 • , σ = 3 • • length: Gaussian, mean = 800m, σ = 50m
• height: Gaussian, mean = 350m, σ = 5m

• number of faults: 3

• probability of existence p e = 0.8

• sinuosity wavelength along strike direction: 500m

• sinuosity wavelength along dip direction: 500m

2. The second fault family corresponds to branching faults:

• strike: uniform, [50 -70] • • dip: Gaussian, mean = 70 • , σ = 3 • • length: Gaussian, mean = 500m, σ = 30m
• height: Gaussian, mean = 300m, σ = 5m

• number of faults: 2

• probability of existence p e = 0.7

• sinuosity wavelength along strike direction: 350m

• sinuosity wavelength along dip direction: 300m

The total number of faults varies in [0,5] according to equation 3. The matrix has a constant porosity of 0.2 and a constant permeability of 100mD.

Faults are considered flow barriers with a permeability set at 0.01mD, hence uncertainty about fault transmissibility is not investigated.

Inversion Scheme

The Metropolis algorithm is designed to sample asymptotically the posterior distribution function so that the starting point or seed to generate the Markov chain of models is meaningless [START_REF] Mosegaard | Monte carlo sampling of solutions to inverse problems[END_REF].

Then, after the burn-in time is over, samples of the posterior can be picked according to a "waiting time" to get independent samples.

To speed up the process, we propose to generate first 2000 samples from the prior and then to use the 100 models that best match flow data as seeds for generating 100 independent Markov chains of models using the acceptance rule described in section 4.2. Each chain stops when the misfit is below a given threshold or 500 models are generated. This process does not guarentee unbiased sampling of the posterior distribution, but corresponds to a practical use of the method in which only models deemed acceptable are retained for further use.

N 1400m

Figure 3: The reference model has 4 faults. The shaded zone represents the fault zone thus ensures faults to be inside this zone.

Reference Model and selected models

The reference model is randomly drawn among the 2000 prior models and is presented in figure 3. Then, the misfit is computed for other prior models using equation 4 and water-cut data. Results are presented in figure 4.

The misfit threshold to accept posterior models is set to 10 (from 51 for the worst of the 100 selected models); only four models are below that threshold among the 2000 models.

Sensitivity to Perturbation Amplitude Parameter

The perturbation amplitude parameter p controls the exploitation versus exploration of the parameter space. To test the sensitivity of our methodology to parameter p, the best model is perturbed using different values of p.

100 models are generated for each value of p, and their flow-response misfit is computed and summarized in figure 5. Statistics are almost constant for small perturbations ranging from 0.0001 to 0.03 maybe due to the impact of the fault sinuosity seed that generates a new sinuosity whatever the perturbation value. Moreover, maximum values are unexpectedly constant, maybe due to a fault configuration in unstable state, e.g. a fault center coordinate close to a neighbor fault may switch of side even for small changes. Minimum misfit values are also almost constant due to the perturbation process that does not exclude small perturbations even for large p values. Consequently, we focus on centile values that do not depend on extrem values.

For p ranging from 0.04 to 1, the median misfit increases linearly with a correlation coefficient of 0.997.

These results greatly depend on prior information, thus cannot be generalized. Indeed, as described in section 4.1.1, the parameters evolution speed depends on the parameters distribution, e.g. a Dirac distribution always generates the same value thus the speed is null whatever the perturbation amplitude is. Given these results, the perturbation amplitude p is set to 0.03 for the case study which ensures small misfit changes while enabling exploration of the parameter space.

Results

Any perturbation among the fault parameters may entail a fault network topology change and impact the flow response. Moreover, parameters are connected, e.g. a change of fault position increasing the distance between two faults may be balanced by a fault growth decreasing the distance between the two faults, thus cannot be considered individually. In case the fault size is small as compared to other faults, changes of fault location or orientation may entail unsignificant changes in flow response. Consequently, we propose to compare fault occurence probability maps (figures 6 and 7) instead of posterior distributions.

Only three Markov chains among the 100 chains did not manage to reduce the misfit under the input threshold. An example of misfit evolution for a seed model is shown in figure 8. The average acceptance rate is about 43% for all Markov chains, which is relevant with prior studies [START_REF] Gelman | Efficient Metropolis jumping rules[END_REF][START_REF] Mosegaard | Monte Carlo analysis of inverse problems[END_REF]. Both prior probability maps do not show any particular feature but a large zone with medium probability. Posterior probability maps show specific features that are present in the two subsets which suggests relative convergence and stability of results.

Posterior probability maps relatively well resolve fault location uncertainty. A typical realization would have two faults of family A on both sides of the main flow path and one fault of family B playing the role of flow barrier along the flow path, which is relatively similar to the reference model.

Other faults located elsewhere in the model do not impact the flow response, which may explain why these areas are not very well resolved on the posterior probability maps. Indeed, any hard data constrains fault location and the impact of conditioning is still to be evaluated.

Analysis of Structural Scenarios

The interpretation of subsurface data may lead to different interpretations [START_REF] Bond | What do you think this is? "conceptual uncertainty" in geoscience interpretation[END_REF] and thus different geological scenarios. Fluid-flow data represent additional data that may help to analyze different interpretations and to falsify those that do not match observations [START_REF] Tarantola | Popper, bayes and the inverse problem[END_REF]. In this Only models that match flow data are retained and other interpretations are falsified. The study is not meant to be statistical but practical, hence scenarios that do not have any neighbor model with an acceptable misfit to water-cut data after a large number of iterations are falsified.

Data Overview

Eight structural scenarios have been built from four interpreted cross-sections, one of them is considered as the reference scenario (figure 9). There are one injector in the center and four producers in the corners of the studied area of size 1400x1400x200m. As in section 5.1, the misfit is the root mean squared error (equation 4) of water-cut data. Matrix porosity is set at 0.2, matrix permeablity is 5mD and fault permeability is 0.01mD.

Results

Initial interpretations are perturbed until an acceptable misfit is reached.

This arbritary choice is motivated by a practical need to obtain a match in limited time, a more rigorous version of MCMC could be run in principle to better sample the search space, at the expense of computational efficiency.

The evolution of the misfit from one model to another is presented in fig-

ure 10. The average acceptance rate for all scenarios is about 33%. Scenarios The smooth evolution of fault parameters, even for the number of faults coupled with size evolution, enables to generate correlated flow responses from one sample to another. Only the fault sinuosity does not evolve smoothly from one sample to another but gradual deformation could help to obtain more continuous transitions [START_REF] Hu | Gradual deformation and iterative calibration of gaussian-related stochastic models[END_REF]. Other fault-related uncertainties could be considered to complete the method, such as fault trans-missibility which may have a significant impact on the flow response. In the presented examples, faults occur in a homogeneous geological layer but rock heterogenetities could also be taken into account. In such cases, fault displacement should also be modeled for it may connect or disconnect transmissive layers. This point will be addressed in a future contribution.

In our first simple MCMC example, the forward problem, i. [START_REF] Sambridge | Geophysical inversion with a neighbourhood algorithm-I: Searching a parameter space[END_REF]Scheidt and Caers, 2009;[START_REF] Suzuki | Dynamic data integration for structural modeling: model screening approach using a distance-based model parameterization[END_REF] could also be used either separately or in conjunction with the proposed method for more efficient sampling.

Interpretations of inversion results are difficult to analyze due to the amount of data and correlation between parameters. As argued by [START_REF] Mosegaard | Monte carlo sampling of solutions to inverse problems[END_REF], there is still a need for methods to summarize information and re-

Figure 1 :

 1 Figure1: Typical modeling and inversion workflow. 3D structural modeling is based on the interpretation of subsurface data, which may entail uncertainties due to model approximations, sparse/low resolution data and possible interpretive bias. Once geological structures have been modeled, petrophysical heterogeneities are computed using geostatistical methods conditionning to subsurface data. Fluid-flow simulation then enables to compute flow predictions that can be compared to observations. The proposed parameterization opens the way to inversion strategies by updating not only petrophysical properties, as proposed by classical inverse methods, but also 3D geological structures.

  stochastic fault model sequentially simulates faults in the domain of interest. Each step during the simulation of a fault object stochastically set a fault parameter. Finally, each fault can be described by six parameters: dip, strike, length, height, sinuosity and center. During the fault simulation, one or more random numbers are used by the stochastic steps to generate each fault parameter, i.e. only one random number is used per dip, strike, height, length parameters whereas three random numbers are required to draw the fault center position. Consequently, each fault can be parameterized by a vector with eight components F (r 1 , ..., r 8 ). A model m with n faults is fully parameterized by a random vector with 8n components m = {F 1 , ..., F n } or m = {(r 1 1 , ..., r 1 8 ); ...; (r n 1 , ..., r n 8 )}.

  Metropolis algorithm was originally designed to sample Gibbs-Boltzmann distributions in chemical science and has been since successfully applied in several fields such as mathematic, physic and biology. It is a Markov Chain Monte Carlo method, Markov Chain implicitly meaning first-order Markov chain, i.e. current state only depends on previous state, and Monte Carlo meaning that it is a random sampling method. The Metropolis algorithm is an importance sampling algorithm designed to solve highly nonlinear inverse problems by sampling the posterior probability density function defined by the product of the likelihood function and the prior probability density function. A general description of Monte Carlo methods as solutions to inverse

  . Compute b low = max(0, (r i k (1 -p))) with p a perturbation amplitude given in input; 2. Compute b up = min(1, (r i k (1 + p))); 3. Define a uniform distribution D k U[b low , b up ]; 4. Randomly draw r j k from D k .Consequently, the neighborhood N i corresponds to an infinite set of models.However, any point B in parameter space can be reached from any point A provided a sufficient number of steps, i.e. the random walk covers the whole parameter space. Moreover, the proposed transition process from model i to model j is symmetric, i.e. the probability p(m i |m j ) of moving from model m j to model m i is equal to the probability p(m j |m i ) of moving from model

Figure 2 :

 2 Figure 2: Fault length distribution coupled with existence probability parameter. At the top, an initial gaussian CDF representing the fault length with mean m. At the bottom, the distribution has been modified to account for the existence probability p e . Parameter ϵ should be carefully set in function of the perturbation amplitude p.

Figure 4 :

 4 Figure 4: Mifsit repartition computed from equation 4. Only the 100 selected models are used as seeds to generate Markov chain of models and four of them have a misfit below 10 thus are considered as posterior models.

  section, we propose to analyze different structural scenarios built from the interpretation of seismic sections by generating a family of models around each interpretation. The sinuosity and size of each deterministic fault of the initial interpretation are perturbed. Other fault parameters such as fault center location and orientation are frozen. The initial interpretations are used as seeds to generate Markov chains of models as in section 5.1 and conditioning to seismic interpretations is made during neighborood search.

I

  1 , I 4 , I 5 and their corresponding chain of models have a flow-response that remains far from the reference thus are falsified. Scenarios I 0 , I 2 , I 3 , I 6 have a neighbor model within an acceptable misfit range, respectively I ′ 0 , I ′ 2 , I ′ 3 , I ′ 6 (figure 11), which are considered consistent with available water-cut data. Conclusions would have been much different by only considering the initial interpretations misfit, e.g. interpretations I 2 , I 3 , I 6 have a larger initial misfit than falsified interpretations I 1 , I 4 , I 5 . The low-misfit models I ′ 0 , I ′ 2 , I ′ 3 , I ′ 6 can then be used for further modeling 6 Discussion and Conclusions The proposed parameterization of faults enables to sample posterior fault networks with various number of faults and fault connections. The modeling method uses a realistic 3D representation of fault surfaces that goes beyond planar elliptical surfaces and enables to model realistic fault arrays. It considers faults as a random vector corresponding to fault parameters and opens the way to inversion of geological structures using realistic descriptions. It goes beyond traditional inversion schemes by pertubing both fault geometry and fault network topology. The first presented example shows how fault position uncertainty can be reduced using our parameterization and a version of the Metropolis algorithm to sample posterior models. The sampled 100 posterior models suggest specific features about fault location. The second example shows the interest of the methodology when conditioning data is available, by analyzing structural scenarios and falsifying those that do not have any neighbor model consistent with observations. The simulation method enables to permanently honor data points during the perturbation process, which reduces the neighborhood around selected scenarios and ensures the consistency of the proposed models.

  e. the flow simulation, represents about 91% of the computation time (about 9% for fault simulation, the time for discretization and generation of flow input files is negligible) for an average of 1min20s per MCMC iteration on a 3.4Ghz PC with about 8000 flow elements and 16000 connections per discretized model. When complexity increases as in the second example with a larger number of faults and wells, this computation time is even more important (about four times longer than first example) even if time proportions remain relatively constant: the inversion spends 88% of the computation time to solve flow equations, 12% to simulate faults (maybe due to extra simulation steps needed for data conditioning, other steps are insignificant). These performance issues and the large number of parameters do not enable to rigorously sample the posterior. The presented inversion strategy is parallel by nature since Markov chains of models can be computed independently, which reduces the total computation time of the presented examples. Other methodologies could also be used in order to reduce the computation time, e.g. Mariethoz et al. (2010) propose a stochastic stopping criterion to reduce the number of forward runs. Distance-based optimization

Figure 5 :Figure 9 :Figure 10 :Figure 11 :

 591011 Figure 5: Misfit values in function of perturbation amplitude p in logarithmic scale. Vertical bars represent interquartile ranges with median value. For p values ranging from 0.0001 to 0.03, median values are almost constant, then increase with larger interquartile ranges. Minimum (red squares) and maximum (green squares) values are almost constant for all p values due to the perturbation process and maybe unstable fault position.

Acknowledgments

The authors would like to acknowledge the Gocad Consortium for funding this research as well as Stanford University for providing the GPRS flow simulator. This work is partially supported by the European Research Council (ERC-StG-205693).