Véronique Cova
email: veronique.cova@iae-aix.com

Diego Rinallo
email: diego.rinallo@kedgebs.com

Revisiting the separation between sacred and profane: Boundary-work in pilgrimage experiences

. In this paper, we contribute to such literature with some epistemological reflections on the complex relationships between religions and markets, which we support theoretically and with empirical illustrations from an ongoing study of two Christian pilgrimages in Europe (Santiago de Compostela and Lourdes). Our reflection is based on the concept of boundary-work which brings an alternative perspective to the sacred/profane e opposition.

Introduction

The concept of boundaries is playing an increasingly important role in a wide range of literatures (Lamont 1 ). In social science, it is useful to understand several social relations as the study of nationalism, citizenship, immigration, ethnicity and race, gender and sexual differences. We have referenced two major papers presenting a literature review and a research agenda. Lamont and Molnar (2002) synthesized disparate strains of social science research by analyzing the interaction between symbolic and social boundaries and discussing the many benefits of focusing on fundamental processes of boundary work across these areas. [START_REF] Pachucki | Boundary processes: Recent theoretical developments and new contributions[END_REF] paper takes stock of the most recent scholarship on symbolic boundaries and shows how these interact with social boundaries.

The concept of boundary work was proposed originally by [START_REF] Gieryn | Boundary-work and the demarcation of science from non-science: strains and interests in professional ideologies of scientists[END_REF] who attempts to delineate between what is "science" and what is "non-science". [START_REF] Gieryn | Boundary-work and the demarcation of science from non-science: strains and interests in professional ideologies of scientists[END_REF] understands the practices of scientists as strategic actions to (re)claim scientific authority and jurisdiction over a contested field. He demonstrated how involved actors politically manage the symbolic and social borders of science.

In the same vein, this paper focuses on the problem of demarcation between the sacred and the profaned. We want to show how the boundary between sacred and profane can be managed in one hand, by the institutions and the market, and in the other hand, by the individuals. Boundaries are to be found in interaction between people who identify themselves collectively in different ways; this interaction can in principle occur anywhere or in any context. « Identity is a matter of boundary processes rather than boundaries » (Jenkins, 1996, p.127). Concerning The Way, individuals create symbolic boundaries for separating pilgrims from hickers and for generating feelings of similarity and group membership.

In the same vein, [START_REF] Gieryn | Boundary-work and the demarcation of science from non-science: strains and interests in professional ideologies of scientists[END_REF] analyzes the way boundaries between "science" and "nonscience" are constructed and defended, we identified how similar strategies can be deployed to create or to bolster support for religious approaches. The aim for religious institutions could be to fund religious beliefs and to nuture faith but in the same time, to take opportunity of the marketisation of religion [START_REF] Mittelstaedt | Pilgrimage as Exchange: The Economy of Salvation and the Theory of Exchange[END_REF]. For their part, the pilgrim, as a consumer [START_REF] Belk | The Sacred and the Profane in Consumer Behavior: Theodicy on the Odyssey[END_REF], could aim to live a sacralized experience or to perform it.

The sacred and profane separation

The 'received view', stemming from the work of Émile [START_REF] Durkheim | The Elementary Forms of the Religious Life[END_REF] and Mircea [START_REF] Eliade | The Sacred and Profane: The Nature of Religion[END_REF], is one of complete separation between sacred and profane: religion is about 'sacred things' which are set apart and protected from profaning, mundane forces such as those represented by commerce. According to the separation metaphor, religion and markets are two mutually exclusive spheres of human life. In this paper, we propose that such dualistic, and in a certain sense romantic view of religion separated from the market is just one of many possible discourses about the relationship between the sacred and profane. As also suggested by [START_REF] Belk | The Sacred and the Profane in Consumer Behavior: Theodicy on the Odyssey[END_REF], with the secularization of society, individuals are increasingly finding the sacred in the spheres of politics, science, art, and consumption. Consumer researchers have also frequently found signs of the sacred in various nature or community-based mundane consumption activities. It does not come as a surprise that consumers can ascribe sacred properties to mundane (i.e., non religious or spiritual) ways of consuming the pilgrimage.

Non-dualistic discourses about the sacred have however appeared in interpretive consumer research. O'Guinn and Belk (1989) were perhaps the first to highlight the merger of religion and consumption in the context of a religious theme park. Bilgin (2012) examined how religion and ideology intertwine in the marketplace, influencing in turn the meaning of global brands among religious consumers. [START_REF] Mcalexander | The Marketization of Religion: Field, Capital, and Consumer Identity[END_REF] proposed that religion and consumption are completely interwoven and not at all distinct in the context of the Church of Latter Days Saints. We contribute to such literature by resituating the separation of sacred and profane and suggesting that religious institutions, together with marketplace actors and consumers, at least in some circumstances engage in boundary work to keep the sacred separated from the profane -a task which is however difficult to achieve because of abundant evidence to the contrary. The examination of the discourses and practices that are at the base of such boundary work contributes to a better understanding of the sacred and profane in consumer culture, that is one of the key debates in interpretive consumer research.

Pilgrimage as a research site

We support our theoretical claims with introspections and selected observations from an ongoing ethnographic study on the intersection of market and religion in contemporary Christian pilgrimages. Our main research site is the Way of Saint James (or Camino de Santiago), which is a pilgrimage route to the shrine of the apostle St. James in the Cathedral of Santiago de Compostela in Northwestern Spain, where the remains of the saint are believed to be buried. The first author has conducted participant fieldwork walking as pilgrim the Santiago route for two to three weeks per year from 2010 to 2014. Our secondary sites is the Lourdes sanctuary, in Southwestern France, at a site where the Virgin Mary is believed to have appeared in 1858 to a local peasant girl, Bernardette Soubirous. The second author and a research assistant conducted participant observation in Lourdes for a combined eight days in 2013. Both authors kept detailed field notes and built a database of pictures documenting various aspects of their experiences. Observation data was supplemented with interviews to pilgrims and netnographic analyses of pilgrim blogs ??? and the web sites maintained by religious authorities and marketplace actors.

Our analysis seeks to understand how marketplace dynamics have transformed the age-old institution of pilgrimage, which predates the advent of capitalism, mass tourism, and consumer culture. While focused on the present day, we also immersed ourselves in the history of pilgrimage in general, and the two sites in particular. Early in our analysis, we realized that the vision of a rigid separation between sacred and profane is not supported by historical sources: countless examples can be found of the intertwining of religious and commercial activity in pilgrimages routes already in the Middle Ages (see for example Bell and Dale, 2011). Despite such evidence, pilgrim discourse, popular media, and academic writings often complain that consumerism and commercial forces are desecrating forces undermining the 'true', 'authentic', 'sacred' nature of pilgrimage. Also scholars have often diffused such romantic view, displaying a tendency to demarcate what pertains to religion in a strict sense (the sacred) from mundane aspects of existence, including those pertaining to commerce. This is also evident in theoretical attempts to distinguish 'true pilgrims' from tourists (e.g., [START_REF] Cohen | Pilgrimage and Tourism: Convergence and Divergence[END_REF]. More recent research has however observed that commerce and pilgrimage are intertwined. In perhaps one of the best treatment of the subject, [START_REF] Reader | Pilgrimage in the Marketplace[END_REF] propose that "the dynamics of the marketplace, with its themes of pilgrimages being promoted, reshaped, invented and exhibited to increase their custom, along with issues of consumerism and the acquisition of material goods and souvenirs, are not antithetical to pilgrimage (or to 'religion'), but crucial to its successful functioning, development, appeal and nature" (p. 15). Still, the separation of sacred and profane remains a powerful narrative affecting consumer experience and the promotion of pilgrimage sites by both religious and lay actors.

The sacred and profane in pilgrimage

In this section, we provide some examples of how agents engage in boundary work to keep religion and market separated in ways that are sometimes peaceful and uncontested, sometimes conflict-laden. This is perhaps best exemplified by Lourdes. Outside the Sanctuary, hundreds of merchants sell religious objects (e.g., rosaries, statues of the Virgin Mary, candles). Inside, in contrast, there is no apparent sign of commercial activity. The candles that are sold by outside vendors, inside the shrine are available in exchange of a free offer (not a price -even though 'suggested' offers, depending on candle size, are listed in various currencies). More importantly, the water from the spring that emerged in connection to the Marian appearances is freely available to everybody at any time of day and night. Outside, vendors sell containers of various sizes (often in plastic: glasses, virgin-shaped bottles, flasks, etc.) but the water per se is not sold. For those who want to receive the water of Lourdes at home, the Sanctuary offers a shipping service -also in this case, however, the water is free: "The applicant only pays the postage and packing". The rigid boundary between the Sanctuary -free from commercial activity -and the outside world, regulated by a commercial logic, does not however occur naturally: it is created, maintained and enforced by Sancturary authorities that maintain a 'vigilante' approach and act swiftly in case of transgressions (see below) In other words, Sanctuary authorities prevents the sacred from being profaned by marketplace forces and actors at the discursive and material level.

In the context of Santiago de Compostela, marketplace actors, working in tandem with tourism development authorities, re-packaged pilgrimage as heritage tourism, downplaying its religious significance in order to make appeal to a wider audience. Despite the traditional emphasis on pilgrim devotion and contemplation, which is still evident in the hardships suffered by those who walk the route by feet, burdened by backpacks, sleeping in Spartan communal lodging, and avoiding alcohol, and frequently attending Mass, many tour operators now sell hyper-real, glamorous version of the Camino experience (luxury hotels, fine dining, bags sent from one stop to the other, etc.). The Way of Saint James today attract a very diverse group of consumers for a variety of reasons -often attributing limited importance to religion or spirituality. Still, some consumers describe their experience as "true pilgrimage" and, despite being immersed in a market environment, do of their best to avoid what they perceive as inauthentic, watereddown, profane version of the Camino experience. Their sacralizing boundary work often entails resentment and criticism of those consumer who are tempted by the market and "do not understand the true nature" of the pilgrimage. Also marketers are very frequently criticized for being animated by excessively commercial goals and often treating 'true pilgrims', who are less profitable than other market segments, disrespectfully.

Boundaries-work and actors strategies

Religious institutions are caught between two positions: An activist approach to evangelize and make converts to the religion (to run after the lost sheep) or a marketisation of their practices. One of the sign of the Catholic Church's declined religious authority is the re-signification of the sacred nature of pilgrimage by marketplace actors and consumer alike. Some consumers, while still ascribing sacred meanings to pilgrimage, do so in ways that are very far from the Catholic Church's "official" teachings. As suggested by critics of the Turners' influential work [START_REF] Turner | Image and Pilgrimage in Christian Culture[END_REF], pilgrimage is indeed a realm of competing discourses (Eade and Sallnow, 1991), and those emanating by the Catholic Church are not necessarily the most relevant for the actors envolved.

The Santiago route is often seen as a pathway to self-transformation; a rite of passage; a way to get into contact with nature or the authentic self; and an experience brotherhood/sisterhood or community thanks to authentic encounters with fellow-minded individuals. This way of experiencing the sacred (which, for lack of a better alternative, we might term the 'mundane' sacred) is just one of many. The Santiago route can also be experienced as sacred based on religious or spiritual meaning systems distinct from those supported by the Catholic Church. For example, best-selling author Paulo Coelho's (1987) autobiographical book The Pilgrimage (also known as Diary of a Mage) describes the route in terms of its significance for the Western mystery traditions and esoteric Christianity, whereas some of our informants undertake the route with Neopagan or New Age sensibilities. For example, in the movie The Way, Tom (Martin Sheen) honors his son's desire who died on the Camino, by finishing the journey. Inexperienced as a trekker, Tom soon discovers that he will not be alone on the Way. He meets other pilgrims from around the world, each with their own issues and looking for greater meaning in their lives. Through Tom's unresolved relationship with his son, he discovers the difference between the life we live and the life we choose.

In contrast, in Lourdes many consumers experience the sacred in ways mostly consistent with the doctrine of the Catholic Church. Some of our atheist and agnostic informants told us "you can feel there's something there". Others, raised as Catholics but turned non-believers, suggested that the elements of the Church rituals brought them childhood memories and back to a time where they 'used to believe', as long-forgotten words and gesture were spontaneously, albeit mechanically, repeated. At the many shops surrounding the Sanctuary, only very occasionally we observed the sale of merchandise not aligned with official Church teachings: in one case, it was crystals, thought to bring luck or have healing properties according to various New Age teachings; in an other, a statue of the Archangel Uriel who, while still very important in contemporary occultism, unlike Michael, Gabriel and Raphael is not mentioned in the Bible and, as a result, is not considered eligible for veneration.

It is interesting to mention that on the one hand, whereas Catholic Church states clearly that it must exist a formal separation between the sacred and the profaned (to fit with the Holy Scriptures), it allows market organizations around its religious practices. For it, the boundary between the sacred and the profane is blurred, leaving space for marketisation. It needs the market. On the other hand, pilgrims deplore that religion increasingly has introduced profaned activities and call for the resacralization of the religion. This double opposite movement is at play in each f of our research fields as the below The pilgrims who strive for less secularity

The institutions which ghettoize the shrine

Protecting the sacred from market profanations

We observed limited Church reactions to pilgrimage marketization in the case of Santiago. Those in charge of running pilgrim hospices routinely ask questions to ascertain whether those who demand lodging are 'true pilgrims' or simply profane tourists who want to exploit the low prices of these facilities. On the contrary, in Lourdes -where the city has grown exponentially thanks to the Marian appearances -religious institutions still hold significant power and have been able to exert their authority on marketplace actors. For example, we were told by hotel managers that in the 1960s it was not possible to equip rooms with TV sets, because of a concern that television might distract customers from a proper pilgrimage experience. Similarly, hawkers are strictly forbidden inside the Sanctuary.

In Lourdes, we also observed an interesting case of marketplace appropriation and resignification. A commercial operator created in 2009 an organic cosmetics brand (Crème de Lourdes) that employs the spring water from the grotto of Lourdes. Four years later, shrine authorities took notice and, while not bringing the case to courts (the claim 'water from Lourdes' would apply to tap water from Lourdes, making it unlikely to win the battle), publicly condemned the commercial exploitation of the water and warned the brand owner to remove from his website all pictures and videos from, and links to, the official shrine website. Additionally, the Church denies that the grotto water has any healing properties; the many miracles experienced by the faithful in Lourdes are, based on Catholic theology, to be ascribed to Divine intervention. Without an official religious sanction, Crème de Lourdes has to justify the water's healing properties (and the premium price applied to its products) to a competing meaning system, based on New Age views and unorthodox science on the electromagnetic properties of water.

Conclusion

Pilgrimage exists in symbiosis with the marketplace. We are not the first to highlight this and to raise the more general point that the sacred and profane are and have always been interwoven. However, boundaries between them fluctuate differently according to the cases. In Lourdes, the Catholic Church has created a strong border. It drives away « the merchants from the Temple » and separates the sacred from the secular. In Santiago de Compostela, the hirelings have invaded the religious shrine and quite all places ; all activities are marketed. So, in the first case, Lourdes, the sacred is limited in a dedicated space, on specific objects (the water), in a closed organization. In the second case, Santiago, the sacred loses its central position and becomes transformed in profane services (to pay for Confession or for stamp on Credential), in marketed products (religious trinkets), in commercial places (La Cruz de Ferro de Foncebadon for touristic destination). On the other side, the consumer either accepts and even contributes to the marketisation of the holy pilgrimage, or he/she wants to perceive the sacred as separate from the profane and actively aspires to this split (Santiago). Consumers, as religious institutions and marketplace actors alike similarly engage in boundary work.

Both Lourdes and Santiago de Compostela are (to a varying extent) dominated by the Roman Catholic Church; had we observed Jerusalem, which is a holy city for the three Abrahamic religions, or Glastonbury [START_REF] Maclaran | Magic and Merchandise: Spiritual Shopping in Glastonbury[END_REF], which has special significance for Christianity, Neopaganism and the New Age, such variety would have been more evident. By "opening the black box" of the sacred/profane boundary, we hope to advance literature on the marketing and consumption of religion and spirituality.

table shows

 shows

		:
		Santiago (The Way)		Lourdes
	Religious institutions…	develop a marketisation of the	advocate	for	a	clear
		pilgrimage	separation from the profane
	Pilgrims …	take offence and regret a lost of
		spirituality
	Who enact the boundary-
	work ?

http://educ.jmu.edu//~brysonbp/symbound/papers2001/LamontEncyclo.html


François Bonneau

Guillaume Caumon

Philippe Renard

Impact of a Stochastic Sequential Initiation of Fractures on the Spatial Correlations and Connectivity of Discrete Fracture Networks

Stochastic discrete fracture networks (DFNs) are classically simulated using stochastic point processes which neglect mechanical interactions between fractures and yield a low spatial correlation in a network. We propose a sequential parent-daughter Poisson point process that organizes fracture objects according to mechanical interactions while honoring statistical characterization data. The hierarchical organization of the resulting DFNs has been investigated in 3D by computing their correlation dimension. Sensitivity analysis on the input simulation parameters shows that various degrees of spatial correlation emerge from this process. A large number of realizations has been performed in order to statistically validate the method. The connectivity of these correlated fracture networks has been investigated at several scales and compared to those described in the literature. Our study quantitatively confirms that spatial correlations can affect the percolation threshold and the connectivity at a particular scale.

Introduction

Fracture Networks in rock masses significantly affect hydrodynamic and geomechanical processes. Due to incomplete observations and measurements at all relevant scales, numerical modeling is widely used to describe the geometry and understand the behavior of fractured rocks. However, efficiently honoring both statistical descriptions of the fractured medium and all the physical principles driving fracture genesis remains a challenge. This paper reports on investigations to integrate physical and geological principles in stochastic discrete fracture network (DFN) simulations and assessing the impact of this integration on emerging DFN properties such as network connectivity and fractal dimension. Spatial characteristics of subsurface fracture networks are the result of complex processes. These processes have been described by linear elastic fracture mechanics and studies about the initiation, growth, interaction and termination of fractures [START_REF] Olson | Joint pattern development : effects of subcritical crack growth and mechanical crack interaction[END_REF][START_REF] Renshaw | Numerical simulation of fracture set formation: a fracture mechanics model consistent with experimental observations[END_REF][START_REF] Tuckwell | The control of stress history and flaw distribution on the evolution of polygonal fracture networks[END_REF][START_REF] Jing | A review of techniques, advances and outstanding issues in numerical modelling for rock mechanics and rock engineering[END_REF][START_REF] Welch | A dynanic model for fault nucleation and propagation in mechanically layered section[END_REF]. Fractures initiate at rock flaws [Pollard andAydin, 1988, Cosgrove and[START_REF] Cosgrove | The initiation, propagation, and arrest of joints and other fractures[END_REF]. Heterogeneities such as fossils, grains, cavities, micro-cracks and other objects having elastic properties different from those of the surrounding rock modify the stress field in such a way that the magnitude of local stresses may exceed the strength of the rock, thereby initiating a fracture. Fractures create local heterogeneities in both the stress field and the bulk properties of the rock. [START_REF] Pollard | Progress in understanding jointing over the past century[END_REF] describe the decrease of mechanical stress along cracks (in the shadow zone) and its increase at fracture tips (in the stress concentration zone). Fractures then grow as long as the effective stress locally exceeds the strength of the rock [START_REF] Griffith | The phenomena of rupture and flow in solids, Philosophical transactions of the royal society of london[END_REF][START_REF] Griffith | The theory of rupture[END_REF]. Numerical modeling of this process requires to evaluate the stress state everywhere and at all times in the simulation domain. It is generally based on boundary conditions and uses finite element, boundary element or discrete element methods which are computationally intensive for large scale 3D problems and are often subject to gridding and discretization problems. Therefore, stochastic approaches are often applied for efficiently simulating very large fracture networks (for reviews, see e.g. [START_REF] Jing | A review of techniques, advances and outstanding issues in numerical modelling for rock mechanics and rock engineering[END_REF], [START_REF] Chilès | Stochastic modeling of natural fractured media: A review[END_REF], [START_REF] Dershowitz | Advances in discrete fracture network modeling[END_REF], [START_REF] Dowd | A comparison of methods for the stochastic simulation of rock fractures[END_REF]; see also De [START_REF] De Dreuzy | Synthetic benchmark for modeling flow in 3D fractured media[END_REF] for a benchmark). Stochastic approaches do not solve the mechanical problem but aim at defining a mathematical proxy that reproduces the final geometry of natural fracture networks. Those methods rely on the statistical inference and characterization of the fracture network from field and analog observations. Classical stochastic discrete fracture network simulations consider fractures as simple 2D objects (rectangles or ellipses) in 3D space. The position of rock flaws stimulated by the stress field is usually generated either by a homogeneous or heterogeneous Poisson point process [START_REF] Stoyan | Fractals, random shapes and point fields[END_REF]Stoyan, 1994, Lantuéjoul, 2002]. The dimension and orientation of the fracture are simulated for each object independently by sampling from prior probability distribution functions with a Monte Carlo process. This approach generally simulates fractures set by set following their chronology in relation to tectonic events (e.g., [START_REF] Mace | Integration of fracture data into 3D geomechanical modeling to enhance fractured reservoirs characterization[END_REF], [START_REF] Yamaji | Clustering of fracture orientations using a mixed Bingham distribution and its application to paleostress analysis from dike or vein orientations[END_REF], [START_REF] Lamarche | Relationships between fracture patterns, geodynamics and mechanical stratigraphy in carbonates (South-East Basin, France)[END_REF]). However, the geometry and the position of each fracture in each set is randomly and independently selected. Even if the statistics describing each set are honored, the flow behavior and emerging fracture network features such as connectivity and spatial correlations may be biased. This is attested for instance by [START_REF] Odling | A "conductance" mesh approach to the permeability of natural and simulated fracture patterns[END_REF] who highlight in two dimensions the differences of global permeability between natural and simulated fracture networks sharing the same length and orientation statistics. This problem can be addressed by several improvements to Poisson-based DFN simulation methods. An avenue is to refine the statistical models used in the DFN simulation see for instance [START_REF] Chilès | Fractal and geostatistical methods for modeling of a fracture network[END_REF], [START_REF] Xu | A new computer code for discrete fracture network modelling[END_REF]. Another is to incorporate the effects of mechanical interactions which control fracture initiation, growth and arrest. This has been done by reproducing the effects of these interactions on the fracture geometry [START_REF] Josnin | A three-dimensional model to simulate joint networks in layered rocks[END_REF], and mimicking the fracture growth using planar fractures [START_REF] Cladouhos | Are fault growth and linkage models consistent with power-law distributions of fault lengths?[END_REF][START_REF] Swaby | An Interactive 3d Fracture Modelling Environment[END_REF][START_REF] Davy | A model of fracture nucleation, growth and arrest, and consequences for fracture density and scaling[END_REF] or non planar fractures [START_REF] Cacas | Nested geological modelling of naturally fractured reservoirs[END_REF][START_REF] Srivastava | Geostatistical simulation of fracture networks[END_REF][START_REF] Bonneau | A methodology for pseudo-genetic stochastic modeling of discrete fracture networks[END_REF]. Among these methods, [START_REF] Cladouhos | Are fault growth and linkage models consistent with power-law distributions of fault lengths?[END_REF] relate a simplified two-dimensional linkage model to the emergence of power law distributions. Similarly, [START_REF] Davy | A model of fracture nucleation, growth and arrest, and consequences for fracture density and scaling[END_REF] establish theoretical and numerical links between fracture seeding and development processes and emerging fracture statistics in simulated DFNs (density and power law length distribution). In this paper, we use similar ideas to evaluate numerically the link between pseudo-mechanical DFN simulation and emerging characteristics such as the connectivity and fractal dimension. For this, we propose a simulation method which emulates the stress perturbation effects around fractures during their seeding. This numerical method aims at both reproducing the statistics describing natural fracture networks and the geometrical hierarchy between fractures. The main idea is to use a sequential parent daughter Poisson point process that favors the initiation of small fractures at the tips of larger fractures (Section 2). A large number of realizations is performed on synthetic data to validate the method in a statistical sense. We discuss and quantify the impact of this sequential fracture seeding on the DFNs organization in terms of spatial correlation (Section 3), connectivity and percolation (Section 4).

Sequential Stochastic Discrete Fracture Network Simulation

Several pseudo-genetic approaches have been proposed to simulate fracture geometry by propagating fractures from an initial state using simple mechanical concepts on growth and arrest [START_REF] Josnin | A three-dimensional model to simulate joint networks in layered rocks[END_REF][START_REF] Srivastava | Geostatistical simulation of fracture networks[END_REF][START_REF] Bonneau | A methodology for pseudo-genetic stochastic modeling of discrete fracture networks[END_REF][START_REF] Davy | A model of fracture nucleation, growth and arrest, and consequences for fracture density and scaling[END_REF]. Instead, this work focuses on the fracture nucleation to capture hierarchical organization of natural fracture networks. We are making the assumption that fracture radial growth and coalescence is correctly described by the length distribution law. Objects are progressively implanted in their final state and so we neglect the fracture growth, interaction and coalescence that leads to the final fracture geometry. Large fractures are simulated first because they have grown far from discontinuities that may interact and stop their growth. We use an approximate local stress perturbation model (Figure 1, Section 2.1) to later organize small fractures around these major objects with a probabilistic acceptance of fracture seeds (Section 2.2).

Simplified Fracture Object

Produced DFNs are made of planar objects characterized by their position, extension and orientation. These objects may represent either a single or several coalesced discontinuities. , 1996[START_REF] Srivastava | Geostatistical simulation of fracture networks[END_REF][START_REF] Bonneau | A methodology for pseudo-genetic stochastic modeling of discrete fracture networks[END_REF]. Most of them used a spherical isotropic model. In this work, we use a simple model constituted by two ellipsoids in order to take into account the anisotropy often observed in stress concentration around fracture tips (Figure 1, [START_REF] Lyakhovsky | Scaling of fracture length and distributed damage[END_REF]).

Marrett

These ellipsoids separate volumes where the probability to observe a crack is increased (respectively decreased) because of the stress concentration (respectively relaxed). The geometry of each ellipsoid is defined by factors weighting the fracture size ; (α s , β s ) for shadow zone and (α p , β p ) for stress concentration zone. This model, inspired by the shape of symmetric stress perturbation around tensile fractures, aims at approximating the first order stress perturbation of all types of fractures. This may be locally inconsistent (e.g., if fractures are parallel to the main compressive stress direction). However far-field stress is generally not known during DFN simulation, so this model is expected to emulate geomechanical effects on average (Section 5).

The next Section describes a stochastic simulation process that uses these ellipsoids to drive the seeding of fractures.

Sequential Parent-Daughter Poisson Point Process

Our global sequential seeding process is described in Appendix A and Figure 2. It uses a heterogeneous Poisson point process that generates N potential fracture seeds according to a prior fracture density map (d (x, y,z) ) that may be obtained by strain analysis [START_REF] Cacas | Nested geological modelling of naturally fractured reservoirs[END_REF][START_REF] Mace | Integration of fracture data into 3D geomechanical modeling to enhance fractured reservoirs characterization[END_REF], seismic methods [START_REF] Sayers | Seismic characterization of reservoirs containing multiple fracture sets[END_REF][START_REF] Serrano | Distribution of crack density parameter in Central Betic Cordillera (Southern Spain)[END_REF][START_REF] Rodriguez-Herrera | Stress-induced signatures in 4d seismic data: Evidence of overburden stress arching[END_REF] or other ancillary observations. Then, the nucleation of fractures uses a probabilistic acceptance/rejection test that takes into account the fractures previously simulated. Simulation is then performed in a given number of sequences (s ∈ [1, N ]). During the first seeding sequence, the probability (P (x, y,z) ) to simulate a fracture at a position generated by the Poisson point process is uniform and equal to: N s 1 N = 1/s. Planar fractures are simulated by a Monte Carlo sampling of orientation and dimension distribution laws. Around each simulated fracture, a simplified geomechanical model (Figure 1) divides the volume of interest into several parts in which the prior probability (P (x, y,z) ) to activate a fracture at a given seed (x, y, z) has to be locally adjusted considering:

1. The number of stress concentration zones overlapping the (x, y, z) location (N acc (x, y, z)) and its associated effect g: the probability to stimulate a rock flaw inside stress concentration zones is increased by a factor g ≥ 1.

2. The number of shadow zones defined at the (x, y, z) location (N shad (x, y, z)) and its associated effect h: the conditional probability to stimulate a rock flaw inside the shadow zone is decreased by a factor h ≥ 1.

h and g are input parameters standing for the impact of shadow and stress concentration zones. As a consequence, we define the adjusted probability P s (x, y,z) at each step of a process divided in s sequences as follows:

P s (x, y,z) = N s 1 + g.N acc (x, y, z) 1 + h.N shad (x, y, z) / N i=1 1 + g.N acc (x, y, z) i 1 + h.N shad (x, y, z) i
(1) For each seeding sequence, the Poisson point process simulates N seeds that are scanned twice (algorithm in Appendix A). The first iteration evaluates the denominator of Equation ( 1) by counting N acc (weighted by g) and N shad (weighted by h) for each seed. The second iteration uses the previous result that quantifies the fracturing state to compute and normalize P (x, y,z) according to the impact of already simulated fractures.

The proportion of stress concentration/relaxation zones (Figure 1) depends on both the fracture length and on the number of simulated fractures. In other words, the proportion of these zones depends on the fracture density. In the early stage, only a few fractures have been activated and the fracturing process is not strongly impacted by previous fractures. In this case, every fracture center proposed by the Poisson point process has the same probability (P (x, y,z) = 1/s) to be selected for simulation because N acc (x, y, z) ≈ N shad (x, y, z) ≈ 0. During the simulation, the number of fractures increases and consequently their impact on the simulation process increases as well. Finally, the method naturally reproduces the "poorly" and "well-developed" fracture networks [START_REF] Wu | An experimental study of the relationship between joint spacing and layer thickness[END_REF] and the transition regime between these two states [START_REF] Davy | A likely universal model of fracture scaling and its consequence for crustal hydromechanics[END_REF]. "Poorly" and "well-developed" fracture networks have also been described in terms of "unsaturated" and "saturated" [START_REF] Josnin | A three-dimensional model to simulate joint networks in layered rocks[END_REF] or "dilute" and "dense" fracture networks [START_REF] Davy | A likely universal model of fracture scaling and its consequence for crustal hydromechanics[END_REF].

This approach organizes small fracture objects at the tips of longest ones. It neglects the interaction between fractures during their growth and the fractures that are nucleated during the fracture growth. As a result, our method may underestimate the number of fractures that initiate from older fractures, which may lead to an underestimation of the network connectivity. Strategies to incorporate the effects of interactions during fracture growth should therefore be considered in future research.

Impact on DFN Spatial Correlation

In this section, we analyze the DFN spatial correlations emerging from the above sequential seeding process. A quantitative characterization of spatial correlations is also possible using fractal theory. This aspect has been much studied in the literature for its ability to describe natural fracture systems [START_REF] Mandelbrot | The fractal geometry of nature[END_REF][START_REF] Velde | Fractal patterns of fractures in granites[END_REF][START_REF] Ouillon | Hierarchical geometry of faulting[END_REF][START_REF] Bour | Clustering and size distributions of fault patterns: theory and measurements[END_REF][START_REF] Bour | A statistical scaling model for fracture network geometry, with validation on a multiscale mapping of a joint network (hornelen basin, norway)[END_REF], Darcel et al., 2003]. The measurement of the fractal dimension of fracture networks can be related to several definitions that quantify the singularities of the multifractal spectrum. [START_REF] Bonnet | Scaling of fracture systems in geological media[END_REF] present a review for characterizing and estimating the capacity dimension, the information dimension and the correlation dimension that can be considered as the three first moments of the multifractal spectrum.

In this work, we quantify the correlation dimension D c of fracture network models using a pair correlation function (C 2 (r)) computed from the fracture center repartition, similarly to [START_REF] Davy | Some consequences of a proposed fractal nature of continental faulting[END_REF] and [START_REF] Bour | A statistical scaling model for fracture network geometry, with validation on a multiscale mapping of a joint network (hornelen basin, norway)[END_REF] :

C 2 (r) = 2 × N (r) N t × (N t -1) ≈ r D c (2)
where N (r) is the number of pairs of fracture centers whose separation distance is less than r. N t is the total number of fracture centers in the considered volume of rock. [START_REF] Ouillon | Hierarchical geometry of faulting[END_REF] and [START_REF] Bonnet | Scaling of fracture systems in geological media[END_REF] have shown that computing the correlation dimension using Equation ( 2) provides more accurate spatial correlation estimates than the classical boxcounting method. This method assumes an isotropic distribution of fracture centers which is a reasonable approximation in simple cases where we start from intact rock and uniform orientation distributions. We compare the spatial correlations that emerge from a hundred DFNs simulated with either our sequential Poisson point process or the homogeneous Poisson process. To get a consistent correlation dimension, we simulate DFNs made of fractures whose lengths range over two orders of magnitude (1 to 30 meters). Consistently with our hypothesis that the fracture length distribution should reflect the growth and the coalescence process, we choose a power law length distribution of exponent 1.8 [START_REF] De Dreuzy | Hydraulic properties of twodimensional random fracture networks following a power law length distribution: 2. permeability of networks based on lognormal distribution of apertures[END_REF], truncated between 1 and 30 m. Each DFN contains 10, 000 fractures simulated from a uniform density map in a volume of dimensions 90×90×90 m. The correlation dimension estimation focuses on the central volume (30 × 30 × 30 m) to reduce edge effects due to fracture truncation. For the sequential Poisson point process, 0.1% of the total number of fractures was simulated during each sequence. Other parameter values for our pseudo-geomechanical model (Figure 1) were arbitrarily chosen as follows:

β s = 1, α s = 1, h = 100; β p = 0.5, α p = 2,
g = 100 (see Section 3 for a sensitivity analysis related to these parameters).

We studied the impact of the sequential parent-daughter Poisson point process on DFN correlation dimension for three kinds of DFNs:

• DFNs made of one set of vertical fractures (strike uniformly distributed between 80 • and 100 • ) Figure 3 • DFNs with fractures that can have any orientation (fracture dip are uniformly distributed between 0 • and 90 • ; and strikes uniformly distributed between 0 • and 360 • ) Figure 3 (c) and Figure 4 (c).

The homogeneous Poisson point process does not generate spatial correlations (Figure 3). The average fractal dimension computed over a hundred DFN realizations is very close to 3 in all types of DFN (Table 1). On the contrary, the sequential seeding of fractures (Figure 4), creates correlated DFNs with a correlation dimension between 2.47 ± 0.09 for DFNs made with a unique fracture set and 2.58 ± 0.09 for DFNs made of randomly oriented fractures (Table 1). This demonstrates that defining simple rules, inspired from geomechanical concepts, to constrain progressive stochastic fracture nucleation impacts the correlation dimension of produced networks. This also suggest that prior knowledge of fracture set orientation slightly affects the correlation dimension. However, the correlation dimension defined Equation 2 is making the assumption an isotropic positioning of fractures. This assumption is difficult to make with our method if fractures are not randomly oriented. In the following, we focus on density and connectivity of randomly oriented fractures and we isolate the effect of the sequential seeding strategy.

Impact of the Sequential Seeding Process on Fracture Density

This section studies the ability of the produced DFNs to honor the input fracture density map that quantifies the number of fracture centers per unit volume (P30). We quantify the average number of fracture centers that are simulated per volume unit by mapping fracture centers on the simulation volume (90 × 90 × 90 Cartesian grid of 1 × 1 × 1 m voxels, Section 2.3). We compute maps of fracture center likelihood (often called E-types, see [START_REF] Journel | Nonparametric estimation of spatial distributions[END_REF]) over a large number of realizations and we compare them to input fracture density maps in order to detect and quantify the mismatch due to the simulation method. We simulate 1,000 DFNs using the parameters described in Table 2. E-types are first computed for uniform fracture density maps and using a homogeneous Poisson point process. Input density is 1.37 × 10 -2 cent er.m -3 , corresponding to DFNs with 10, 000 fractures (Figure 5a). The homogeneous Poisson point process perfectly reproduces the input density map. It produces a uniform E-type map with an average value of 1.37 × 10 -2 cent er.m -3 and a standard deviation of 2.08 × 10 -4 m -3 (Figure 5b). The sequential seeding process also produces an E-type map almost uniform with an average values of 1.37 × 10 -2 cent er.m -3 but with a higher standard deviation (1.30 × 10 -3 m -3 , Figure 5c). The higher standard deviation indicates that the convergence of E-type is more difficult due to the clustering of each realization. The sequential process also introduces a small bias in the fracture repartition by simulating more objects in the center of the volume. We tried to reduce the intensity factors from h = g = 100 to h = g = 10. This increased the DFNs correlation dimension from 2.65 ± 0.09 to 2.78 ± 0.06. The E-type map obtained in this case is still very close to the one presented in Figure 5c. Additional tests on larger Cartesian grids suggest that the acceptance/rejection process produces subtle edge effect, which could be further investigated.

We performed the same study with a heterogeneous density map (Figure 6a). The E-type map corresponding to a heterogeneous Poisson point process (Figure 6b) is correlated to input density. The linear correlation coefficient linking the two properties is equal to 0.98. Note that producing DFNs that follow a heterogeneous spatial repartition does not create significant spatial correlations between objects (D c ≈ 3). The E-Type obtained by mapping correlated DFNs (D c = 2.66) is correlated to the reference map but with a lower coefficient (0.76, Figure 6c). As in the homogeneous case, we note a border effect that biases the global fracture repartition. It seems that objects are more likely to be simulated at the center of the volume. Reducing the intensity factors from h = g = 100 to h = g = 10 increases the DFNs correlation dimension to D c = 2.79, reduces spatial correlations and increases the correlation factor between the reference map and E-type map to 0.85. Because the intensity of spatial correlations may depend on the parameters of the sequential seeding process, we propose a sensitivity analysis in Section 3.

Sensitivity Analysis on the Sequential Poisson Point Process

The parameters of the method are (1) the relative extents of stress concentration zone (α p and β p ), (2) the relative extents of the shadow zone (α s and β s ), (3) the intensity factors (g and 5a) 3.01 ± 0.06 Figure 6b Classical heterogeneous (Figure 6b) 2.98 ± 0.06 Figure 5c Sequential1 homogeneous (Figure 5a) 2.65 ± 0.09 Figure 6c Sequential 1 heterogeneous (Figure 6b) 2.66 ± 0.08 h, Equation (1)) and (4) the number of seeding sequences (s). All these parameters relate to the mechanical parameters of the rock, but no direct or statistical relationships have yet been established. This sensitivity study aims at quantifying the relationships between input parameters and spatial correlations in DFNs.

Simulations are performed with similar parameters as in Section 2.3:

1. Simulations are performed on a Cartesian grid with 90× 90 × 90 cubic voxels.

2. Fracture seeding uses an uniform density map generating 10, 000 fractures.

3. The distribution law describing the fracture length follows a truncated power law with an exponent of 1.8 and values bounded between 1 and 30 meters. In this section, we describe and quantify the impact of both the geomechanical model and the chronological settings on the DFNs spatial correlations. The results of this sensitivity analysis are plotted in Figures 7, 8 and 11 and discussed in the next sections. For each set of tested parameters, we show a 2D section of a 3D DFN (red section, Figure 3 and Figure 4).

For each cases, we associate an average estimation of the correlation dimension computed over 100 realizations using the method described in Section 2.3.

Sensitivity Study on Geomechanical Model Geometry Factors

To test the sensitivity of the method, we simulated 10, 000 fractures in the same domain as in previous tests, by implanting ten fractures per sequence (0.1% of the final DFN). We have varied the extent of stress concentration and shadow zones changing the values of β s and α p in the range [1-3] and [2-4], respectively, while keeping the other input factors constant (β p = 1, α s = 0.5 and intensity factors g = h = 10).

The sensitivity analysis shows that the increase of shadow or stress concentration zones extents increases the correlation dimension of DFNs produced. In spite of the relatively large standard deviation of the correlation dimension, we observe a progressive dilution of spatial correlations from the clustered DFN (Figure 7a) to the uncorrelated DFN (Figure 7i).

Sensitivity Study on the Geomechanical Model Intensity Factors

In this sequence, we use the same chronology as in Section 3.1 (0.1% of the DFN is simulated per sequence). 100 simulations have been run using the geometric parameters described in Figure 1 (α s = β s = 1, α p = 2 and β p = 0.5). The intensity factors in the stress concentration and shadow zones (g and h) vary between 10 and 100.

Results show a progressive decrease of the fractal dimension from D c = 2.79 (Figure 8a) to D c = 2.58 (Figure 8i) when g and h increase. Figure 8 also suggests that the impact of g is not as significant as the impact of h. We note a strong variation in the correlation dimension when we increase the value of h from 10 (D c = 2.8, Figure 8a, b andc) to 50 (D c = 2.6, Figure 8d, e and f). Several other simulations have been made in order to better characterize this trend (Figure 9).

More quantitatively, we found that the correlation dimension D c mainly varies as a function of h (for g = 10) following the experimental model :

D c (h) ≈ -0.18 log h + 2.96 ± 0.09 (3)
The correlation dimension D c , however, is a quasi-constant function of g. Considering h = 10 the experimental model of D c (g) is described by : D c (g) ≈ -0.03 log(g) + 2.81 ± 0.09 (4)

The low impact of the parameter g can be explained because the correlation dimension is a parameter that quantifies the spacing between fracture centers. Because g mainly attracts small fractures at the tip of larger ones, it affects very marginally the distance between fracture centers. Does this finding hold in the case of preferential fracture orientation? To answer that question, we repeated this study in the case of single set DFNs. Preferential fracture orientation does not change the behavior of D c (h) because we used an isotropic shadow zone model. Then, we focused on the D c (g) and established the following experimental model:

D c (g) ≈ -0.05 log(g) + 2.83 ± 0.08 (5)
The similarity between Equations 4 and 5 suggests that the relation between the correlation dimension D c and the intensity factors g and h is marginally affected by the orientation distribution of the fracture set.

Sensitivity Study on the Fracture Seeding Chronology

Theoretically, an idealized simulation should generate fractures one by one to take into account the impact of every newly simulated fracture in further simulation steps. In practice, simulating only one fracture per sequence is time consuming. For example, in our implementation on a Intel(R) Core(TM) i7-2600 CPU @ 3.4GHz, building a DFN with 10,000 fractures takes only few seconds when done in one sequence; but it takes approximately ten minutes when simulating fractures 10 at a time. In our simulation method, the proportion of fractures that are simulated per sequence p varies with the total number of sequences (s) following:

p(s) = 100 s (6)
This means that for high values of s, we need a strong increase of s to observe only a small decrease of p(s). A high increase of s leads to a high increase of the computation time. The current section aims at testing the sensitivity of the method concerning the chronological seeding of fractures.

For this analysis, we simulated 10, 000 fractures with a number of sequences (s) that varies between 1 and 1, 000 (Figure 10). Simulations have been done using the parameters summarized in Table 3.

The results suggest that simulating 10% of fractures at each sequence is sufficient to converge to a constant correlation dimension (D c ). Furthermore, the value of D c stabilizes as soon as 0.4% of the DFN is simulated at each sequence. Another result of the sensitivity analysis is that several DFNs may have the same correlation dimension but different spatial organizations (Figure 11). In other words, the hierarchical organization of small fractures that clustered around largest fracture tips may not be described by D c . Indeed, the correlation dimension only quantifies the clustering of fracture centers and does not quan- and 10 sequences; it means that 0.1%, 1% and 10% respectively of the fracture network have been simulated per sequence. Whereas the correlation dimension is comparable (2.8), we note some differences in the global organization (X-shapes, hierarchy). This impact is quantified in Section 4.

tify the position of fractures according to their size. Figure 11 shows the organization of fracture traces obtained with various number of sequences. An increase of the number of sequences seems to cluster small fractures at large fracture tips and also seems to decrease the number of X-shaped contacts. This may have an impact on the DFNs connectivity that we propose to study in Section 4.

Impact of the Sequential Poisson Point Process on 3D DFNs Connectivity

Connectivity is a key property of fracture networks because it is directly linked to physical behavior of fractured rocks. Numerous studies have used percolation theory to characterize flow and transport in porous and fractured media (see [START_REF] Berkowitz | Percolation theory and its application to groundwater hydrology[END_REF] and [START_REF] Berkowitz | Characterizing flow and transport in fractured geological media : A review[END_REF] for reviews). In this section, we investigate the DFN percolation threshold and local connectivity model that emerge from the sequential Poisson point process. We describe and quantify the impact of spatial correlations and more generally the impact of the chronological seeding of fractures on percolation properties in the isotropic case. We propose a statistical analysis over 1,000 3D DFNs that are simulated considering the following parameters:

1. The 30 × 30 × 30 cells Cartesian grid sets the scale of the study and provides a uniform prior fracture density map.

2. Fracture strike and dip follow a uniform distribution.

3. Power law fracture length distribution is defined as follows:

n(l) = l -a ; l ∈ [1, 30].
Darcel [2003] shows that the connectivity at percolation is mainly driven by the largest fractures and not by spatial correlations for fracture networks with a small power law exponent (a ≈ 1). A larger power law exponent increases the proportion of small fractures in the network; then, spatial correlations may have a stronger impact on the percolation study. Therefore, we chose to study the impact of sequential seeding process on the connectivity of DFN of power law exponent a = 3. The sequential seeding process uses the parameters described in Table 3. We made the proportion of fractures simulated at each sequence (p(s), Section 3.3) varying in the following range :

1. p(s) = 100% (equivalent to the classical Poisson point process.)

2. p(s) = 10% 3. p(s) = 1% 4. p(s) = 0.1%
During the simulation process, the uniform fracture density was progressively increased until the percolation state would be reached. The percolation state is characterized when a cluster of fractures gets through the whole volume of rock. It means that each of the six faces of the volume have to be intersected by the percolating cluster (Figure 12).

The Percolation Parameter

Many studies have been made to identify a parameter that presents an invariant threshold that characterizes the percolation state of fracture networks. [START_REF] Robinson | Connectivity of fracture systems -a percolation theory approach[END_REF][START_REF] Robinson | Numerical calculations of critical densities for lines and planes[END_REF] first linked the percolation parameters to the fracture density in the case of unit length fracture randomly distributed and randomly oriented in 2D space. He found that the average number of intersection per fracture describes the percolation independently of the fracture orientation distribution law. [START_REF] Balberg | Excluded volume and its relation to the onset of percolation[END_REF] investigated the estimation of the excluded volume around every fracture in a network. The excluded volume is the average volume surrounding an object into which the center of another object cannot lie without intersecting it. This work allows to define a very general expression of a factor describing the percolation of 3D fracture networks. Finally, [START_REF] Bour | Connectivity of random fault networks following a power law fault length distribution[END_REF][START_REF] Gringarten | FRACNET : Stochastic simulation of fractures in layered systems[END_REF][START_REF] De Dreuzy | Percolation parameter and percolation-threshold estimates for 3d random ellipses with widelyscattered distributions of eccentricity and size[END_REF] defined the percolation factor p f , for planar fractures as:

p f = ( n i=1 < l 3 i >)/V (7)
where n is the number of fractures in the fracture network, < l 3 i > is the third moment of the length distribution law, and V is the considered volume of rock.

We generated 1, 000 percolating DFNs for each of the four different sequential seeding processes previously described. Figure 13 shows the probability for these DFNs to percolate against the percolation factor. We observe that the percolation threshold decreases when we decrease the proportion of fractures simulated at each seeding sequence (p(s); Figure 13). This corroborates the findings of [START_REF] Darcel | Corrélations dans les réseaux de fractures : caractérisation et conséquences sur les propriétés hydrauliques[END_REF] that spatial correlations impact the connectivity and the percolation threshold (for networks with a high power law exponent). We observe here that the sequential seeding process creates spatial correlations that can be measured by the correlation dimension D c and allows DFNs to reach the percolation state with a smaller factor p f as compared to uncorrelated DFNs. We also show that the sequential process may create DFNs with a hierarchical organization (Section 3.3) that cannot be discriminated by the correlation dimension (D c = 2.8 for p(s) = 0.1, 1 or 10) but still reduces the percolation threshold. This is probably due to the local organization of fractures induced by the process.

Average Number of Intersections per Fracture

This section investigates the local connectivity of DFNs simulated in Section 4.1. The average number of intersections per fracture normalized by the fracture intensity (P32: fracture area per unit volume in m 2 .m -3 ) has been plotted in function of the fracture length (Figure 14). Simulating all fractures in the same sequence produces fracture networks where the number of intersections per fracture increases with their length. For such DFNs, where fractures are randomly and independently simulated, the local connectivity model follows a single power law model (black line, Figure 14). Decreasing the number of fractures simulated at each sequence changes the global organization of fractures and modifies the model describing the local connectivity. Also, the number of intersections between small objects tends to increase whereas it decreases for large ones (Figure 14). It seems to naturally reproduce the two fracturing regimes described by [START_REF] Wu | An experimental study of the relationship between joint spacing and layer thickness[END_REF], [START_REF] Josnin | A three-dimensional model to simulate joint networks in layered rocks[END_REF], [START_REF] Davy | A likely universal model of fracture scaling and its consequence for crustal hydromechanics[END_REF]. The early initiation of fractures appears in a "poorly" fractured network. It leads to little interactions with other objects that create large fractures with few connections (large fractures are simulated first). The progressive apparition of additional fractures intensifies the interactions between objects and increases the number of intersections per object. This explains why small fractures, that appear when the fracture network is "well-developed", have relatively more intersections.

Figure 14 shows how progressive the fracture apparition should be to observe a smooth transition from the "poorly" to the "well-developed" fracturing regime. If the number of fractures simulated per sequence (p(s)) is too high, a discontinuity appears in the local connectivity model (dot lines, Figure 14). In such cases, the local connectivity model of medium fractures follows the same model as the one observed for uncorrelated fracture networks. The best way to capture the behavior of the transition scale and to mimic the actual fracturing process would be to simulate the fractures one by one. However, simulating 0.1% of fractures per sequence is much more efficient and gives a good idea of the continuous local connectivity model.

Discussion and Conclusion

The presented method uses a stochastic, sequential and non stationary nucleation process to simulate hierarchical and spatially correlated DFNs. It can be conditioned by statistics that characterize the natural fracture network, but also uses rules inspired from mechanical concepts to progressively seed fractures in a 3D model. As compared to the "likely Universal Fracture Model" [START_REF] Davy | A likely universal model of fracture scaling and its consequence for crustal hydromechanics[END_REF][START_REF] Davy | A model of fracture nucleation, growth and arrest, and consequences for fracture density and scaling[END_REF]] that randomly nucle- ates and propagates fractures until their length is statistically equal to the distance to the nearest fracture, we individualize a shadow zone and a stress concentration zone in order to drive the progressive fracture nucleation process. Large fractures, simulated first, organize the nucleation of later fractures.

We show that such a chronological fracture seeding produces a hierarchical organization of fractures that leads to correlated networks. Reducing the proportion of fractures simulated at each sequence likely gets closer to the natural fracturing process and intensifies the spatial correlations between fractures. We identify a threshold proportion of fractures simulated at each sequence from which the estimated correlation dimension stops changing (Section 3.3). This value is experimental and probably depends on parameter settings selected in this study. The computed correlation dimension does not completely capture the spatial features of the simulated networks. Further study need to investigate new a measure of spatial correlations that take into account the anisotropy in the positioning of fracture center. Indeed, the correlation dimension (Equation 2 describes the organization of fracture centers and does not fully describe the hierarchical organization of objects of various size. Decreasing the proportion of fractures simulated at each sequence tends to cluster small fractures at the tips of large ones. It visually decreases the number of X-shaped intersections that may be inconsistent according to the relative chronology of fracture growth (Figure 11). Section 4 shows that DFNs characterized by the same correlation dimension may have significant differences in their local connectivity model due to their spatial organization. The best case would be to consider every fracture previously simulated in the current fracture simulation. However, it would lead to impractically feasible computational costs for large 3D simulation domains.

Spacing and clustering of DFN are important aspects that impact the connectivity of the network. [START_REF] Masihi | A correlated fracture network: modeling and percolation properties[END_REF] proposed a calibration process based on the minimization of the elastic energy dissipated by fractures. They suggest that such optimized networks present a positive correlation between the length and the spacing of objects and they show that the global organization of fractures has an impact on the percolation threshold. Our sequential parent-daughter Poisson point process tends to impose a similar organization by defining shadow zones extents proportionally to fracture size. We have shown that our fracture simulation strategy significantly impacts DFN properties at percolation (Section 4). The sequential seeding produces hierarchical DFNs that favor fracture connections and produce clusters. It allows to reach the percolation with a reduced percolation threshold compared to classical Poisson simulations. The main limitation of the method is that we simulate fractures in their final state and we neglect mechanical interactions and fracture nucleation during the growth. This simplification reduces the complex geometry of fractures that result from the growth and coalescence processes into planar objects. These can lead to an underestimation of the network connectivity. Further studies should take into account an efficient model to reproduce the fracture growth and initiation during the growth by considering ramified and curved objects [START_REF] Gringarten | FRACNET : Stochastic simulation of fractures in layered systems[END_REF][START_REF] Srivastava | Geostatistical simulation of fracture networks[END_REF][START_REF] Bonneau | A methodology for pseudo-genetic stochastic modeling of discrete fracture networks[END_REF]. Incorporating curved fractures in the simulation process could have an additional effect on connectivity properties. Generating realistic geometry for non planar fracture intersections in 3D is an issue hardly discussed in the literature. But, because the proposed method impacts the relative position of fractures according to mechanical considerations, it may facilitate the representation of such 3D branching contacts.

The hierarchy and the spatial correlation of fractures in DFNs are set by the mechanical model associated to each fracture. We propose parameters that describe a symmetrical shape that roughly approximates the first order stress anomalies around fractures. Such a parametrization needs to be simple in order to keep the simulation computationally efficient. It would be interesting to investigate other shapes in order to reflect the asymmetry of stress concentration near sliding fracture tips [START_REF] Kim | Fault damage zones[END_REF] or to use explicit mechanical models [START_REF] Paluszny | Numerical fracture growth modeling using smooth surface geometric deformation[END_REF]. In the same spirit, new rules may be implemented in order to impose spatial correlation on neighboring fracture orientations. This could be suitable to produce "wing-cracks" structures often observed in nature. More generally, a comparison between simulated DFNs and independent data (field, analogs) would be a valuable work to validate the simulation method. In this study, we modified the parameters that set the stress concentration zone (or respectively the shadow zone) extent and intensity to quantify their impact on the fractal dimension of DFN models (Section 3). However, statistics on the fractal dimension that emerges from hundred realizations, show a relatively large standard deviation. As a consequence, it seems difficult to constrain precisely the expected correlation dimension according to natural network characterization. The shadow zone intensity factor (h) is the parameter that most impacts the value of the emerging fractal dimension. An experimental model that links the value of the fractal dimension to the value of h has been identified (Section 3.2). This model is an interesting tool to set practically the fractal dimension for output DFNs. Further work has to be done in order to better understand this emerging behavior and relate it to the physical fracturing process.

acknowledgments

The authors are grateful to the academic and industrial sponsors of the RING-Gocad Consortium managed by ASGA (Association Scientifique pour la Géologie et ses Applications) for funding this work, and to Paradigm Geophysical for providing the Gocad software and API. We would also like to thank Philippe Davy and all reviewers for constructive comments which helped improve this paper. This work was performed in the frame of the "Investissements d'avenir" Labex RESSOURCES21 (ANR-10-LABX-21). The data for this paper are available by contacting the corresponding author at francois.bonneau@univ-lorraine.fr.

Figure 1

 1 Figure 1 Simplified 3D stress model around a fracture. Two ellipsoids individualize the stress concentration zone (in blue) and the shadow zone (in red). Note that the intersection between the blue and the red ellipsoids is a part of the shadow zone. The extent of each zone is deduced from fracture dimensions by applying scaling factors (α s , α p , β s and β p ). In this figure, we set α s = 1, β s = 1, α p = 2 and β p = 0.5.

  (a) and Figure 4 (a). • DFNs made of two sets of vertical fractures (strike of the first set uniformly distributed between 25 • and 35 • and strike of the second set uniformly distributed between 145 • and 155 • ) Figure 3 (b) and Figure 4 (b).

Figure 2

 2 Figure 2 One seeding sequence follows several steps. (a) Fractures simulated previously and their associated influence zones. They will drive the activation of new fractures during further steps. The classical Poisson point process generates possible fracture center locations (b, following a uniform density in this example). Activated fractures are then selected following Equation (1) (c).

Figure 3

 3 Figure 3 3D and 2D views of DFNs obtained with the Poisson point process for (a) a single fracture set DFN; (b) a two fracture set DFN and; (c) a DFN where fractures are randomly oriented. 2D views show the traces of fractures on the red slice. Fractures are uniformly positioned in space. Graphs show for each kind of DFN the correlation function estimation C 2 (r) in black and the estimation of the DFN correlation dimension D c in grey. As the fractal dimension is approximately 3, the experiments show that the classical Poisson point process generates uncorrelated DFNs.

Figure 4

 4 Figure 4 3D and 2D views of DFNs obtained with the Sequential parent-daughter Poisson point process for (a) a single fracture set DFN; (b) a two fracture set DFN and; (c) a DFN where fractures are randomly oriented. 2D views show the traces of fractures on the red slice. Fractures seems organized. The seeding process entails a spacing between fractures and attracts small fractures at larger fracture tips. Graphs show for each kind of DFN the correlation function estimation C 2 (r) in black and the estimation of the DFN correlation dimension D c in grey. The presented sequential activation process builds correlated DFN characterized by a correlation dimension inferior to 3.

Figure 5

 5 Figure 5 E-types computed from 1,000 realizations. (a) shows homogeneous reference maps used. (b) shows E-Type obtained from classical Poisson point process that produces uncorrelated DFNs (D c is the average fractal dimension of DFNs generated). (c) shows E-Type map obtained from sequential Poisson point process that produces spatial correlations. It shows a moderate bias that is interpreted to be related to a border effect.

Figure 6

 6 Figure 6 E-types computed from 1,000 realizations. (a) shows heterogeneous reference maps used. (b) shows E-Types obtained from classical Poisson point process that produces uncorrelated DFNs (D c is the average fractal dimension of DFNs generated). (c) shows E-Type map obtained from sequential Poisson point process that produces spatial correlations.

4.

  Fracture dip and strike are randomly selected in [0 • , 90 • ] and [0 • , 360 • ], respectively.

Figure 7 Figure 8 Figure 9

 789 Figure 7 Sensitivity study on the fracture geomechanical model geometry factors. Several models (a), (b), (c), (d), (e), (f), (g), (h) and (i) have been tested and associated to their average correlation dimension. Increasing the size of shadow or stress concentration zone dilutes the impact of the sequential process and progressively produces non fractal geometries.

Figure 10 Figure 11

 1011 Figure 10Fractal dimension vs proportion of DFN simulated per sequence. Correlation dimension of 2.8 is obtained when 10% or less of the fracture network is simulated per sequence. Then, in the context of our example, a more robust result is obtained as soon as 0.4% of the fracture network is simulated per sequence.

Figure 12

 12 Figure 12 3D DFNs realizations at percolation state for extreme investigated cases. (a) Geometry of percolation DFN obtained simulating all fractures in the same sequence. (b) Geometry of percolation DFN obtained simulating 0.1% of the fractures at each sequence. Percolating clusters are shown in red.

Figure 13

 13 Figure 13Impact of sequential seeding on the percolation threshold. It describes the probability for DFNs to percolate as a function of the percolation factor (p f ). This experiment shows that the percolation threshold decreases with the proportion of fractures initiated at each sequence.

Figure 14

 14 Figure 14Impact of the proportion of fractures simulated per sequence (p(s)) on the average number of intersections per fracture. Uncorrelated DFNs follow a simple power law model that only depends on the fracture length (black line model). Spatial correlations between fractures are obtained by reducing the proportion of fractures simulated per sequence and quantified by the correlation dimension (D c ). It impacts the connectivity model by (1) increasing the number of intersections for small fractures and (2) reducing it for large ones. Reducing the proportion of fractures (p(s)) simulated at each sequence allows to better capture the behaviors of the transition scale (even if the correlation dimension does not change).

Table 1

 1 Average correlation dimension (D c ) over a hundred DFN realizations.

	Poisson point Process	Single set	Two sets	Random orientation
	Classical	3.01 ± 0.07 3.01 ± 0.06	3.00 ± 0.01
	Sequential	2.47 ± 0.09 2.55 ± 0.09	2.58 ± 0.09

Table 2

 2 DFN simulations parameters for E-Types generations

	E-Type map Poisson point Process	Input density map	DFNs Fractal Dimensions
	Figure 5b	Classical	homogeneous (Figure

Table 3

 3 Input parameters for sequential DFN simulations used in Section 3.3 and 4.

We simulate 0.4% of the DFN at each sequence using the same parameters as described Section

2.3.

A Sequential seeding process algorithm Parameters : (1) d (x, y,z) , d ims, or ients: The fracture density map, dimension and orientation distribution laws for the considered fracture set (deduced from field observations and analog or regional knowledge).

(2) s: Total number of sequences to perform.

(3) g: Intensity factor for stress concentration zones.

(4) h: Intensity factor for shadow zones.

(5) DF N i : Discrete Fracture Network from the i th simulation stage.

Output

: DF N i+1