Paul Baville
email: paul.baville@univ-lorraine.fr

Marcus Apel

Silvan Hoth

Dirk Knaust

Christophe Antoine

Cédric Carpentier

Guillaume Caumon

Computer-assisted stochastic multi-well correlation: Sedimentary facies versus well distality

Keywords: Stratigraphy, Well Correlations, Sedimentary Facies, Well Distality, Walther's Law of Facies, Dynamic Time Warping

To reflect interpretation uncertainties, an increasing number of interpreters may generate different sets of correlations from the same data set. To reduce human interpretation bias, we use a computer-assisted method for generating stochastic multi-well correlations to be used in uncertainty studies. This method relies on a translation of sedimentary concepts into cost functions that are related to each possible correlation. All these correlation costs are used to populate a cost matrix in order to apply the Dynamic Time Warping algorithm and to compute the n-best correlation sets having the n-least cumulative costs. The proposed cost function reflects prior knowledge about sediment transport direction, and it is tested on two wells from the North Sea, which penetrate a Middle Jurassic reservoir interval. Well markers are described by two parameters: (1) the sedimentary facies representing the depositional environment, and (2) the relative distality of the well computed from its position along the sediment transport direction. The main principle of correlation used in this article follows Walther's Law of Facies: a well marker (described by facies and distality) cannot be correlated with another well marker described by a depositionally deeper facies at a more proximal position, or a depositionally shallower facies at a more distal position. This approach produces consistent stratigraphic correlations, and highlights the sensitivity of the solution to the facies zonation and to the relative disability. Therefore, the proposed rule offers a way to coherently consider chronostratigraphic correlation and the associated uncertainties at a scale smaller than generally considered in deterministic correlation.

Introduction

The stratigraphic correlation of sedimentary rock series observed in borehole data or outcrops is a common procedure in basin studies. Correlation affects the understanding of sedimentary deposits, the quantitative analysis of subsurface resources and source-to-sink models. In order to predict the subsurface geometry, stratigraphy is used to correlate sedimentary horizons from sparse wells (e.g., [START_REF] Nagy | Stratigraphy[END_REF]. These correlations may be done at a large scale from reflection seismic data, or at a smaller scale from well data. In both cases, the main objective of the stratigraphic correlation is to determine the surfaces bounding the stratigraphic units, where geometry can then be predicted away from observations. In this work, we focus on stratigraphic correlation generated from well data (small vertical meter-scale and large horizontal kilometer-scale). There are two types of well data: (1) well logs, corresponding to quasi-continuous (regular sampling) geophysical measurements along the well path (e.g., gamma ray, sonic, neutron porosity), and (2) well zones, corresponding to categorical reservoir properties and defined by their top and bottom depths along the well path (e.g., biozones, coarsening-up intervals, structural zones, sedimentary facies).

Well correlation may be achieved manually by experts using well log and core sample interpretations (e.g., [START_REF] Bourquin | High-resolution sequence stratigraphy of an alluvial fan-fan delta environment: Stratigraphic and geodynamic implications -An example from the Keuper Chaunoy Sandstones, Paris Basin[END_REF]. The resulting scenario is a correlation set. This process generally involves iteration between the identification of intervals or markers on each well, and the correlation of these features across wells. Although uncertainties in correlation exist, a geologist generally proposes the "best case" scenario based on a careful data analysis and in accordance to their experience, or sometimes more (two to five) scenarios, reflecting some -but not all -interpretation uncertainties. Increasing the number of interpreters may help to gain confidence about a particular interpretation, or conversely to produce more scenarios.

Ideally, a quantification of correlation uncertainties would call for generating every possible scenario and to evaluate their likelihood. However, the very large number of possible scenarios makes their manual validation almost impossible. To overcome this limitation, a numerical method to provide automatically multiple stratigraphic scenarios is proposed (e.g., [START_REF] Pels | Automated biostratigraphic correlation of palynological records on the basis of shapes of pollen curves and evaluation of next-best solutions[END_REF][START_REF] Lallier | Uncertainty assessment in the stratigraphic well correlation of a carbonate ramp: Method and application to the Beausset Basin, SE France[END_REF][START_REF] Edwards | Uncertainty management in stratigraphic well correlation and stratigraphic architectures: A training-based method[END_REF]. The scenarios are consistent with the input data and with some explicitly defined stratigraphic concepts, which guarantee that the solutions are acceptable.

In this paper, we consider chronostratigraphic principles (Section 1), which are considered as the most predictive to model the subsurface. Among chronostratigraphic approaches, sequence stratigraphy describes the stratigraphic structure of sedimentary basins based on depositional sequences (e.g., [START_REF] Mitchum | Seismic stratigraphy and global changes of sea level: Part 6. Stratigraphic interpretation of seismic reflection patterns in depositional sequences: Section 2. Application of seismic reflection configuration to stratigraphic interpretation[END_REF][START_REF] Galloway | Genetic stratigraphic sequences in basin analysis I: Architecture and genesis of flooding-surface bounded depositional units[END_REF][START_REF] Catuneanu | Sequence stratigraphy: Methodology and nomenclature[END_REF]. Depending on the well distality (well position along a proximal-to-distal transect), depositional processes vary and lead to different isochronous sedimentary facies and vertical stacking patterns (e.g., [START_REF] Borgomano | Stratigraphic well correlations for 3-D static modeling of carbonate reservoirs[END_REF][START_REF] Bhattacharya | Practical problems in the application of the sequence stratigraphic method and key surfaces: Integrating observations from ancient fluvial-deltaic wedges with Quaternary and modelling studies[END_REF][START_REF] Lallier | Uncertainty assessment in the stratigraphic well correlation of a carbonate ramp: Method and application to the Beausset Basin, SE France[END_REF][START_REF] Knaust | Bay-head deltas as hydrocarbon reservoirs: The Middle Jurassic Hugin Formation in block 15/3 of the South Viking Graben, Norway[END_REF].

To retrieve such isochronous surfaces from stratigraphic series, we rely on the Dynamic Time Warping algorithm (Section 2). This numerical tool computes the likelihood of a correlation set among all possible correlation sets and returns the best correlation set, i.e., the most likely correlation set [START_REF] Sakoe | Dynamic programming algorithm optimization for spoken word recognition[END_REF][START_REF] Waterman | The match game: New stratigraphic correlation algorithms[END_REF][START_REF] Hale | Dynamic warping of seismic images[END_REF][START_REF] Wheeler | Simultaneous correlation of multiple well logs[END_REF] or the n-best correlation sets [START_REF] Pels | Automated biostratigraphic correlation of palynological records on the basis of shapes of pollen curves and evaluation of next-best solutions[END_REF][START_REF] Lallier | Uncertainty assessment in the stratigraphic well correlation of a carbonate ramp: Method and application to the Beausset Basin, SE France[END_REF][START_REF] Edwards | Uncertainty management in stratigraphic well correlation and stratigraphic architectures: A training-based method[END_REF]. The likelihood of one correlation set is computed by minimizing the cumulative cost of correlations and transition between correlations. In this paper, the key methodological contribution is the explicit integration of distality in the elementary correlation costs of the Dynamic Time Warping. This enables to compute global correlation likelihood and to rule out impossible lateral facies variations along the generated isochronous layers.

In this method, all input data (well logs and core samples of the interval of interest) must be acquired within unfaulted stratigraphic parasequence successions. This consideration may be justified in the presence of seismic data or dipmeter data analysis. We focus on relative chronostratigraphic well marker correlation based on (1) the sedimentary facies, that represents the depositional environment interpreted as zones along the well paths, and (2) the relative distality of the well along the sediment transport direction. The proposed method is first validated with a simple synthetic data set (Section 2.3) and secondly applied on an unfaulted (validated by dipmeter data) sedimentary coastal deltaic data set from the South Viking Graben in the North Sea [START_REF] Knaust | Bay-head deltas as hydrocarbon reservoirs: The Middle Jurassic Hugin Formation in block 15/3 of the South Viking Graben, Norway[END_REF] (Section 3).

Conceptual stratigraphic correlation framework

Lithostratigraphy versus chronostratigraphy

Considering the empirical principle of superposition stating that "strata have been deposited above older strata than themselves" [START_REF] Steno | De solido intra solidum naturaliter contento dissertationis prodromus[END_REF], stratigraphy can be interpreted in, at least, two different ways (Figure 1) and may lead to very distinct correlations with direct implication for paleogeography and reservoir connectivity:

-Lithostratigraphy correlates strata by comparing their composition or rock type. Well markers are correlated if they are similar (or, even better, identical) (Figure 1.A) to generate layers or beds containing the same facies (characteristics) (Figure 1.B). However, this method does not consider possible lateral shifts of facies along a proximal-to-distal direction in the basin. Sedimentary reservoir correlation is usually based on deterministic well-data analysis. Correlating these data is essential in subsurface modeling prior to defining the geometry of geobodies and to model their petrophysical properties. In quantitative subsurface modeling and sedimentological studies, chronostratigraphic time line interpretation is generally preferred to lithostratigraphic principles because it best reproduces the lateral and vertical facies variations derived from depositional processes and the continuity of thin beds (Figure 1) (e.g., [START_REF] Ainsworth | Correlation techniques, perforation strategies, and recovery factors: an integrated 3-D reservoir modeling study, Sirikit field, Thailand[END_REF]. However, lithostratigraphy is often used to complement chronostratigraphy when observations are closer than the lateral size of sedimentary bodies and the chronostratigraphic resolution.

Sequence stratigraphy

Sequence stratigraphy provides several concepts and methods to analyze, characterize and predict stratigraphic and sedimentary architecture in the subsurface. The main idea is to subdivide the sedimentary strata into several stratigraphic sequences based on genetic considerations about the depositional history [START_REF] Sloss | Integrated facies analysis, in: Sedimentary facies in geologic history[END_REF]. A sequence is defined as "a relatively conformable succession of genetically related strata bounded by unconformities or their correlative conformities" [START_REF] Mitchum | Seismic stratigraphy and global changes of sea level: Part 11. Glossary of terms used in seismic stratigraphy: Section 2. Application of seismic reflection configuration to stratigraphic interpretation[END_REF]. The key idea of sequence stratigraphy is to relate stratigraphic units to cycles of accommodation and sedimentation rate, which are controlled by global (eustatic) seal-level changes, local paleotopography, subsidence, and climate [START_REF] Johnson | Time-rock model for Siluro-Devonian continental shelf, Western United States[END_REF][START_REF] Galloway | Genetic stratigraphic sequences in basin analysis I: Architecture and genesis of flooding-surface bounded depositional units[END_REF][START_REF] Ruiz | Sequence stratigraphy and facies analysis -Computer aided interpretation and correlation[END_REF][START_REF] Catuneanu | Sequence stratigraphy: Methodology and nomenclature[END_REF]. This concept leads to several depositional system tracts subdividing these sequences [START_REF] Catuneanu | Sequence stratigraphy of clastic systems: Concepts, merits, and pitfalls[END_REF]:

-The Transgressive system tract develops during a period of positive accommodation a, which is not compensated by sedimentation rate s (a/s > 1). The vertical stacking of facies is fining-upward for both marine and alluvial environments, coeval with a backstepping of coastal marine facies in a landward direction. The system tract is bounded at the top by a maximum flooding surface.

-The Regressive system tract includes all regressive deposits corresponding to two successive system tracts:

1. The Highstand system tract is formed when the sedimentation rate becomes larger than the accommodation (normal regression, 0 < a/s < 1) [START_REF] Catuneanu | Towards the standardization of sequence stratigraphy[END_REF]. A part of the sediment volume is stacked vertically (aggradation), but the excess material is deposited ahead of the depositional system (progradation). This system tract starts just above the maximum flooding surface and includes, in marine settings, coarsening-upward (prograding) deposits.

2. The Lowstand system tract is formed after a period of negative accommodation (forced regression, a < 0) at the beginning of a renewal in accommodation, when the sedimentation rate is larger than the rate of accomodation (normal regression, 0 < a/s < 1). Its upper limit is the maximum regressive surface, which can merge with the transgression surface [START_REF] Embry | Practical sequence stratigraphy[END_REF]. In marine settings, it is characterized by coarsening-upward deposits, which are not necessarily entirely preserved due to ravinement erosion during transgression.

These different system tracts correspond to successions of regressive and transgressive periods recorded in sediments (Figure 2.A). [START_REF] Wheeler | Time-Stratigraphy[END_REF] proposed combined diagrams, which show the same stratigraphic section using two different representations: (1) a stratigraphic section as function of depth, and (2) a stratigraphic section as function of geological age (Figure 2.B) [START_REF] Driscoll | Three-dimensional quantitative modeling of clinoform development[END_REF][START_REF] Mallet | Space-Time Mathematical Framework for Sedimentary Geology[END_REF]Holbrook and Bhattacharya, Figure 2: Successive transgressive and regressive parasequences corresponding respectively to transgressive and regressive system tracts. Inspired by [START_REF] Kedzierski | 3D Marine Sedimentary Reservoir Stochastic Simulation Accounting for High Resolution Sequence Stratigraphy and Sedimentological Rules[END_REF]. (A) Depth cross-section and (B) time cross-section (Wheeler diagram;[START_REF] Wheeler | Time-Stratigraphy[END_REF]). (single column -colored)

2012; [START_REF] Amosu | WheelerLab: An interactive program for sequence stratigraphic analysis of seismic sections, outcrops and well sections and the generation of chronostratigraphic sections and dynamic chronostratigraphic sections[END_REF][START_REF] Qayyum | The Wheeler diagram, flattening theory, and time[END_REF].

Sedimentary facies description helps in the interpretation of the paleoenvironment and depositional event successions along the well path. While the herein proposed method is in general applicable to various sedimentary environments, it has been applied to marginal-marine deltaic wave-dominated deposits.

Principle of correlation

In this work, relative chronostratigraphy is used to associate well markers. The objective of this work is to correlate time windows. We can assume that over a given area, some deposition events are quasi-synchronous depending on basin geometry, subsidence, distribution and depositional system. Therefore, along a given isochronous surface, all deposits have (approximately) the same geological age (Figure 2.B). In this case, correlating these isochronous deposits is the same as correlating time windows and vice versa.

A sedimentary domain can be characterized by two main directions: proximalto-distal direction (dip), given by the regional sediment transport direction, and the direction perpendicular to it (depositional strike). In practice, these two directions may change laterally (e.g., in fan-shaped deltas), but, for simplicity the case of a Cartesian geometry is considered here to schematize the internal structure of a sedimentary basin margin (Figure 3). In the direction of sediment transport (BB'), there is a superposition of sigmoidal bodies whereas along the depositional strike direction (AA'), a superposition of monoclinal bodies occurs. These sedimentary bodies are defined as clinoforms [START_REF] Rich | Three critical environments of deposition, and criteria for recognition of rocks deposited in each of them[END_REF] and are sets of facies (from the most proximal to the most distal region) during the same depositional event. Depending on the context (distance to fluvial sediment source, dominance of waves / tides in the depositional process), clinoforms are observed in every direction, even though they have different shapes along depositional strike or sediment transport direction.

Figure 3 illustrates two successive clinoforms formed during a normal regression and bounded by isochronous surfaces (blue lines). These two clinoforms may be divided into four sedimentary facies classified by their relative depositional depth. These clinoforms may be defined at a large time interval (4th orderparasequences), in the range of 150 000 to 200 000 years [START_REF] Knaust | Bay-head deltas as hydrocarbon reservoirs: The Middle Jurassic Hugin Formation in block 15/3 of the South Viking Graben, Norway[END_REF]. And at a lower scale, clinoforms can be subdivided in several stratigraphic system tracts.

A clinoform may be subdivided into several sequences and the purpose of this work is to correlate well markers in order to reproduce these subdivisions and to generate a chronostratigraphic model of the subsurface.

Figure 4 illustrates that along the two principal directions, well marker correlation likelihood may be qualitatively evaluated by applying Walther's Law of Facies. This leads to several configurations, which are not plausible if we assume that the depositional depth increases with distality:

-A sedimentary facies in a well cannot be associated with a depositionally shallower facies in a more distal well.

-A sedimentary facies in a well cannot be associated with a depositionally deeper facies in a more proximal well.

Considering these two conditions, it is possible to postulate a principle of correlation stating that a well marker (described by facies and distality taken at the center of an interval having a constant facies and a constant distality) cannot be correlated with another well marker described by a depositionally deeper facies at a more proximal position, or a depositionally shallower facies at a more distal position.

This principle allows us to evaluate every single well marker correlation line independently and not the entire correlation set. Therefore, the next section focuses on the computation of likely correlation sets.

Correlation likelihood computation -Dynamic Time Warping

The Dynamic Time Warping algorithm is a dynamic programming method which can be used to quantify the likelihood of a large number of scenarios by finding the best one [START_REF] Sakoe | Dynamic programming algorithm optimization for spoken word recognition[END_REF]. In our method, we use the two-dimensional Dynamic Time Warping algorithm, implying that we compare two wells composed of N markers and M well markers respectively. Each pair of markers (i, j) has a correlation cost c(i, j), and each transition between two correlations (i, j) and (k, l) has a transition cost t j,l i,k as illustrated in Figure 5.

Correlation and transition costs are computed for every possible well marker correlation and correlation transitions. Assuming that extremities of the two wells are correlated, for each scenario, the global cost is computed iteratively from the top correlation (1,1) to the bottom correlation (N, M ) given by:

∀ (i, j) ∈ 1, N × 1, M , C(i, j) = c(i, j) + min    t j,j i-1,i + C(i -1, j) t j-1,j i-1,i + C(i -1, j -1) t j-1,j i,i + C(i, j -1)    . (1)
Figure 5 illustrates the two-dimensional Dynamic Time Warping algorithm.

Considering two wells composed of N markers and M markers respectively, it is possible to generate a two-dimensional matrix where each cell c(i, j) corresponds to the correlation (i, j). Each cell is filled by the correlation cost c(i, j) and is linked to its three up-and-left neighbors by the transition cost t j-1,j i-1,i

(red-arrows). The Dynamic Time Warping algorithm iteratively minimizes the global cost and keeps the previous path with the lowest cost C(i -1, j -1) (green-arrows). The algorithm ends when it reaches the last correlation (N, M ) (black-arrow).

Dynamic Time Warping was initially introduced in stratigraphic applications to correlate facies [START_REF] Smith | New Stratigraphic Correlation Techniques[END_REF][START_REF] Waterman | The match game: New stratigraphic correlation algorithms[END_REF] and can also be used to correlate well-logs [START_REF] Hale | Dynamic warping of seismic images[END_REF][START_REF] Wheeler | Simultaneous correlation of multiple well logs[END_REF][START_REF] Wu | Incremental correlation of multiple well logs following geologically optimal neighbors[END_REF]. For relative chronostratigraphy, the main difficulty is to define appropriate correlation costs c and transitions costs t. For example, Lallier et al. ( 2016) used paleobathymetry and depositional slope constraints to define carbonate facies correlation rules. [START_REF] Edwards | Uncertainty management in stratigraphic well correlation and stratigraphic architectures: A training-based method[END_REF] derived correlation and transition costs by scanning three-dimensional models obtained by processbased stratigraphic simulation. These approaches can incorporate depositional concepts into automatic correlation, but they call for determining quantitative parameters which may be difficult to infer in some depositional environments.

In the next subsection, we use the more qualitative notion of distality discussed in Section 1.3 to define new correlation costs c(i, j).

Correlation cost computation -c(i, j)

If we consider Figure 3 as a simple approximation for a sedimentary basin margin geometry along the depositional strike direction (AA') and depositional dip direction (BB'), it is possible to qualitatively estimate the well marker correlation likelihood as illustrated in Figure 4.

Input data are (1) sedimentary facies interpretation defined along the entire well by their top and bottom depths, and (2) well distality defined by a piecewise constant function of the depth according to the depositional context. The proposed method may be applied on two wells, referred to as well 1 and well 2. For each well, well markers are taken at the center of each sedimentary facies, and assigned to f directly given by the sedimentary facies encoding (facies codes increase from the deepest (f = 1) to the shallowest sedimentary facies). Well distality d may vary along a vertical well, but is assumed to be constant within a clinoform (distality is increasingly defined from the most distal position within the basin (d = 1) to the most proximal position within the basin), i.e., a parasequence. Note that this ordering can be reversed as long as facies depositional depth and distality are encoded consistently.

Considering that well 1 is composed of N markers m 1 (i), i ∈ 1, N and well 2 is composed of M markers m 2 (j), j ∈ 1, M , each marker m 1 (i) and m 2 (j) can be expressed as:

m 1 (i) = {f 1 (i), d 1 (i)} , m 2 (j) = {f 2 (j), d 2 (j)}. (2)
In order to make every correlation (i, j) comparable, their generation are based on adimensional terms f (i, j) and d(i, j) computed from normalization terms f 0 and d 0 computing respectively from the maximum facies and distality variations observed between the two wells:

f 0 = max i,j (|f 1 (i) -f 2 (j)|) , d 0 = max i,j (|d 1 (i) -d 2 (j)|). ( 3
)
The terms f (i, j) and d(i, j) are given by

f (i, j) = |f 1 (i) -f 2 (j)| f 0 , d(i, j) = α × |d 1 (i) -d 2 (j)| d 0 , ( 4
)
where α is a scaling coefficient representing how the lateral size of the depositional is deemed to scale with the inter-well distance.

Based on a conceptual sedimentary basin geometry illustrated in Figure 3 and on the principle of correlation considering that a facies cannot be associated with a depositionally shallower (respectively, deeper) facies in a more distal (respectively, proximal) well, it is finally possible to define two possible cases:

1. Facies and distality interpretations are inconsistent with the principle of correlation, this way the correlation cost is given by c(i, j) = +∞.

(5)

2. Facies and distality interpretations are consistent with the principle of correlation, then the correlation cost is given by (Figure 6)

c(i, j) = (d(i, j) -f (i, j)) 2 . ( 6
)
This cost function may be analyzed for several cases:

-d(i, j) = 0 (wells are orthogonal to the transport direction). In this case, facies are preferentially associated if their depositional depths are similar and the cost function is assimilated to the normalized facies variance given by:

c(i, j) = f (i, j) 2 = |f 1 (i) -f 2 (j)| f 0 2 . (7)
-d(i, j) = 1 (wells are parallel to the transport direction and have extreme distality values). In this case, facies are preferentially associated if their depositional depths are different, and the cost function is assimilated to the opposite of the normalized facies variance:

c(i, j) = (1 -f (i, j)) 2 = 1 - |f 1 (i) -f 2 (j)| f 0 2 . ( 8
)
-In all the other cases, the cost function ( 6) entails that the closer the distality, the higher the correlation likelihood between facies having similar depositional depths (Figure 6). 4) and c(i, j) is given by Equation ( 6). (single columncolored)

2.2. Transition cost computation -t j,l i,k

According to the cost matrix illustrated in Figure 5, each well marker correlation is preceded by one well marker correlation among three possible correlations (red arrows) and their transitions are more or less likely. Transition costs are included within the Dynamic Time Warping algorithm to constrain the result as expressed in Equation ( 1).

Currently, our method does not compute transition costs from input data. These vertical and horizontal transitions within the cost matrix correspond to gaps in correlation sets. In Dynamic Time Warping, gaps are naturally penalized as the diagonal path in the cost matrix minimizes the number of transitions. However, to further penalize stratigraphic gaps, we decide to define a non-zero constant gap cost given by t j-1,j i,i = t j,j i-1,i = 0.1 and t j-1,j i-1,i = 0.

In this case, 0 ≤ c(i, j)+t j-1,j i-1,i ≤ 1.1, and the current cost is no more normalized. A new correlation cost c * (i, j) is therefore given by c * (i, j) = 0.9 × c(i, j), so that 0 ≤ c * (i, j) + t j-1,j i-1,i ≤ 1 is attributed to admissible correlations.

Synthetic data set

The proposed correlation cost function is first tested on synthetic data illustrated in Figure 7. Each scenario is composed of five wells (1,2,3,4,5), so we apply the Dynamic Time Warping method incrementally by defining a correlation path and computing the costs between each added well and the closest well already processed [START_REF] Lallier | Uncertainty assessment in the stratigraphic well correlation of a carbonate ramp: Method and application to the Beausset Basin, SE France[END_REF].

In our examples, wells are ordered and correlation sets are computed between two successive wells, that is (1-2), then (2-3), then (3-4) and then (4-5). We generated the most likely solutions for each synthetic data set.

This synthetic data set is composed of three simple scenarios based on simple sedimentary basin geometries:

A. Sedimentary basin margin transversal cross-section. All the considered wells have the same distality, so the correlation cost c(i, j) is computed using Equation ( 7). The best correlation set having the lowest cumulative cost C = 0.00 is unique and corresponds to a layer cake correlation ( C. Sedimentary deltaic transversal cross-section. The distality increases from the center (1 ≤ d n (i) ≤ 3) and the correlation cost c(i, j) is computed using the distality costs ( 5) and ( 6). There are nine correlation sets which reach the lowest cumulative cost C = 2.58. Figure 7.C illustrates these realizations in black and red dashed lines. Outcomes are more difficult to evaluate than in scenarios A and B, but they are relatively close to the simplified subsurface model shown in the background.

The reference correlation set is recovered, but it has the same cost as some correlation sets which cross the reference time lines. This illustrates the indeterminacy of the interpretation based solely on incomplete data and on the proposed cost. In cases B and C, it is important to note that the number of realizations costing the less is smaller than expected because of the sequential nature of the correlation. Alternative correlations could be obtained by randomizing the correlation order.

3. Application: Hugin Formation, Gudrun-Sigrun Field, North Sea

The proposed method for computing the likelihood of a well marker correlation from sedimentary facies and well distality yields promising results on simple synthetic data sets and is now applied to a real data set from the North Sea. This data set consists of several wells drilled in the South Viking Graben targeting a Middle Jurassic sandstone reservoir. The interval of interest is the Hugin Formation (Early Callovian part) within the Gudrun-Sigrun Field area (North Sea Block 15/3) (Knaust and Hoth, 2020).

Stratigraphy and sedimentary context

Based on dipmeter data analysis and biostratigraphic interpretation, the Hugin Formation within the Gudrun-Sigrun Field area appears as an unfaulted chronostratigraphic unit. Its top and bottom are highlighted in seismic images by coal bed reflectors. However, the three-dimensional geometry of the Hugin Formation cannot be interpreted from seismic data and can only be determined from well data [START_REF] Knaust | Bay-head deltas as hydrocarbon reservoirs: The Middle Jurassic Hugin Formation in block 15/3 of the South Viking Graben, Norway[END_REF].

From drilled and partially cored wells in the zone of interest, well logs and core samples were analyzed to define eight sedimentary facies (Figure 8) and described by [START_REF] Knaust | Bay-head deltas as hydrocarbon reservoirs: The Middle Jurassic Hugin Formation in block 15/3 of the South Viking Graben, Norway[END_REF] as the bottom "part of idealized shallowingupward successions that are bounded by marine flooding surfaces interpreted and therefore are interpreted as parasequences [. . . ] and prograding delta lobes": The stratigraphic context suggests deposition in distributary channels" [START_REF] Knaust | Bay-head deltas as hydrocarbon reservoirs: The Middle Jurassic Hugin Formation in block 15/3 of the South Viking Graben, Norway[END_REF].

These facies codes are globally ordered by decreasing distality from 1 to 8, so they can be used directly in Equation ( 6). Moreover, facies can also be grouped before correlation to reduce complexity, and assess the impact of the level of interpretation detail on the correlation outcomes. To this end, the eight initial facies are grouped in three possible facies associations according to their depositional positions along a distal-to-proximal transect:

A. Sedimentary facies are divided into three groups. ( 1) Distal sedimentary facies 1-3 corresponding to depositional environments from restricted lagoons to patch reefs. (2) Intermediate sedimentary facies 4-6 corresponding to deltaic environments from the prodelta to the upper delta front.

(3) Proximal sedimentary facies 7-8 corresponding to mouth bars and distributary channels.

B. Sedimentary facies are divided into three groups. ( 1) Distal sedimentary facies 1-3 corresponding to depositional environments from restricted lagoons to patch reefs. (2) Intermediate sedimentary facies 4 corresponding the prodelta. (3) Proximal sedimentary facies 5-8 corresponding to depositional environments from the lower delta front to distributary channels.

C. Sedimentary facies are divided into two groups. ( 1) Distal sedimentary facies 1-4 corresponding to depositional environments from restricted lagoons to the prodelta. (2) Intermediate sedimentary facies 5-8 corresponding to depositional environments from the lower delta front to distributary channels.

Well data

The proposed correlation cost function based on sedimentary facies and well distality is applied on two wells, well 1 and well 2 from the Gudrun-Sigrun Field area. These two wells are fully cored in the bottom part of the Hugin Formation and biostratigraphic analysis performed. It confirms that the lower part of the Hugin Formation in wells 1 and 2 have the same geological age (i.e., Early Callovian) (Figure 9).

During the Early Callovian, deposits in wells 1 and 2 were part of the same bayhead delta sourced with sediments from East to Northeast [START_REF] Knaust | Bay-head deltas as hydrocarbon reservoirs: The Middle Jurassic Hugin Formation in block 15/3 of the South Viking Graben, Norway[END_REF]. Based on the depositional environment classification, sedimentary facies are interpreted along the lower part of the Hugin Formation in well 1, referred to as the distal well because of its western position, and in well 2, referred to as the proximal well because of its eastern position (Figure 9). It is assumed that distality is constant in each well (d 1 (i) = 1 and d 2 (j) = 2 in Equation ( 4)).

This assumption is based on a stable direction of sediment transport during the considered time interval, which is supported by dipmeter data.

Additionally, transgressive and regressive parasequences (T-R parasequences) were interpreted along these two wells. This correspondence is supported by biostratigraphic interpretation and by the identification of major marine flooding events characterized by reworked sediments and scattered peebles above these ravinement surfaces [START_REF] Knaust | Bay-head deltas as hydrocarbon reservoirs: The Middle Jurassic Hugin Formation in block 15/3 of the South Viking Graben, Norway[END_REF]. However, it is important to note that no absolute evidence exists about how these regressive sequences are correlated between both wells. The deepest parasequence identified in well 1 and the shallowest parasequence identified in well 2 (Figure 9) are analyzed by [START_REF] Knaust | Bay-head deltas as hydrocarbon reservoirs: The Middle Jurassic Hugin Formation in block 15/3 of the South Viking Graben, Norway[END_REF] and dipmeter data acquired in these parasequences confirm a sediment transport direction towards West-Southwest, i.e., from well 2 to well 1.

Results and interpretations

Simulating computer-assisted correlations on this data set involves computing a correlation cost c(i, j) for each possible well marker correlation. First, a correlation cost minimizing the facies variance is computed from the interpreted well markers using Equation ( 7), corresponding to lithostratigraphic correlation. Then, the distality correlation cost function is applied on the data using Equation (8), and on the three different facies associations defined in Section 2.3.

Figure 10: The correlation set having the lowest cumulative correlation cost computed using Equation ( 7), corresponding to a lithostratigraphic concept. Red correlation bands highlight well marker correlation, inconsistent with the principal sediment transport direction in the considered interval. (single column -colored)

Sedimentary facies variance minimization

In order to minimize the sedimentary facies variance (e.g., [START_REF] Hale | Dynamic warping of seismic images[END_REF], the well distality is not taken into account. Figure 10 shows the realization having the lowest cumulative correlation cost computed with Equation ( 7), i.e., considering that wells have the same distality, and using only lithostratigraphy.

Simulation outcome highlights inconsistent well marker correlations (red correlation bands on Figure 10), for example a well marker interpreted as a patch reef sedimentary facies in the well 1 (between depths z = 17.7 m and z = 18.59 m) is correlated with a well marker interpreted as an open lagoon in the well 2 (between depths z = 10.91 m and z = 11.03 m). These inconsistent well marker correlations go against the principle of correlation postulated in Section 1.3.

Sedimentary facies versus well distality

Well correlation computed by minimizing the sedimentary facies variance did not return convincing outcomes. Indeed, lithostratigraphic concepts do not, in general, yield isochronous stratigraphic surfaces. Using Equations ( 5) and ( 6) to account distality significantly changes the correlation outcomes as illustrated in Figure 11. The correlation panel illustrated in Figure 11 eases the comparison between the result of the two most likely correlation sets among the n-most likely scenarios computed using lithostratigraphic cost function (7) (dashed red lines) and the distality correlation rule ( 5) and ( 6) (full black lines).

Green correlation bands highlight that well markers which were inconsistently correlated using the sedimentary facies variance minimization (red correlation bands in Figure 10) are now consistently correlated using the proposed correlation cost function. This comes from the cost function ( 5). Moreover, we observe that several facies in one well are correlated with only one facies in the other well. The computer-assisted correlation process generates a gap due to erosion or non-deposition. The depths of these gaps are not fixed a priori by the algorithm which only considers a gap cost penalty for each transition, so they emerge from the method in order to maximize the correlation set likelihood. However, if the number of well markers in the different wells differs, multiple correlations to one single marker must occur, and gaps must be simulated.

Another way to visualize a well-to-well correlation problem is to plot a twodimensional cost matrix representing all well marker associations (Figure 12). 5) and ( 6). The lines correspond to the four most likely scenarios, i.e., the correlation sets with the four lowest cumulative costs computed using Equation (1). Note that the transition costs t j,l i,k associated to the edges or corners of the cost matrix are not visible in this display. (single column -colored)

On this matrix, we can observe that a high number of well marker correlations are excluded thanks to the principle of correlation implemented by Equation ( 5). Indeed, an infinite correlation cost means that the well marker correlation is not consistent. For example, the large mouth bar observed in the well 1 (between depths z = 9.75 m and z = 15.95 m) can only be associated with mouth bars (between depths z = 1.33 m and z = 2.54 m) and distributary channels (between depths z = 3.36 m and z = 4.17 m) in the well 2, which is more proximal. The principle of correlation constrains the computer-assisted well-correlation process and generates only admissible scenarios.

In Figure 12, the most likely correlation set is represented by the red line through the cost matrix and each bullet corresponds to each well marker correlation belonging to this correlation set. It is possible to display several scenarios, for example the four scenarios having the lowest cost in Figure 12. The cost matrix view eases the analysis of differences between correlation sets because we can directly target well marker correlations which are rarely or often simulated.

Moreover, the two-dimensional cost matrix highlights simulated gaps (vertical and horizontal lines) by applying the proposed correlation cost function on the North Sea data set returns two specific types of gaps in this particular scenario (Figure 12):

1. One well marker interpreted as a distributary channel in well 2 (between depths z = 9.75 m and z = 15.95 m) is correlated with a group of well markers in well 1 (between depths z = 2.52 m and z = 21.50 m). This significant gap stems from the lack of a proximal facies (distributary channel or mouth bar) below 4.5 m on well 2.

2. Two well markers interpreted as a restricted lagoon in well 2 (between depths z = 22.49 m and z = 23.88 m and between depths z = 24.51 m and z = 25.00 m) are correlated with two groups of well markers in well 1 (between depths z = 6.35 m and z = 10.37 m and between depths z = 10.91 m and z = 15.30 m).

These observations are consistent with the interpretation of regressive parasequences during which sediments are preferentially deposited along the basin margin rather than in the proximal or distal positions [START_REF] Galloway | Genetic stratigraphic sequences in basin analysis I: Architecture and genesis of flooding-surface bounded depositional units[END_REF][START_REF] Homewood | Corrélations haute résolution et stratigraphie génétique : une démarche intégrée[END_REF][START_REF] Borgomano | Stratigraphic well correlations for 3-D static modeling of carbonate reservoirs[END_REF][START_REF] Catuneanu | Sequence stratigraphy: Common ground after three decades of development[END_REF]. Hence, this realization, which is defined by the highest likelihood and based on the principle of correlation, yields consistent results for the sedimentary facies correlation, that is constrained by the well distality.

However, these results are not consistent with the stratigraphic parasequence interpretation (red triangles in Figure 9). Therefore, this method can help targeting zones of uncertainty in order to make alternative interpretations of core samples and well logs. Alternatively, some method parameters such as the relative well distality (α in Equation ( 4)) could be tuned to match ancillary constraints, before using the approach to explore uncertainties.

Impact of the sedimentary facies associations

The correlation cost function is very strict and excludes well marker correlation according to their distality and sedimentary facies interpretation. However, some depositional environments may be laterally equivalent and this equivalence may not be captured by the rules. Moreover, relative well distality may vary through time. This is why we tested it on different sedimentary facies associations (Section 2.3) in order to evaluate the impact of sedimentary facies lateral equivalence. Figure 13 plots one of the most likely computed correlation sets using the facies association A as input against the most likely computed correlation set without facies associations. The main difference between this outcome and the most likely correlation set based on the eight sedimentary facies (dashed red correlation lines) is highlighted by the purple correlation band. Without the facies association A, the well marker interpreted as a restricted lagoon in well 1 (between depths z = 8.65 m and z = 8.91 m) seems to pinch out between the wells (outcome is a gap), but with the facies association A, this well marker is correlated with a well marker interpreted as an upper delta front in well 2 (between depths z = 2.54 and z = 3.36 m). Correlation sets obtained with the facies associations B and C are illustrated in Figure 14.

The main difference between the outcome computed with 8 facies (Figure 11) and outcomes based on facies associations A, B and C is the number of simulated gaps. Assigning the same value to several facies presenting very close depositional conditions (lateral equivalence) leads to a smaller number of gaps and amounts to a more relaxed interpretation of the principle of correlation. Indeed, two facies belonging to the same group have the same code value even if their association goes against the principle of correlation.

For example, considering the facies association A, a well marker interpreted as patch reef (distal facies) in well 1 can be associated to a well marker interpreted as restricted lagoon (distal facies) in the proximal well (Figure 13). The number of plausible well marker correlations (number of painted cells within the cost matrix) should increase directly by associating facies according to their depositional conditions. Increasing the number of plausible correlation may increase the number of diagonal correlation transitions (no gap) and should lead to outcomes with less gaps because of their cost penalty (Section 2.2).

Discussion

The proposed method has been applied to two wells of the Gudrun-Sigrun Field area, but it can be extended to n wells as described in Section 2.3. In this case, wells must be classified by their relative positions. Normalization terms defined in Equations ( 3) and ( 4) are computed on the entire data set to ensure the consistency of the results. However, the order of well correlation may have a strong impact on the outcomes because only a limited set of scenarios can be propagated through the multi-well correlation process [START_REF] Edwards | Uncertainty management in stratigraphic well correlation and stratigraphic architectures: A training-based method[END_REF].

Independently of the well order, the application of the proposed correlation rule on the Gudrun-Sigrun data set shows that the definition of well distality and sedimentary facies values have a strong impact on the results:

Variation of the principal sediment transport direction. In sedimentary basins, the principal sediment transport direction may vary between two overlapping deltaic lobes, and the presence of several sediment sources may lead to a variation of the principal sediment transport direction in the system. In both cases, the variation of the principal sediment transport direction may be bracketed with reasonable accuracy with dipmeter data. Figure 10 shows the most likely correlation set considering that both wells 1 and 2 have the same distality, i.e., they are located along the strike direction. The associated correlation cost matrix to the correlation given by Equation (7) (i.e., lithostratigraphic correlation) is entirely colored by correlation cost because no well marker correlation is excluded. Figure 11 shows the most likely correlation set considering that wells 1 and 2 are located at the two extremities of the basin along the principal sediment transport direction, and the corresponding correlation cost matrix (Figure 12) illustrates that several well marker correlations are excluded thanks to the principle of correlation. Considering two wells, and considering two possible principal sediment transport directions in a circle of 360 degrees, a high number of correlation sets which are plausible in one case may become not possible in the other case. For example, if the relative well position along the distality axis switches, the outcomes might significantly change because of well marker association exclusion changes, thanks the principle of correlation.

Variation of the depositional system size with respect to well spacing. Results are also sensitive to the distality variation along a proximal-to-distal transect, which cannot be determined with certainty. To visualize the effects of changing 7) and ( 8) with α = 2/7. The black line corresponds to the most likely correlation set and the red line corresponds to the most likely correlation set computed with α = 1 (Figure 12). (single column -colored)

the relative well distality, we computed the cost matrix assuming a depositional system of lateral size equal to 3.5 km, corresponding to a scaling parameter α = 2/7 (Figure 15). Comparing the costs of Figures 11 and15, we see that the distality rule is less discriminant and tends to generate correlation paths closer to the diagonal, i.e., corresponding to a more steady preservation rate in each well in the considered interval.

Both uncertainties about the variation of the principal sediment transport direction and the relative depositional system size can be addressed by playing with the scaling coefficient α in Equation (4):

-A variation of the principal sediment transport direction may be modeled by a scaling coefficient varying between α = -1 and α = 1. A scaling coefficient α = -1 corresponds to the inversion of the principal sediment transport direction and a scaling coefficient α = 0 corresponds to the principal sediment transport direction being perpendicular to the wells.

-A variation of the depositional system size with respect to well spacing may be modeled by a scaling coefficient varying between α = 0 and α = 1.

A scaling coefficient α = 1 corresponds to the extreme position of wells within the positional system and a scaling coefficient α = 0 corresponds to two wells having the same depositional position within the system.

Modification of sedimentary facies associations. Figures 11, 13 and 14 show chronostratigraphic correlation outcomes computed with Equation (8) using three different facies associations and no facies association. Considering groups of lateral equivalent facies, the correlation cost matrix will be modified because two well markers interpreted as lateral equivalent sedimentary facies are allowed to be correlated. The correlation cost function takes as input parameters the sedimentary facies variation f (i, j) and the well distality variation d(i, j). In the case that two sedimentary facies are laterally equivalent, the variation f (i, j) is null and the correlation cost function will always compute a correlation cost between the two given well markers. Figure 16 shows the correlation cost matrix associated to the facies association A. This highlights that more plausible well marker correlations exist than when considering the initial 8 facies. For instance, well markers interpreted as prodeltas in the proximal well can be correlated to well markers interpreted as upper or lower delta fronts in the distal well (Figure 16).

This computer-assisted method aims to retrieve possible chronostratigraphic lines, and may identify facies association and gaps which are consistent with sedimentary concepts. For example, very thin distal facies or proximal facies may be associated with a thicker succession of intermediate facies during marine regression [START_REF] Homewood | Corrélations haute résolution et stratigraphie génétique : une démarche intégrée[END_REF]. Considering the consistency of our method applied on sedimentary coastal deltaic data set, it may be used to confront manual parasequence interpretations from well-logs and core-samples, and highlight where uncertainties exist.

However, the correlation set simulation is solely based on correlation cost functions used by the Dynamic Time Warping, so the results are only as good as the input rules [START_REF] Hale | Dynamic warping of seismic images[END_REF][START_REF] Wheeler | Simultaneous correlation of multiple well logs[END_REF][START_REF] Lallier | Uncertainty assessment in the stratigraphic well correlation of a carbonate ramp: Method and application to the Beausset Basin, SE France[END_REF][START_REF] Edwards | Uncertainty management in stratigraphic well correlation and stratigraphic architectures: A training-based method[END_REF], and in our case, the likelihood of the principle of correlation. For example, the correlation process does not take into consideration the erosion surfaces observed at the top of regressive parasequences. Therefore, additional rules could be integrated in the method. Based on sequence stratigraphy, and considering that during a regressive period sediments are mainly deposited in the intermediate area of the basin, and that during a transgressive period sediments are eroded in the distal part and mainly deposited in the proximal part of the basin as qualitatively illustrated in Figure 2, a transition cost t j-1,j i-1,i may be implemented based on the facies vertical transitions. The thickness of facies intervals could also be considered to further constrain the preservation rate 8) (facies versus distality). The black line corresponds to the most likely correlation set (Figure 13) and the red line corresponds to the most likely correlation set computed without facies association (Figure 12). (single column -colored) emerging from the correlations [START_REF] Lallier | Management of ambiguities in magnetostratigraphic correlation[END_REF][START_REF] Hagen | Dynamic time warping of palaeomagnetic secular variation data[END_REF].

More generally, the proposed cost normalization, 0 ≤ c(i, j) ≤ 1, makes it possible to define multi-criteria correlation costs. Indeed, several independent normalized costs may be combined in the Dynamic Time Warping algorithm by a simple multiplication. As another example, detailed biostratigraphic interpretations are available in the South Viking Graben data set and may be used to constrain the correlation by adding a binary cost function, which allows or excludes a correlation considering the biostratigraphy [START_REF] Caumon | Stochastic stratigraphic correlation of multiple wells seen as a directed acyclic graph creation problem[END_REF].

Additional data or concepts could also be integrated in the method, preferably as cost function to be combined to the proposed correlation cost function, as an ancillary geological likelihood data to evaluate the results (e.g., [START_REF] De La Varga | Structural geologic modeling as an inference problem: A Bayesian perspective[END_REF]. Posterior analysis on the generated stratigraphic model may be achieved to evaluate the physical likelihood of the correlation outcomes (layer connectivities) like hydraulic or tracer simulation tests (e.g., [START_REF] Baville | From well logs to stratigraphic layering: Automation, uncertainties and impact on reservoir behavior[END_REF]. Another way to scrutinize and improve the results could be to analyze the geometry of the chronostratigraphic units and their consistency with dipmeter data and subsidence analysis.

Conclusion

In this paper, the purpose of the proposed correlation cost function is to generate several plausible higher-order chronostratigraphic scenarios between lower-order interpreted time lines (e.g., from seismic imaging and biostratigraphic analysis). The obtained correlation results are simulated based on sedimentary facies and well distality interpretations in several wells, and have implications in terms of detailed paleogeography in the considered interval.

The proposed correlation rule is based on the principle of correlation: a facies cannot be associated with a depositionally shallower (deeper) facies in a more distal (proximal) well. This method requires sedimentary facies interpretations along well paths and the distality of all wells computed from the well position along a proximal-to-distal transect, to compute a correlation cost for each possible well marker correlation given by Equations ( 5) and ( 6).

The overall goal of the method is to generate time lines which reflect both allo cycles and auto cycles. Given the incomplete information available to recover such time lines, ambiguity is clearly expected and it should not be hidden, as it may affect the understanding and the way applied geoscientists model facies, flow units and thin horizontal permeability barriers in subsurface reservoirs [START_REF] Lallier | Uncertainty assessment in the stratigraphic well correlation of a carbonate ramp: Method and application to the Beausset Basin, SE France[END_REF]. The approach proposed in this paper makes it possible to include the concept of distality in such uncertainty quantification and to generate alternative layer connectivities between linear stratigraphic sections.

CRediT author statement

Paul Baville: Conceptualization, Investigation, Methodology, Software, Visualization, Writing (original draft). Marcus Apel: Conceptualization, Data Curation, Funding Acquisition, Resources, Validation, Writing (reviewing and editing). Silvan Hoth: Conceptualization, Data Curation, Funding Acquisition, Resources, Validation, Writing (reviewing and editing). Dirk Knaust: Conceptualization, Data Curation, Resources, Validation, Writing (reviewing and editing). Christophe Antoine: Software. Cédric Carpentier: Supervision, Writing (reviewing and editing). Guillaume Caumon: Conceptualization, Methodology, Supervision, Funding Acquisition, Validation, Writing (reviewing and editing).

Figure 1 :-

 1 Figure 1: Lithostratigraphic (left) and Chronostratigraphic (right) facies correlation and geological layers interpolation from the same data set. Inspired by Ainsworth et al. (1999). (single column -colored)

Figure 3 :

 3 Figure 3: Three-dimensional illustration of a simplified sedimentary basin margin geometry along the depositional strike direction (AA') and the principal sediment transport direction (BB'). Parasequence boundaries (blue lines) are isochronous lines separating clinoforms. (single column -colored)

Figure 4 :

 4 Figure 4: Well marker correlation likelihood based on facies association within a sedimentary basin along the depositional strike direction (AA') and the principal sediment transport direction (BB'). Parasequence boundaries (blue lines) are isochronous lines separating clinoforms. Facies associations are represented by dashed or full green or red lines according to their likelihood. (single column -colored)

Figure 5 :

 5 Figure 5: Illustration of the Dynamic Time Warping algorithm at step (i, j) with two wells containing N and M markers. (single column -colored)

Figure 6 :

 6 Figure6: Cost function computing the correlation cost between two well markers based on relative sedimentary facies and relative well distality variations. f (i, j) and d(i, j) are normalized parameters given by Equation (4) and c(i, j) is given by Equation (6). (single columncolored)

Figure 7 :

 7 Figure 7: Application of the proposed cost function on three synthetic data sets. (A) Sedimentary basin margin transversal cross-section. (B) Sedimentary basin margin or deltaic longitudinal cross-section. (C) Sedimentary deltaic transversal cross-section. For each scenario, we consider only scenarios having the lowest cumulative cost represented by dashed lines. Black lines are correlations occurring in all realizations, and red lines are correlations which are not systematic. (single column -colored)

  Figure 7.A). The outcome is consistent because all facies are similar along correlation lines. Hypothetical parasequence boundaries may be interpreted as horizontal lines. B. Sedimentary basin margin or deltaic longitudinal cross-section. The distality increases along the distal-to-proximal transect (1 ≤ d n (i) ≤ 5) and the correlation cost c(i, j) is computed using Equation (6). There are two correlation sets which reach the lowest cumulative cost C = 0.79 (Figure 7.B). Black dashed lines are correlation lines which occur in the two scenarios, and red dashed lines show differences between the two equivalent realizations. Outcomes are consistent because along the sediment transport direction, depositionally deep facies in distal wells are correlated with depositionally shallow facies in proximal wells. These two realizations return three separated correlation groups (connected correlation lines), and hypothetical parasequence boundaries may be interpreted as lines separating two successive chronostratigraphic layers.

Figure 8 :

 8 Figure 8: Depositional environment of different sedimentary facies interpreted by Knaust and Hoth (2020) from core samples of the Hugin Formation (Bathonian-Callovian) in two wells in the South Viking Graben, North Sea. (single column -colored)

Figure 9 :

 9 Figure9: Sedimentary facies, parasequence and biostratigraphic interpretation along wells 1 and 2 in the lower part of the Hugin Formation within the Gudrun-Sigrun Field area[START_REF] Knaust | Bay-head deltas as hydrocarbon reservoirs: The Middle Jurassic Hugin Formation in block 15/3 of the South Viking Graben, Norway[END_REF]. (single column -colored)

Figure 11 :

 11 Figure11: The correlation set having the lowest cumulative correlation cost computed using Equation (8). Green correlation bands highlight well marker correlations which were inconsistent with the sedimentary facies variance minimization and consistent with the proposed rule. Dashed red correlation lines correspond to the most likely correlation set based on the variance minimization (Figure10). (single column -colored)

Figure 12 :

 12 Figure12: Cost matrix C[i, j] corresponding to the two considered wells. Each cell representing a correlation between two well markers is colored by its correlation cost c(i, j) computed with Equations (5) and (6). The lines correspond to the four most likely scenarios, i.e., the correlation sets with the four lowest cumulative costs computed using Equation(1). Note that the transition costs t j,l i,k associated to the edges or corners of the cost matrix are not visible in this display. (single column -colored)

Figure 13 :

 13 Figure13: One correlation set having the lowest cumulative correlation cost from the facies association A. The purple correlation band highlights the main difference between this outcome and the outcome obtained without the facies association. Dashed red correlation lines correspond to the most likely correlation set computed without facies association (Figure11). (single column -colored)

Figure 14 :

 14 Figure 14: One correlation set having the lowest cumulative correlation cost from the facies association B (A) and from the facies association C (B). Dashed red correlation lines correspond to the most likely correlation set without facies association (Figure 11). (single column -colored)

Figure 15 :

 15 Figure15: Correlation cost matrix corresponding to chronostratigraphic correlation given by Equation (6) (facies versus distality). This case is intermediate between Equations (7) and (8) with α = 2/7. The black line corresponds to the most likely correlation set and the red line corresponds to the most likely correlation set computed with α = 1 (Figure12). (single column -colored)

Figure 16 :

 16 Figure16: Correlation cost matrix corresponding to chronostratigraphic correlation (Classification A) given by Equation (8) (facies versus distality). The black line corresponds to the most likely correlation set (Figure13) and the red line corresponds to the most likely correlation set computed without facies association (Figure12). (single column -colored)

Acknowledgments

This work was performed in the frame of the RING project at Université de Lorraine (http://ring.georessources.univ-lorraine.fr). We would like to thank the industrial and academic sponsors of the RING Consortium managed by ASGA and Equinor for their support. The code corresponding to this paper was implemented as a stand-alone C++ and Python algorithm (WeCo) available to RING Sponsors. The authors are grateful to Equinor and the PL025 licence partners Neptune, OMV and Repsol for providing permission to use the Gudrun-Sigrun Field area data set in this work and to publish it.


A Yvan

Guifo Fodjo
email: a.guifofodjo@etu.upmc.fr

Jerry Lacmou Zeutouo
email: jerry.lacmou-zeutouo@inria.fr

Samuel Bowong
email: sbowong@gmail.com

Separation of Concerns in an Edge-Based Compartmental Modeling Framework

Keywords: Separation of concerns, Compartmental models, Contact network, Epidemiology modeling tool, Edge-based compartmental network

A well-known framework with strong potential for epidemic prediction and the ability to incorporate realistic contact structures is edge-based compartmental modeling (EBCM). However, models built from this framework lead to a multiplication of ordinary differential equations and many parameters to be estimated, which make the models complex and difficult to extend or to reuse. The Kendrick approach has shown promising results in generalizing compartmental models to take into account aspects of contact networks while preserving the separation of concerns, thus allowing to define modular, extensible and reusable models. But this generalization of compartmental models to contact network aspects is still limited to a few contact networks. In this paper, we present an attempt to extend Kendrick's approach from an approximation of EBCM models to further support aspects of contact networks, thereby improving the predictive quality of models with significant heterogeneity in contact structure, while maintaining the simplicity of compartmental models. This extension consists of an integration of the basic reproductive number R 0 into the compartmental SIR framework. This attempted is validated using Miller's mass action and the approximation of EBCM configuration model.

Introduction

Mathematical modeling and computer simulation have been widely used in epidemiology to find appropriate control strategies and means of control [START_REF] Levin | From individuals to epidemics[END_REF]. One of the most common models in mathematical modeling using the compartmental framework is the mass action model (see Figure 1). The basic assumption in this model is that susceptible and infectious individuals meet at random and can spread the disease [START_REF] Cuddington | Ecological paradigms lost: routes of theory change[END_REF][START_REF] Keeling | Modeling infectious diseases in humans and animals[END_REF]. Susceptible individuals (S) can become infected from the force of infection 1 λ(S, I, N) = βI where β is the contact transmission rate and N is the population size. Infectious individuals (I) can be recovered at the recovery rate γ. Compartmental models are typically first defined as ordinary differential equations (ODEs) such as Equation 1. These models can be studied analytically and/or simulated using algorithms such as RungeKutta. However, it is considered more realistic to adopt a stochastic viewpoint on these models considering them as Continuous-Time Markov Chains (CTMCs). The latter can be derived from the ODEs modulo some widely accepted, albeit simplifying, probabilistic assumptions.

Although the mass action model is known largely for its conceptual and mathematical simplicity [START_REF] Kermark | Contributions to the mathematical theory of epidemics. part i[END_REF], it has been found to have some shortcomings, including the fact that it incorrectly assumes that all individuals have the same contact rate and implicitly assumes that all partnerships are infinitely short [START_REF] Miller | Edgebased compartmental modelling for infectious disease spread[END_REF][START_REF] Wang | Edgebased epidemic dynamics with multiple routes of transmission on random networks[END_REF]. [START_REF] Miller | Edgebased compartmental modelling for infectious disease spread[END_REF][START_REF] Kiss | Mathematics of epidemics on networks[END_REF][START_REF] Wang | Edgebased epidemic dynamics with multiple routes of transmission on random networks[END_REF] built on these limitations of the mass action model to introduce the edge-based compartmental modelling (EBCM) framework. It is a compartmentbased model extension approach where the contacts between individuals are materialized by the edges [START_REF] Miller | Edgebased compartmental modelling for infectious disease spread[END_REF][START_REF] Wang | Edgebased epidemic dynamics with multiple routes of transmission on random networks[END_REF]. The main idea is to integrate the heterogeneous mixture by considering the population as a network of individuals where the contacts follow a given probability distribution. Thus, it consists in incorporating into the ordinary differential equations of the mass action model, a kind of social heterogeneity (heterogeneous contact rate), while taking into account the impact of the partnership duration2 .

The results obtained from this incorporation have enabled realistic and more predictable contact structures to be considered. However, the modification of the compartmental framework of mass action results in more and more multiple and complex differential equations, especially when new concerns (age, sex, control strategy, etc.) and parameters are taken into account [START_REF] Balde | Contact vaccination study using edge based compartmental model (ebcm) and stochastic simulation: An application to oral poliovirus vaccine (opv)[END_REF]. Moreover, another difficulty of the EBCM framework is that the models it contains are not easily scalable, extensible and therefore not reusable.

Kendrick's approach [START_REF] Bui | Separation of concerns in epidemiological modelling[END_REF][START_REF] Bui | The Kendrick modelling platform: language abstractions and tools for epidemiology [+ correction art[END_REF] has shown its ability to define concerns (age, sex, spatial heterogeneity, etc.) as independent (possibly incomplete) models that are then combined into stochastic automata networks (SANs) [START_REF] Plateau | Stochastic automata networks[END_REF] using a tensorial sum operator. Stochastic dependencies between concerns are then introduced in a second phase so that independent concerns can be reused and combined in other models much more easily. In [START_REF] Fodjo | Separation of concerns in extended epidemiological compartmental models[END_REF], the authors showed that this approach, based on compartmental models, could be extended to take into account certain aspects of contact networks while enabling the building of reusable models. But this generalization of compartmental models to contact network aspects is still limited to a few contact networks (Poisson, Exponential and Scale free).

While it is well known that taking into account as-pects of contact networks in compartmental models leads to more realistic and predictive models, there is also a large body of work that recognizes the crucial importance of the basic reproductive number R 0 in predicting epidemics [START_REF] Meyers | Network theory and sars: predicting outbreak diversity[END_REF][START_REF] Danon | Networks and the epidemiology of infectious disease[END_REF][START_REF] Molina | Modelling the spread of diseases in clustered networks[END_REF][START_REF] Heesterbeek | Modeling infectious disease dynamics in the complex landscape of global health[END_REF][START_REF] Zhang | Modeling epidemics spreading on social contact networks[END_REF][START_REF] Trapman | Inferring r 0 in emerging epidemics-the effect of common population structure is small[END_REF][START_REF] Yang | The computational approach for the basic reproduction number of epidemic models on complex networks[END_REF]. This parameter refers to the number of new cases caused by a single randomly infected individual in a completely susceptible population. When R 0 < 1, epidemics are impossible, while when R 0 > 1, they are possible. [START_REF] Aparicio | Building epidemiological models from r0: an implicit treatment of trans-mission in networks[END_REF] suggest modifying the SIR compartmental framework to incorporate the parameter R 0 into the ordinary differential equations. In this study, in order to take the aspects of contact networks, R 0 is approximated to Poisson, Exponetial and Scale free contact networks.

In this work, we propose an attempt to extend Kendrick's approach from an approximation of the EBCM approach while maintaining the separation of concerns and preserving the simplicity of compartmental models. This approximation consists in constructing a concern in the sense of Kendrick's approach (i.e. as a stochastic automaton that can then be combined) from the R 0 of each EBCM model. The extension of the Kendrick approach is done by incorporating the base reproduction rate R 0 into the compartmental SIR framework. The simulation results obtained are similar to those of the EBCM configuration model presented in [START_REF] Miller | Edgebased compartmental modelling for infectious disease spread[END_REF].

Miller et al.'s mass action model and edge-based compartmental configuration models

In this section, and for the purposes of this work, we will restrict ourselves to the standard mass action model of [START_REF] Miller | Edgebased compartmental modelling for infectious disease spread[END_REF] and the edge-based compartmental configuration model.

Miller et al.'s mass action model

This model is constructed like the system of Equations 1) with some modifications. Such as the fact that the authors of [START_REF] Miller | Edgebased compartmental modelling for infectious disease spread[END_REF] assume that an infected individual causes new infections at the rate βS(t), where β is the transmission rate per infected individual. Recovery occurs at the rate of γ.

   S ′ = -βIS I ′ = βIS -γI R ′ = γI (2)
Figure 2: Flow diagram of the EBCM configuration model of [START_REF] Miller | Edgebased compartmental modelling for infectious disease spread[END_REF].

Moreover, when building the mass action model (see system of Equations 1), many works [START_REF] Cuddington | Ecological paradigms lost: routes of theory change[END_REF][START_REF] Keeling | Modeling infectious diseases in humans and animals[END_REF][START_REF] Martcheva | An introduction to mathematical epidemiology[END_REF] interpret the transmission rate β as the product of the contact rates and the transmission probability. But in the specific case of the mass action model of Miller et al.'s, the authors assume that the transmission rate per infected person β = β⟨k⟩ where ⟨k⟩ is the average degree of the contact network considered.

EBCM configuration model

The flow diagram of the EBCM configuration model is shown in Figure 2. The compartments S, I, and R represent the proportions of susceptible, infected, and recovered, respectively, as in the case of the mass action model in section 2.1. To calculate S(t), I(t), and R(t), the authors note that these are the probabilities that a random test node u is in each compartment. We calculate S(t) noting that this is also the probability that none of u's partners has yet transmitted to u. The probability that a randomly selected partner v has not yet transmitted the infection to u.

For large networks, [START_REF] Miller | Edgebased compartmental modelling for infectious disease spread[END_REF][START_REF] Wang | Edgebased epidemic dynamics with multiple routes of transmission on random networks[END_REF] assumed that the neighbors of the node test u are independent. Given a degree k, u is susceptible at time t with probability s(k, θ(t)) = θ(t) k . Thus, S(t) = ∑ k P(k)s(k, θ(t)) = ψ(θ(t)). The ODEs of the EBCM configuration model are obtained in (3).

S = ψ(θ) I = 1 -S -R (3)
The φ S , φ I and φ R compartments of Figure 2 represent respectively the probabilities that a partner v is susceptible, infected and recovered but has not yet transmitted the infection to u. The compartment 1-θ, represents the probability that there is infection (it is done with the rate βφ I between the compartment φ I and 1θ).

In order to determine θ, consider that (θ = φ S + φ I + φ R ). At time t = 0, θ = φ S ≃ 1. Because we consider that there has not yet been an infection, so

φ I = φ R ≃ 0.
The central parameter of this calculation is the de-termination of φ I . Thus,

φ I = θ -φ S -φ R (4)
It remains to calculate φ S and φ R explicitly to obtain φ I . If we finally consider that the v is infected and that there was an infection, then

θ ′ = -βφ I (5)
To determine φ R , the authors use the fact that in Figure 2 the fluxes from φ I to φ R and from φ I to (1θ) are proportional to one another. Both φ R and (1θ) are equal to zero at time zero since we assume that no infection or recovery events can occur prior to time zero. By integrating the relation

dφ R dt = γ β d(1 -θ) dt (6)
and using the initial condition

φ R (0) = (1 -θ(0)) = 0 (7) Thus, φ R = γ(1 -θ) β (8)
To determine φ S , the authors rely on the fact that a partner v has a degree k with probability P n (k) = kP(k)/⟨k⟩. Given a degree k, v is susceptible with probability θ (k-1) . This allows us to obtain

φ S = ∑ k P n (k)θ (k-1) = ∑ k kP(k) ⟨k⟩ θ (k-1) = ψ ′ (θ) ψ ′ (1) (9)
Therefore, from the relation 4, they obtain

φ I = θ -φ S -φ R = θ - ψ ′ (θ) ψ ′ (1) - γ(1 -θ) β ( 10
)
While it is easy to recognize that the EBCM configuration model captures much more population structure than the mass action model, it is also important to note that the complexity and number of differential equations increases significantly. As we have seen, to obtain Equation ( 10), the authors of [START_REF] Miller | Edgebased compartmental modelling for infectious disease spread[END_REF] had to develop new Equations 4, 5, 6, 7, 8, 9 in addition to the one of the system of (3).

Generalizing Kendrick's approach

The challenge of this work is to show that Kendrick's approach can be extended to support aspects of contact networks using an approximation of the EBCM approach. This would improve the predictive quality of models with significant heterogeneity in the structure of contacts while preserving the separation of concerns. This would result in modular models, easily extensible and reusable. The idea is to avoid the designers/modelers/developers having to build a model with a multitude of ordinary differential equations, having to build an explicit contact network and having to estimate many parameters. [START_REF] Fodjo | Separation of concerns in extended epidemiological compartmental models[END_REF] showed that it was possible to extend the compartmental framework in order to integrate some aspects of contact network models. Their idea was to make the force of infection a central parameter that could be redefined from the extension points (α gen3 , it gen4 , τ gen5 ) and according to the concerns that one wished to model. Thus, from the idea of [START_REF] Bansal | When individual behaviour matters: homogeneous and network models in epidemiology[END_REF] and the application of the Template Method Design Pattern [START_REF] Gamma | Elements of reusable object-oriented software[END_REF], they proposed a generic definition to the force of infection named λ (see Equation 11). [START_REF] Fodjo | Separation of concerns in extended epidemiological compartmental models[END_REF] insisted on the fact that the index "gen" of the extension points α gen , it gen and τ gen is intended to signify that they are applied to a model based on a network of contacts and to indicate that they are generic points, i.e. variable points. λ = α gen it gen τ gen (11)

For our modeling and simulation purposes, we use the same convention as the authors of [START_REF] Fodjo | Separation of concerns in extended epidemiological compartmental models[END_REF]. Thus we use names with the subscript "gen" in the final generic definitions of the models and in the Kendric code.

The force of infection of the mass action model given by (2) described in section 2.1 is given by the relation 12.

λ = βI = β⟨k⟩I (12) From the relations 11 and 12 the identification of the extension points α gen , it gen and τ gen gives :

   α gen = ⟨k⟩ it gen = I τ gen = β (13)
In section 2.2, we have presented the EBCM configuration model and the array of ODEs and parameters to be estimated from the model. This constitutes more effort of understanding, programming for the model designers/modelers/developers. Our goal is to cope with the complexity of ordinary differential equations in the EBCM framework while defining extensible and reusable models.

In order to cope with the complexity of the ordinary differential equations of the EBCM framework while defining extensible and reusable models, we propose to use an approximation of each EBCM model from the mathematical formulation of R 0 .

Indeed, [START_REF] Aparicio | Building epidemiological models from r0: an implicit treatment of trans-mission in networks[END_REF] showed that it was possible to modify the simple compartmental framework in order to integrate the parameter R 06 into the model equations.

We propose to this effect to approximate the EBCM models, by incorporating each R 0 of EBCM model in the compartmental framework SIR (Susceptible, Infected, Recovered) as shown in [START_REF] Aparicio | Building epidemiological models from r0: an implicit treatment of trans-mission in networks[END_REF], so as to define these EBCM model approximations as concerns in the sense of the Kendrick approach (i.e. a stochastic automaton that can be composed) as illustrated in [START_REF] Fodjo | Separation of concerns in extended epidemiological compartmental models[END_REF] and then to be able to determine the force of infection 11, which it would then be possible to decompose into different points of extensions α gen , it gen and τ gen .

Thus, the results of this incorporation of R 0 into the SIR compartmental framework would give ( 14) : [START_REF] Miller | Edgebased compartmental modelling for infectious disease spread[END_REF] gave the following formulation in the case of the EBCM configuration model :

   S ′ = -R 0 SI I ′ = R 0 SI -γI R ′ = γI (14)
R 0 = β (β + γ) ⟨k 2 -k⟩ ⟨k⟩ (15)
From ( 14), the force of infection of the EBCM approximation is given by the relation

λ = R 0 I = β (β + γ) ⟨k 2 -k⟩ ⟨k⟩ I (16)
Thus from the relation 11 the identification of the extension points allows us to obtain :

           α gen = ⟨k 2 -k⟩ ⟨k⟩ I gen = I τ gen = β (β + γ) (17)
In the relationship 17, the term α gen = ⟨k 2 -k⟩ ⟨k⟩ is obtained from the python code of the EoN module 7 for each degree distribution (Homogeneous, Poisson, Bimodal and PowerLaw) specified for the graph generation. For example, the code to obtain α gen from the configuration model generated from a Bimodal degree distribution 8 . We proceed in a similar way for distributions of degree Homogeneous 9 and Power-Law degree distribution 10 , but for the particular case of the Poisson degree distribution, ⟨k 2 -k⟩ = ⟨k⟩ 2 [START_REF] Miller | Edgebased compartmental modelling for infectious disease spread[END_REF] thus,

α gen = ⟨k 2 -k⟩ ⟨k⟩ = ⟨k⟩ 2 ⟨k⟩ = ⟨k⟩ .

Validation and discussion

To validate our approach, we replicated the EBCM configuration model experiments of [START_REF] Miller | Edgebased compartmental modelling for infectious disease spread[END_REF] using the approximation of the EBCM configuration model on the distributions 11 (Homogeneous, Poisson, Bimodal, and PowerLaw) in Kendrick [START_REF] Bui | The Kendrick modelling platform: language abstractions and tools for epidemiology [+ correction art[END_REF][START_REF] Fodjo | Separation of concerns in extended epidemiological compartmental models[END_REF]. In our implementations, the variable points of the λ parameter of ( 11) are named "alphagen", "itgen" and "taugen". The models are built in different entities that are easy to define in modules, extensible and reusable. Indeed, each model generally includes an entity of definition of the concerns then of composition of model (see Figure 6), an entity of initialization of parameters (see Figure 7) and finally an entity of simulation (lines 36 to 41) and visualization (lines 43 to 49) (see Figure 8).

The challenge was to check if Kendrick's approach to separation of concerns could capture approaches such as Miller et al.'s mass action model or the approximation of the EBCM configuration model while maintaining the familiar compartmental framework.

The implementation of the mass action model (see (2)) allowed us to obtain the black curve of Figure 9b. For this implementation, the code is subdivided into two major parts. The first one is the definition of the basic concern SIR (Susceptible, Infected and Recovered) where λ is defined (see Figure 3). In this implementation, λ is defined in a general way (see lines 9 to 10). The second part of the code, is the implementation of the mass action concern named here "maConcern" as shown in Figure 4. In this definition, lines 16 to 22 allow to redefine the extension points "alphagen", "itgen" and "taugen". Then in line 24 we compose the basic concern and the mass action concern.

As for the implementation of the approximation of the EBCM configuration model in Kendrick, we proceed in a similar way to the mass action model. But in the case of this model (see Figure 6), the extension point "alphagen" is obtained according to the specified degree distribution. In the case of the Bi-modal degree distribution for example, the value of "alphagen" is given in line (17 to 18) of Figure 6. Note that in our implementations, for approximation of the EBCM configuration model what varies from one degree distribution to another is the value of "alphagen". After having defined the concerns of the basic SIR model (see Figure 3) and the approximation of the EBCM configuration model (see Figure 6), we compose the different concerns as illustrated in Figure 5. The results obtained enable us to have in Figure 9b the curves of the various distributions of degrees of the model of configuration EBCM implemented in Kendrick.

It can be seen that, whether it is the mass action model of Miller et al.'s. or the approximation of the EBCM configuration model, Kendrick's approach allows us to easily define the concerns ("maConcern" and "cmBimodal" in the case of the approximation of the EBCM configuration model on a Bimodal degree distribution) in a way that is separate from the basic SIR concern "sirConcern". Therefore, the basic concern "sirConcern" can be reused without "maConcern" or "cmBimodal". Kendrick's approach enabled us to define "myConcern" and "cmBimodal" as independent models that can then be composed.

Regarding the results obtained while implementing the EBCM configuration model on the Power-law degree distributions (see orange curves of Figures 9a and9b), the same dynamics are observed but with a difference regarding the infectious. This can be explained by the fact that some individuals (superspreaders) have an abnormally high number of contacts at the beginning of the epidemic. We often note a strong propensity of individuals to attach themselves preferentially to the individual with more contacts in this type of network. We also note that the moment of the epidemic peak is reached at the same time.

EBCM configuration model simulations on homogeneous, Poisson, and bimodal degree distributions using our EBCM model approximation approach (see Figure 9b) yield improved curve heights (lower curves) compared to a typical configuration model approach on homogeneous, Poisson, and bimodal degree distributions simulated with EoN (see Figure 9a). Furthermore, we note that the curves of Figure 9b have the same dynamics as those seen in Figure 9a and significantly improve the predictive quality that one would expect from a mass action model (black curve). We also note that the timing of the epidemic peak of the homogeneous, Poisson and bimodal curves in Figure 9b is slightly slower compared to the respective homogeneous, Poisson, and bimodal curves in Figure 9a.

Our approach to approximating EBCM models can be applied to mean field social heterogeneity (MFSH) algorithms12 and a fixed-degree dynamic model (DFD)13 [START_REF] Miller | Edgebased compartmental modelling for infectious disease spread[END_REF][START_REF] Istvan | Mathematics of Epidemics on Networks: From Exact to Approximate Models[END_REF].

Conclusion

In this paper, we proposed to generalize Kendrick's approach to the consideration of realistic contact structures from the edge-based compartmental modeling (EBCM) framework while preserving the separation of concerns in compartmental epidemic models [START_REF] Bui | Separation of concerns in epidemiological modelling[END_REF]. To do this, we applied the solution of [START_REF] Fodjo | Separation of concerns in extended epidemiological compartmental models[END_REF] which involved defining the usual λ parameter of epidemic models as a kind of model method with three extension points, which allowed us to easily capture aspects of contact network models.

In order to validate this integration, we applied the approach of [START_REF] Fodjo | Separation of concerns in extended epidemiological compartmental models[END_REF] to the mass action and the approximation of the EBCM configuration model of [START_REF] Miller | Edgebased compartmental modelling for infectious disease spread[END_REF] and we were able to obtain similar results close to those of [START_REF] Miller | Edgebased compartmental modelling for infectious disease spread[END_REF]. Both models (mass action model and the approximation of the EBCM configuration model) were defined as a very simple and distinct concern of the basic SIR model in Kendrick.

If it is obvious that we were able to obtain similar results from the EBCM configuration model from the approximation of the EBCM configuration model in Kendrick, it is still a bit limiting. This is because it requires to have beforehand a mathematical formulation of R 0 of each EBCM model that we wish to implement as a concern in the Kendrick sense and to be able to find the appropriate decomposition of the extension points of the relation 11. However, the approach we propose avoids having to build an explicit contact network, and moreover we do not have to define a multitude of differential equations. Another encouraging aspect is the fact that most of the EBCM models proposed in the literature have a mathematical formulation of R 0 , which gives us a wide range of models to explore for future work. In the same way, we also think that in the future, it would be interesting to be able to implement typical EBCM models (without model approximation approach) as concerns that could be composed when adding new concerns.

  orcid.org/0000-0003-4414-7453 1 This is the rate at which susceptible individuals become infected.

Figure 1 :

 1 Figure 1: Flow diagram of the mass action mathematical model    S ′ = -βIS I ′ = βIS -γI R ′ = γI

8

  Figure 3: Definition of the classical concern SIR where λ and the extension points are defined in lines 9 to 10.

Figure 4 :

 4 Figure 4: Definition of the concern of mass action in lines 16 to 22. Then composition of the classic concern and the mass action concern in line 24.

Figure 5 :

 5 Figure 5: Composition of the basic SIR concern and the approximation concern of the EBCM configuration model on a Bimodal degree distribution.

Figure 6 :

 6 Figure 6: Definition of the concern of the approximation of the EBCM configuration model on a Bimodal degree distribution.

Figure 7 :

 7 Figure 7: Initialization of simulation parameters.

Figure 8 :

 8 Figure 8: Simulation and visualization of the configuration model.

  (a) EoN simulations. (b) Kendrick simulations.

Figure 9 :

 9 EoN and Kendrick simulations of the mass action model (black curve) and the EBCM configuration model on four different degree distributions : Homogeneous (blue curve), Poisson (red curve), Bimodal (green curve) and Truncated Power Law (orange curve). Each contact network has 500,000 nodes and an average degree of 5. The Poisson degree distribution has an average of 5, half of the nodes have a degree of 2 and the other half a degree of 8 for the bimodal distribution, and finally a truncated Powerlaw distribution in which P(k) ∝ k -v e -k/40 where v = 1.418.

The notion of partnership, materializes here the existence of a contact (edge), between two individuals. In this case, we say that these individuals are partners.

It is the average number of individuals with whom a susceptible individual is in contact or the average degree of nodes in a contact network.

This is the proportion of contacts that are infectious.

The rate per contact at which disease is transmitted between an infectious individual and a susceptible individual.

This parameter was an approximation of the different contact networks considered.

https://epidemicsonnetworks.readthedocs. io/en/latest/GettingStarted.html

It's model where contact rates in the population are assigned using a P(k) or ρ(k) distribution, while considering that the contact duration is negligible.

 13 It's model in which each node in the dynamic network has a constant-valued degreed, assigned using the P(k) distribution.