Pierre Kedzierski

Guillaume Caumon

Jean-Laurent Mallet

Jean-Jacques Royer

Pauline Durand-Riard

3D MARINE SEDIMENTARY RESERVOIR STOCHASTIC SIMULATION ACCOUNTING FOR HIGH RESOLUTION SEQUENCE STRATIGRAPHY AND SEDIMENTOLOGICAL RULES

Quantifying the proportions and distribution of rock types is a key step in modeling sedimentary formations, since it greatly determines the quality of petrophysical and dynamic flow models. While forward stratigraphic models well reproduce deposition processes, they rely on inversion for data conditioning, which critically raises the empty space problem. Therefore, several stochastic methods have been suggested to simulate proportions and distribution of rock types in reservoirs respecting the observation data while reproducing realistic spatial patterns.

However, these methods are always limited by the stationarity decision. We propose a method to generate 3D facies probability cubes which account for well and seismic data, stratigraphic interpretation, sedimentological rules describing the spatial distribution of rock types, and sequence stratigraphy principles. Different probability cubes are computed by integrating various information controlling the lithofacies occurrence, for instance the stratigraphic control of the shoreline migration, by defining linear or inequality constraints at the spatial estimation stage. The generated P-field cubes can then be combined considering the redundancy of the data they express such as the tau model. This methodology is demonstrated on typical synthetic data sets and on a North Sea reservoir.

INTRODUCTION

Facies information is a critical ingredient of reservoir characterization studies. Indeed, facies are often used to regionalize the geostatistical modeling of petrophysical properties like porosity or permeability. The modeling of facies is generally done at reservoir scale by conditional geostatistical techniques like indicator kriging [START_REF] Goovaerts | Geostatistics for Natural Resources Evaluation[END_REF][START_REF] Deutsch | GSLIB: Geostatistical Software Library and User's Guide[END_REF], sequential indicator simulation [START_REF] Alabert | Stochastic imaging of spatial distributions using hard and soft information[END_REF], truncated Gaussian simulation [START_REF] Chilès | Simulation of a nickel deposit: problems encountered and practical solutions[END_REF][START_REF] Galli | The pros and cons of the truncated gaussian method[END_REF], multi-point simulation [START_REF] Strebelle | Conditional simulation of complex geological structures using multiple-point statistics[END_REF]. Although these approaches account for both well and seismic data and a given spatial variability model, they rely only in part on geological principles. One major challenge is therefore to introduce more geological concepts within stochastic simulations in order to increase the realism of facies models and to reduce their uncertainty [START_REF] Tamhane | The role of geology in stochastic reservoir modelling: The future trends[END_REF].

In this paper, we present an attempt to fill the gap between sedimentology, sequence stratigraphy and facies simulation by providing numerical rules, which are driven by geological processes and can be accounted for as 3D spatial trends in stochastic simulation algorithms. The notion of facies, as considered in this paper, follows a sedimentological definition, i.e., corresponds to lithological, biological and hydrodynamical criteria [START_REF] Gressly | Observations géologiques sur le Jura soleurois[END_REF]. As such, the spatial variations of facies are mainly controlled by the ratio between the accommodation rate (A) and the sedimentation rate (S) through time, then by diagenetic transformations after sediment deposition. Two main principles have been integrated in our workflow: (1) the A/S ratio describes the seaward and landward migrations of sedimentary bodies in association with the regression and transgression cycles; (2) the sediment volume partitioning characterizes the preservation of sediments in respect to progradations and retrogradations.

INTEGRATING GEOLOGICAL CONCEPTS IN FACIES SIMULATIONS

Several integration strategies have been suggested to incorporate geological concepts into stochastic geostatistical models. For instance, multiple-point geostatistics consider training images as vector of geological knowledge [START_REF] Strebelle | Conditional simulation of complex geological structures using multiple-point statistics[END_REF], and truncated Gaussian simulations [START_REF] Chilès | Simulation of a nickel deposit: problems encountered and practical solutions[END_REF][START_REF] Galli | The pros and cons of the truncated gaussian method[END_REF] use adjacency rules between facies. Most methods can also account for pre-defined local probability density functions (pdf) to reflect spatial trends in the realizations. This last option is critical to obtain non stationary models which reproduce analog observations, but relies on appropriate techniques to estimate the local pdfs. For this, it is possible to use a deterministic mapping between the water depth and facies [START_REF] Massonnat | Breaking of a paradigm: Geology can provide 3D complex probability fields for stochastic facies modeling[END_REF], the accommodation concept [START_REF] Massonnat | Breaking of a paradigm: Geology can provide 3D complex probability fields for stochastic facies modeling[END_REF][START_REF] Leflon | Modélisation des hétérogénéités lithologiques à l'échelle du réservoir pétrolier en milieu marin et fluviatile[END_REF], transition probabilities between the different facies [START_REF] Mallet | Stratigraphic inversion integrating wells, seismic and sedimentology data[END_REF][START_REF] Labat | Simulations stochastiques de faciès par la méthode des Membership Functions[END_REF], and some maps of anisotropy depicting the main directions of sediment transport [START_REF] Mallet | Geomodeling[END_REF][START_REF] Labat | Simulations stochastiques de faciès par la méthode des Membership Functions[END_REF].

We suggest to enlarge this set of concepts to those of genetic stratigraphy. These principles, which characterize the sedimentary structures in response to the accommodation and sedimentation rates, are integrated as numerical constraints in the interpolation of local probability functions of facies. This interpolation is achieved by discrete smooth interpolation (DSI, Mallet, 1992), which finds the optimal values of some unknown variable ϕ(•) by minimizing a quadratic form measuring how much the model violates a set of linear constraints. These constraints may be weighted by the deemed reliability of the data or phenomenon they represent. A DSI constraint c can be expressed as a system of linear equations GEOSTATS 2008, Santiago, Chile defined on the nodes α of some mesh Ω as:

∑ α∈Ω ∑ ν A c (α) • ϕ(α) b c , (1)
where: ∈ { , = , >}, and A c (α) and b c are given scalar coefficients corresponding to the problem at hand. Soft constraints are honored in a least squares sense while hard constraints are strictly honored.

In the case of several pdf's coming from different methods or accounting for different geological concepts, the combination of these information is a way to reduce the model uncertainty if the information is consistent [START_REF] Tamhane | The role of geology in stochastic reservoir modelling: The future trends[END_REF]. Formally, the problem of estimating the probability for some event A i corresponding to the occurrence of the facies i among the n possible facies at a given location conditionally to all M data D m can be expressed as a function ψ of conditional probabilities

P(A i |D m ): ∀i ∈ {1, . . . , n}, P(A i |D 1 , • • • , D M ) = ψ {P(A i ), P(D 1 |A i ), • • • , P(D M |A i )} . (2)
Several exact or approximate formulations exist for ψ functions, depending on the way to aggregate the data [START_REF] Clemen | Combining probability distributions from experts in risk analysis[END_REF][START_REF] Journel | Combining knowledge from diverse sources: An alternative to traditional data independence hypotheses[END_REF][START_REF] Krishnan | The tau model for data redundancy[END_REF]. We suggest to use a model that can account for the redundancy between the various data types when combining several probability cubes coming from several concepts, for instance the tau model [START_REF] Krishnan | The tau model for data redundancy[END_REF].

REPRODUCING LANDWARD AND SEAWARD MIGRATIONS

High-resolution sequence stratigraphy states that the architecture of the sedimentary deposits is controlled by three factors: the sea-level, the tectonics and the sedimentary influx. The sea-level and the tectonic vertical movements are commonly gathered in the accommodation term, defined as the space available for potential sedimentation [START_REF] Jervey | Quantitative geological modeling of siliclastic rock sequences and their seismic expression[END_REF][START_REF] Cross | Controls on coal distribution in transgressive-regressive cycles, Upper Cretaceous, Western Interior, USA[END_REF]. According to the A/S ratio concept, the migration of the system tracts (i.e. transgression and regression) is induced by the variations of the ratio between the accommodation rate A and the sedimentation rate S [START_REF] Cross | Applications of high-resolution sequence stratigraphy to reservoir analysis[END_REF][START_REF] Sonnenfeld | Volumetric partitioning and facies differentiation within the Permian upper San Andres Formation of Last Chance Canyon, Guadalupe Mountains, New Mexico[END_REF]. This genetic concept illustrates the dependence of stratigraphic architecture on rates of change in accommodation and sediment supply (Figure 1):

-A/S > 1: the available space is bigger than the sediment supply, due to a transgression inducing a landward stepping (retrogradation of the sedimentary bodies);

-A/S = 1: the equilibrium between available space and sedimentary influx results in a vertical stacking (aggradation);

-A/S < 1: the accommodation space is smaller than the sediment supply inducing a progradation of the system tracts.

-A/S < 0: available space is reduced (negative accommodation) creating erosion. Eroded sediments are shifted seaward inducing a forced regression. This geological concept can be integrated in the facies probability interpolation when the different facies express topographical attributes.

VIII International Geostatistics Congress

For instance, the Neptune model [START_REF] Massonnat | Breaking of a paradigm: Geology can provide 3D complex probability fields for stochastic facies modeling[END_REF][START_REF] Leflon | Modélisation des hétérogénéités lithologiques à l'échelle du réservoir pétrolier en milieu marin et fluviatile[END_REF] or the Thalassa model [START_REF] Kedzierski | Intégration de connaissances sédimentologiques et stratigraphiques dans la caractérisation 3D des faciès sédimentaires marins[END_REF] associate facies probabilities to the palaeobathymetry or to functions of the water depth. In both cases, facies probabilities are deduced from the function B describing the topography. Moreover, the sequence stratigraphy links the bathymetric variations to the accommodation and sedimentation rates thanks to the accommodation equation [START_REF] Massonnat | Breaking of a paradigm: Geology can provide 3D complex probability fields for stochastic facies modeling[END_REF]:

∂ B ∂t = A -S. (3)
The Thalassa model computes bathymetry curves on wells, then deduces the accommodation rate from the computed water depth and some input sedimentation rate coming from chronostratigraphic interpretation in 3D [START_REF] Mallet | Space-Time Mathematical Framework for Sedimentary Geology[END_REF] (Figure 2). The accommodation rate is then interpolated with DSI in chronostratigraphic space. By combining this computed accommodation rate A and the sedimentation rate S deduced from 3D chronostratigraphic interpretation, we obtain the topographic variation rate ∂ B/∂t as described by (3). This topography derivative is then integrated using paleotopography at the wells using DSI, before being translated into a facies probability cube, still using DSI (Figure 2).

The A/S ratio concept can be integrated at two levels in this methodology: (1) when the bathymetry B is interpolated on wells and, (2) when the accommodation rate is interpolated in the chronostratigraphic space.

Interpolation of bathymetry on wells

The bathymetry is interpolated on wells while: (1) ensuring that the interpolated water depth belongs to the range associated with each facies, and (2) ensuring that available space is created during transgressions and reduced during regressions, as expected from the A/S ratio concept. Indeed, by combining the expression (3) to the A/S ratio and by assuming that S > 0, the following relation comes:

A S 1 ⇐⇒ ∂ B(α,t) ∂t 0, (4)
with ∈ {>, <} for retrogradation or progradation, respectively. This equation is non-linear and can be linearized using a forward finite difference scheme:

∂ B(α,t) ∂t = B(α + 1,t + ∆t) -B(α,t) ∆t , (5)
where the nodes α and α + 1 are two neighboring well samples belonging to the same sedimentary sequence and fault block, and the times t and t + ∆t are the corresponding stratigraphic ages. After linearization and some simplifications, (4) is then translated as a signed gradient DSI constraint of the form:

B(α + 1) -B(α) 0 ⇔ ∑ β ∈Ω A c sg (β ) • B(β ) b c sg , with: (6) A c sg (β ) =    1 if β = α + 1 -1 if β = α 0 otherwise , = > if transgression < if regression , b c sg = 0.

Interpolation of accommodation rate in chronostratigraphic space

From the bathymetry computed on wells, the accommodation rate is estimated by applying equation (3). These local values of the accommodation rate A are then interpolated in chronostratigraphic space. Practically, two criteria are used to VIII International Geostatistics Congress interpolate A. First, the spatial continuity is such that accommodation variations are smooth laterally and driven by the stratigraphic signal vertically [START_REF] Kedzierski | Intégration de connaissances sédimentologiques et stratigraphiques dans la caractérisation 3D des faciès sédimentaires marins[END_REF]. Second, in transgressive (resp. regressive) system tracts, the accommodation rate should be higher (resp. lower) than the sedimentation rate S derived from 3D chronostratigraphic interpretation on seismic data. Then, for any grid cell β where the sedimentation rate is positive, this condition can be formulated as a variable inequality DSI constraint:

A(α) S(α) ⇔ ∑ β ∈Ω A c vi (β ) • A(β ) b c vi , with: (7) A c vi (β ) = 1 if β = α 0 otherwise , = > if transgression < if regression , b c vi = S(α).

Application

These constraints on the bathymetry and on the accommodation rate have been tested to interpolate facies probabilities on a North Sea field. The estimated facies probabilities are consistent with the well data, the 3D stratigraphic interpretation and the expected lateral and vertical facies variations (Figure 3).

PARTITIONING THE FACIES ALONG THE DISTALITY DIRECTION

Facies repartition may be controlled by the physiography of the topographic surface at the deposition time [START_REF] Kedzierski | Stochastic simulation of marine reservoirs through paleotopography modeling[END_REF]. However, the bathymetry and the slope angle are insufficient to completely describe sedimentary systems, and other parameters such as the distance to the sediment source, or the degree of confinement should be considered. For this, the facies proportion estimation should somehow express energy gradients as functions of the sedimentary influx. We are considering in this paper that the facies are distributed along distality axes; the information about distality is carried by orientation and direction vectors, deduced from the stratigraphic record or provided by a trend map.

Stratigraphic record of the distality

The sediment volume partitioning rule states that the sediments are better preserved in distal zones during the regressions and in the proximal ones during the transgressions (Figure 4) [START_REF] Cross | Applications of high-resolution sequence stratigraphy to reservoir analysis[END_REF][START_REF] Sonnenfeld | Volumetric partitioning and facies differentiation within the Permian upper San Andres Formation of Last Chance Canyon, Guadalupe Mountains, New Mexico[END_REF]: during the relative sea level rise (transgression), sediment storage increases landward, reducing the amount of sediments in the distal pole; in a regressive phase (sea level fall), sediments are further transported towards the basin center. This principle means that the gradient of the sedimentation rate, computed on an isochronous surface, points toward the basin center (resp. sediment source) in regressive (resp. transgressive) system tracts. For some distality vector d, we can then write:

∇ uv S = λ • d, (8)
where λ is a positive (resp. negative) scalar in transgressive (resp. regressive) phases. The variables (u, v) represent the palaeogeographic coordinates on the isochron surface, and ∇ uv S the orthogonal projection of the sedimentation rate onto that surface. This equation tells that the distality vector field can be defined from the lateral gradient of S, since the sedimentation rate S is exhaustively known per the chronostratigraphic interpretation. It also provides a valuable self-consistency check for the stratigraphic and sedimentological interpretations.

Constraining the facies probabilities to honor distality vectors

In our work, the local facies proportions are interpolated by membership functions [START_REF] Mallet | Geomodeling[END_REF][START_REF] Labat | Simulations stochastiques de faciès par la méthode des Membership Functions[END_REF], which are similar to the indicator formalism, except that indicator variables need not be ordered. Accounting for distality directions when interpolating these functions can only be performed in a fuzzy manner: whereas both the direction and orientation of variation of the extreme facies can be characterized, only the direction of variation of the intermediate facies can be defined (Figure 5). A first constraint therefore ensures that the probability gradient ∇ϕ is colinear to the distality d (the equation of this constraint, omitted here for brevity, comes from the linearization of ∇ϕ × d = 0). The second constraint forces the interpolated probability gradient ∇ϕ to be oriented in the same (or opposite) direction of a specified vector field g (in which case ∇ϕ • g > 0).

Application

Consider the synthetic data case presented in Figure 6. The distality vector field may be computed from the sedimentation rate using (8). From the observations of five facies on wells, a membership function characterizing the facies probabilities is interpolated with DSI using collinearity between the probability gradients and the trend vectors. For the facies A (resp E), the probability gradient is also constrained to point towards the same (resp. opposite) direction as the trend vectors. The obtained facies probabilities are consistent with the specified trend vectors and the facies realization with a probability field simulation show lateral variations consistent with sedimentological and stratigraphic interpretations.

CONCLUSIONS

We have proposed a method to incorporate sedimentological concepts into the computation of facies probability cubes, based on a mathematical formulation of these concepts. Several probability cubes are generated, and may be recombined as one final facies cube using for instance the Tau model [START_REF] Krishnan | The tau model for data redundancy[END_REF]. Tools assessing the probability of early meteoric diagenesis of carbonates and diagenesis associated with fault systems have also been developed and integrated in the Thalassa code. These are not developed here for space reason, but they are similar in spirit to the method of [START_REF] Henrion | Stochastic approach for modeling karstic systems controlled by fracture network[END_REF]. Results show very good accordance with the sedimentological interpretations, which is extremely valuable. However, the model is sensitive to the input data and the quality of the chonostratigraphic interpretation. Attempts to apply this model without a good seismic quality have failed, because of a lack of control over the sedimentation rate. Another question to ask about this methodology is whether it overfits interpretations, which could be a problem for realistic uncertainty assessment. However, the current practice of modeling, when discounting these concepts, tends to generate models that field sedimentologists are reluctant to accept. Therefore, our method can also be seen as a way to confront information from 3D modeling, seismic stratigraphy and sedimentology: artifacts become visible when the information is inconsistent, while facies probabilities match very faithfully interpretations when information is consistent. Nevertheless, an important avenue for further research lies in the propagation of uncertainties, both inside in our model and in its input parameters (well correlations, stratigraphic uncertainties, sedimentological parameters...)

Figure 1 :

 1 Figure 1: Control by the accommodation rate and the sedimentary flux: the stratigraphic architecture is driven by the A/S ratio changes. Modified from Homewood et al. (1999)

Figure 2 :

 2 Figure 2: Modeling the facies probabilities as the expression of topographic variations constrained by sequence stratigraphic concepts. The A/S ratio interferes in the interpolation of the bathymetry on wells and in the interpolation of the accommodation rate in the chronostratigraphic space.

Figure 3 :

 3 Figure 3: Cross-section of a North Sea stratigraphic reservoir model showing the most likely facies after estimation accounting for the A/S ratio. Note the reproduction of transgressive and regressive cycles. Data courtesy of Total.

Figure 4 :

 4 Figure 4: Stratigraphic expression of the sediment volume partitioning.

Figure 5 :

 5 Figure5: Trend of a probability vector: when the facies probabilities are distributed along trend vectors, the probability gradient of the most proximal and most distal facies has an unique direction while the probability gradients of the other facies is increasing and then decreasing.

  Figure 6: Estimating the facies probabilities according to trend vectors coming from the stratigraphic interpretation of the GeoChron-based sedimentation rate, or coming from an interpreted facies map. Red colors correspond to low probabilities while white ones represent high probabilities.

GEOSTATS 2008, Santiago, Chile

ACKNOWLEDGEMENTS

We thank the members of the Gocad Research Consortium for funding this research, especially PDGM for providing the Gocad software and API. Total E&P UK provided the North Sea data set, and Stewart Clark is greatfuly acknowledged for stimulating discussions.


Vincent Henrion

Jeanne Pellerin

Guillaume Caumon

A STOCHASTIC METHODOLOGY FOR 3D CAVE SYSTEM MODELING

This paper describes a methodology to generate stochastic models of 3D cave systems. It combines existing geostatistical techniques, namely object-based and variogram-based geostatistics. The first is used to simulate Discrete Fracture Networks (DFNs) which in addition of bedding planes drive the development of karst features. The discontinuity network and the matrix are jointly discretized into a graph of connectivities. Preferential flow paths determining the topology of the cave system are then extracted using graph search algorithms. The distance map to the selected paths is used to conditioned the spatial extent of karst features, i.e., the probability to simulate karst decreases moving away from "dissolution paths".

Eventually, multiple realizations of a distance cutoff are generated with Sequential Gaussian Simulation to perturb the distance function. Resulting models efficiently reproduce the spatial variability of karst features in term of geometry and their spatial organization in term of connectivity.

INTRODUCTION

Carbonate formations exposed to dissolution processes are characterized by extreme spatial heterogeneity due to the presence of complex networks of highly permeable features (e.g. karst conduits and cavities) embedded in less permeable fractured rocks. Cave and paleocave system considerably impact the fluid circulation respectively in groundwater and petroleum reservoirs, involving many engineering issues. For water resources management, karstified aquifers are particularly vulnerable to contaminants which can be transported over wide distances with no or little dilution [START_REF] Bakalowicz | Karst groundwater: a challenge for new resources[END_REF]. In oil and gas reservoirs, paleokarsts form extended brecciated and fractured zones which act either as seals or conduits depending on the type of cementation [START_REF] Loucks | Paleocave Carbonate Resvoirs: Origins, Burial-Depth Modifications, Spatial complexity, and Reservoir Implications[END_REF]. Static and dynamic modeling of such reservoirs calls for properly capturing the spatial organization of karst network and the geometry of individual conduits, i.e. generating realistic 3D karst systems. The construction of such inter-connected tortuous conduits conditionally to possibly dense observation data is a difficult task. Pixel-based multi-point simulation could be used for this purpose, provided an adequate training image, but is known to break the continuity of large objects, which is a problem for karst networks where connectivity is paramount. A stochastic methodology has been developed for creating complex and realistic 3D cave system models that honor field data. It can also accommodate many different observations such as sources and sinks location, water level, bed thicknesses and other subsurface observations. These elements provide constraints for karst conduits location and their geometry. This approach generate multiple equally probable models of 3D cave systems which topology and geometry is consistent with structural geology and other constraints.

APPROACH OVERVIEW

Karst originates from the dissolution of rock by groundwaters. Water preferentially circulates in largest apertures usually fractures and bedding planes. Dissolution causes their enlargement and conduits progressively develop around these preferential flow paths and get organized in a connected karst network. Although the proposed method does not quantitatively model the dissolution process governing karst development, it attempts to create realistic 3D cave systems by stochastic simulations of the main factors driving the topology and geometry of karst conduits. Since fractures and bedding planes act as drivers for karst development, the first step of the method consists in generating 3D models of rock discontinuity surfaces consistent with structural geology of the considered region. Karst processes are however highly selective, enlarging only some of the original discontinuities and developing a conduit network progressively organized in the same way as fluvial systems. For this reason, the fracture network and the matrix are discretized into a graph of connectivities allowing extraction of preferential flow paths using graph search algorithms. This screening test attempts to keep only the most favorable paths for karst development depending on their connectivity and transmissibility; some data conditioning may also come into play at this stage, for instance by selecting the cave outlet or sinkhole as a root for the graph search algorithms. Once the "dissolution paths" have been extracted, the main problem is to simulate the spatial extent of open cave passages. For this purpose, the Euclidean distance map related to the selected flow paths is computed, and used as a constraint for karst development. Then, multiple realizations of a correlated random noise, termed distance cutoff, are generated using sequential Gaussian simulation. The composition, given a criterion, of the distance map and of the distance cutoff produce final images of karst systems (Figure 1).

The whole procedure involves the following steps (Figure 2 The implementation of this methodology is detailed in the following sections.

RESULTS AND DISCUSSION

Structural modeling

Karst formation results from near-surface processes where the two major controls are chemical kinetics and discharges. Chemical kinetics are set by a combination of climate, mineralogy, CO 2 concentration, water table and interaction with other dissolved components. Although these factors can vary considerably, both spatially and temporally, it is unlikely that their variations have a first order impact on flow paths in a given carbonate terrain. In contrast, there are great variations in discharge from one flow path to another and this is the main control over which early flow VIII International Geostatistics Congress [START_REF] Palmer | Origin and morphology of limestone caves[END_REF] estimates that 99% of mapped continental cave passages are initiated along fractures and bedding-planes. However, only a few fractures can evolve into open conduits because karst processes are highly selective, enlarging only the most favorable discontinuities. At an early stage, fissures are narrow and the flow is dispersed among many different routes, each of them with its own overall hydraulic gradient, mean fissure width, flow length and mean discharge. The influence of these parameters have been investigated by geochemical and digital models (see e.g. [START_REF] Dreybrodt | Basic processes and mechanisms governing the evolution of karst. In Speleogenesis and Evolution of Karst Aquifers[END_REF]). Results show that the ideal situation for cave passage development is large fissures with short flow paths, high discharge and steep hydraulic gradient, embedded in thick layers. These facts are important because they can aid in the prediction of cave passages distribution and geometry.

For these reasons a structural model of the area of interest is build integrating all available subsurface informations. Discrete Fracture Networks are simulated using object-based simulation (Stoyan et al., 1995) and constitute with bedding planes the discontinuity network around which karst will preferentially develop. A key issue from a geostatistical perspective, is the inference of the parameters of the stochastic models (fracture density, orientation and size). All sources of informations (wells, seismic, analog outcrops) provide a first evaluation about such fracture density and fracture properties. None of these attributes are however well constrained in the inter-wells space due to the poor seismic resolution and the clustering of 1D data along wells. To more realistically model the spatial distribution and the properties of fractures, a variety of numerical techniques have been developed. Basically, we can distinguish geometrical approaches and geomechanical one. The first consist to look for geometrical attributes revealing zone of intense deformation (structural curvature, distance to major fault) and to relate these to fracture density [START_REF] Lisle | Detection of zones of abnormal strains in structures using gaussian curvature analysis[END_REF][START_REF] Fischer | Predicting the orientation of joints from fold shape : Results of pseudo-three-dimensional modeling and curvature analysis[END_REF][START_REF] Ouenes | Practical application of fuzzy logic and neural networks to fractured reservoir characterization[END_REF]. Geomechanical approaches consist to compute stress and strain field which is then combined with a failure criterion to infer fracture orientation and likely locations [START_REF] Bourne | Predictive modelling of naturally fractured reservoirs using geomechanics and flow simulation[END_REF][START_REF] Macé | Integration of fracture data into 3D geomechanical modeling to enhance fractured reservoirs characterization[END_REF]Maerten et al., 2006;Maerten and Maerten, 2006;[START_REF] Moretti | Working in complex areas: New restoration workflow based on quality control, 2d and 3d restorations[END_REF]. Geomechanical approaches compared to geometrical indicators have the quality to be based on physical laws that govern geological deformation, which make them more reliable to predict fracture attributes. The resulting fracture orientation and density models are then used to inform further stochastic simulations of DFN.

Extraction of preferential flow paths

The 3D structural model is discretized into a graph of connectivities, namely a PipeNetwork, using the method developed by Vitel and Mallet (2005). It basically consist in a set of segments linking a set of nodes which structure derives directly from the representation of a grid based on connectivities. The PipeNetwork discretization provide a flexible framework to handle complicated geological configurations as in fractured reservoirs. The matrix and the DFN are jointly discretized. The resulting PipeNetwork is a completely unstructured graph constituted of three SubPipeNetworks : one for the matrix, another for the fractures VIII International Geostatistics Congress and the last one for the connections between matrix and fractures. Fracture nodes and matrix nodes hold their own set of properties (porosity, permeability, fracture apertures) and pipes hold hydraulic properties (conductivity, transmissibility). A screening test with graph search algorithm is then performed on the PipeNetwork in order to extract preferential flow paths where dissolution is more likely to occur. Path research is made between input and output points representing sinks and sources of the karstic system. Graph search algorithms finds the least-cost path between given nodes using a distance-plus-cost function heuristic to determine the order in which nodes of the graph are visited. In the example presented here (Figure 2(e)), the inverse of the hydraulic property stores on the pipes of the graph have been used.

Field data conditioning

Another issue of cave system modeling is the conditioning to field observations. Data can come from many different sources (wells, exploration survey, geophysics) and provide a set of points where karst or matrix appears.

The conditioning to observations revealing karst or matrix occurrence is treated during the search path step. The data are ported on nodes of the PipeNetwork in order to force preferential flow paths to pass close of karst observations and away from host rock observations. To ensure karst conditioning, the distance between a flow path and the data point must not exceed the minimal size of conduits to be simulated. Hence, a node of the PipeNetwork is randomly chosen in the minimal distance neighborhood of the data point and is included as a conditioning node in the path search (Figure 3). On the contrary, to respect matrix data, the distance between a path and data point must exceed the maximal size of karst conduits. For this, the transmissibility of the segments linking nodes in the maximal distance neighborhood of the data are reduced to prevent any path from crossing these nodes (Figure 3). In addition to hard data conditioning, the structure of the PipeNetwork allows to easily integrate secondary information by transposing it in a coherent manner on the nodes and pipes properties. For instance the distribution of karst, and particularly the maximum depth at which karst can develop, is strongly linked to the water level of the phreatic table. Base level of karst system realizations can be constrained by modifying transmissibilities on the pipes with a multiplier coefficient, greater ones along the water level and lower ones below this surface.

Geostatistical simulation of karst conduits

The previously extracted preferential flow paths define the spatial organization of the cave system. The last step of the methodology consist in simulating the geometry of conduits. The distance to the flow network can be seen as a map of the potential for karst development, i.e. the probability to have karst facies decreases with distance to flow paths and becomes null beyond a given threshold. on the probability density function (PDF) and the variogram model of the cutoff used to inform the stochastic simulation. The first control the final dimensions of the conduits (Figure 4) while the parameters of the variogram model control the undulations of the karst surface (Figure 5). Resulting 3D karst models present realistic spatial organization and geometry of individual conduits. Moreover they can be easily conditioned by secondary informations. For instance conduits size usually increase toward the source of the karstic system. This could be rendered by using kriging with a trend model. Another possibility is to use co-kriging to localy correlate the conduit extent to bed thickness. VIII International Geostatistics Congress Bulleetin, vol. 103, pp. 1-21. Stoyan, D, Kendall, W and Mecke, J (1995). Stochastic Geometry and its applications. John Wiley & sons, New York. Vitel, S and Mallet, JL (2005). A generic network for flow grids. In Proc. 25th Gocad Meeting, Nancy, France.

VIII International Geostatistics Congress

  ): -3D modeling of the structural geology including stochastic realization of a GEOSTATS 2008, Santiago, Chile

Figure 1 :

 1 Figure 1: The distance map to "dissolution paths" (a) is computed (b) and compared to a random distance cutoff (c), to produce realistic cave geometry (d)

Figure 2 :

 2 Figure 2: (a) Structural model and (b) Discrete Fracture network. (c) Discretization into a graph of connectivities (PipeNetwork) and (d) definition of input and output point. (e) extraction of preferential flow paths. (f) Distance function to the selected paths and (g) random cutoff distance simulation. (f) Realization of final cave system .

Figure 3 :

 3 Figure 3: The neighboring nodes of each data points are searched in the graph. Left: Conditioning of host rock data point: the transmissibility value is reduced for all the segments. Right: Conditioning of karst data: one point is randomly selected and its karst value is set to 1.

Figure 4 :

 4 Figure 4: Karst realizations with, from left to right, increasing maximum of a uniform PDF.

Figure 5 :

 5 Figure 5: Karst realizations with, from left to right, increasing ranges of a Gaussian variogram model.

ACKNOWLEDGMENTS

This research was performed in the frame of the GOCAD research project. The companies and universities members of the GOCAD consortium are acknowledged for their support.

CONCLUSION

The methodology proposed in this paper enables to render the complexity of cave systems. Their organization in space is addressed through 3D geological modeling and object-based simulation of discrete fracture networks. Special care must be brought to the simulation of multiple fracture networks because it provides the first order variability of the proposed simulation workflow. Indeed preferential flow paths along fractures and bedding planes are retrieved from an appropriate discretization of the structural model and defines the final topology of the cave systems. The simulation of cave passages around the fracture gives the second order variability. The simple idea of the method, based on geostatistical perturbation of the distance to the preferential flow paths, proves its potential to produce different and complex geometries of cave passages.


Alain Jourdain

come

Introduction Data & Methods

Results

Conclusions

Alain JOURDAIN Department SHSC, Arènes Research unit cnrs 6051 alain.jourdain@ehesp.fr Résumé Cette recherche est une contribution à la connaissance des effets de la pollution atmosphérique sur la population. Une mesure de la concentration en dioxyde d'azote et un suivi de la mortalité de la population ont été réalisés sur les mêmes communautés entre1999 et 2011 en France. Il y a un gradient de mortalité croissant avec le niveau de pollution. Le risque relatif de mortalité entre une pollution moyenne et une faible pollution est estimé à 5%. Ces conclusions sont limitées par la taille de l'échantillon de l'enquête.

XXVIIIe Congrès international de la population de l'UIESP Cape Town, South Africa 29 October to 4 November 2017.

This study is based on data from the HID survey conducted in France in 1999 with a national representative sample of households whose vital status was followed up to 2012. (n = 16896).

A sub-sample of the HID survey was created as part of the "DISPARITES" research program. It identifies the living area of people (iris 2000) thus allowing to determine a social deprivation score and a level of pollution. (n = 1283) http://www.equitarea.org/ The vital status of those generations was followed by Civil registration between 1999 and 2011.

Demographic analysis of the effects of pollution on mortality in four French metropoles.

Nitrogen dioxide (NO2) is one of the most hazardous atmospheric pollutants. (Extrapol, 2006). It has adverse health effects because long-term exposure can impair lung function and increase the risk of respiratory problems (Crouse et al., 2015) .

The arguments for an environmental health policy are built on the burden of mortality or morbidity. (European Environment Agency, 2016).

In this study we propose to measure :

-the life expectancies of generations exposed to different levels of NO2 pollution, -the number of deaths of the table -the relative risk of mortality attributable to different levels of pollution. The specific effect of pollution is 5% when one passes from one class of the terciles to another.

There is a gradient of mortality , with a significant difference of 4 years of life expectancy between the exteme groups of pollution.

A generation of French population subject to heavy pollution would know during its history 42 350 deaths annually.

By controlling gender, socioprofessional status, social deprivation and disability, the relative risk of mortality between average and low pollution is estimated at 5%.

Difference in life expectancy between pollution groups is a good indicator of public health but it omits other variables that could influence this difference.

Mortality rates and annual deaths are good variables to communicate on the pollution burden, but there is a possible bias on age.

The Cox regression, which measures a relative risk of mortality, is less relevant in terms of communication with the general public, but it provides better control over interpretation biases.

Great interest of the demographic approach, in a field where data are scarce; but the main limits of this study are:

-the insufficient sample size (n = 1283), -characteristics of the population are measured at the beginning of observation, deaths are observed during 11 years and the level of pollution is measured at the end of the period. -the assimilation of pollution to NO2 to all atmospheric pollution is a final limitation to the results.

Table 1 :

 1 Number and rate of death by level of NO2 pollution in the stationary populations.

Table 2 :

 2 Cox regression applied to the distribution in terciles of pollution

https://www.researchgate.net/publication/320410481 Full paper available in French and Eglish in the conference website and in:


Ghislaine Gueudet
email: ghislaine.gueudet@espe-bretagne.fr

Birgit Pepin²

in mathematics and resources mediation at entry to university

Keywords: epistemic mediation, pragmatic mediation, resources, secondary-tertiary transition

In this paper we study the use of resources by students in their mathematical work at the beginning of university. The institution offers a variety of resources: lecture notes, books, exercises, websites, to name but a few. Leaning on a theoretical framework by Rabardel, we argue that the university teachers expected an epistemic mediation of these resources, as they supported student learning of (higher level) mathematics. However, analysing two case studies (one in the UK and one in France) we observe that the actual use of resources by novice mathematics students corresponded to a pragmatic mediation, as they searched for worked examples and very similar to their use of resources at secondary school.

RESOURCES AT THE SECONDARY-TERTIARY TRANSITION

The secondary-tertiary transition is a moment of important institutional evolutions: e.g. in terms of the mathematics taught; the mathematical practices expected from students; the institutional support for student learning [START_REF] Gueudet | Investigating the secondary-tertiary transition[END_REF][START_REF] Pepin | Using the construct of Didactic Contract to understand student transition into university mathematics education[END_REF]. In this article w with mathematics, and the expected use of these resources (expected by the institution) significantly change from secondary school to university mathematics education. When entering university, it is assumed that students develop new ways to use resources, mathematical texts in particular. Previous research (e.g. [START_REF] Rezat | The Utilization of Mathematics Textbooks as Instruments for Learning[END_REF]) contends that at secondary school even grade 12 students used their mathematics textbook mostly to search for worked examples, in order to learn rules and how to apply these rules to tasks similar to the tasks they worked on with their teachers in class. The same holds true at the beginning of university: [START_REF] Lithner | n university textbook exercises[END_REF] shows that was mostly oriented towards solving exercises; and that students searched for surface similarities between exercises, in order to choose a procedure. Our hypothesis is that this use of textbooks and other mathematical curriculum material by university students derives from similar practices at secondary school, in particular their use of textbooks. Hence, our research question is the following:

-which use of resources is expected by university mathematics teachers, and how does this compare with the actual use of resources by students?

In the next section we explain the theoretical frame we used ) instrumental approach. Subsequently, we present a case study from the UK TransMaths project 1 at a general level; and a case study in France in the area of Number Theory.

RESOURCE MEDIATION: THEORETICAL FRAME AND METHODS

In terms of theoretical frame we refer to the instrumental approach by Rabardel (2002Rabardel ( /1995)). He distinguishes between an artefact -produced by humans, for an aim of human activity; and an instrument -developed by a subject along his/her goaldirected activity with this artefact. The instrument incorporates the artefact (or parts of it) and a scheme of use for this artefact [START_REF] Vergnaud | Toward a cognitive theory of practice[END_REF]. Following [START_REF] Vygotsky | Mind in society[END_REF], [START_REF] Rabardel | From computer to instrument system: A developmental perspective[END_REF] consider that the goal-directed activity of a subject is mediated by instruments, in particular between the subject and the object of his/her activity. This mediation has a pragmatic value: the instrument permits to reach the aim of the activity; it contributes to the production of an outcome. At the same time it has an epistemic value: it contributes to the development of the subject him/herself and of his/her understanding of the object. This pragmatic/epistemic distinction has been stressed in particular by [START_REF] Artigue | Learning Mathematics in a CAS Environment: The Genesis of a Reflection about Instrumentation and the Dialectics between Technical and Conceptual Work[END_REF] in her work about the use of CAS (Computer Algebra Systems). She identified pragmatic and epistemic values in the instrumented techniques developed by students: these techniques both permitted to reach an aim, for example solve an exercise; and to build new mathematical knowledge, enabling the students to solve further exercises. Drawing on the work of Rabardel, we developed the theoretical frame of the documentational approach (Gueudet, Pepin & Trouche 2012). In this approach, instead of artefacts we considered resources of different kinds, following the definition proposed by [START_REF] Adler | Conceptualising resources as a theme for teacher education[END_REF]: everything that is likely to re-source the activity. Our previous research has mostly been concerned with teacher interaction with resources; our focus here is on We distinguish between available/proposed resources and resources-in-use, comparing the resources available for students (resources offered by the institution), and the resources actually used. We already analysed the two cases presented here in terms of links between the Didactic Contract [START_REF] Brousseau | Theory of didactical situations in mathematics 1970-1990[END_REF]) and the use of resources [START_REF] Gueudet | Didactic contract and secondary-tertiary transition: a focus on resources and their use[END_REF]. In this paper we extend the epistemic/pragmatic distinction to the mediation realised by various kinds of resources (Gueudet, Pepin & Trouche 2012) intervening in students work, and investigate these mediations. The resource mediation can be represented by the following figure (figure 1).

with mathematics.

Students at university use resources of various kinds: paper resources presenting the text of the lecture; their own notes; lists of exercises; textbooks; online resources found on the Internet; etc. These resources mediate the interaction between the student and the object of his/her activity; it has both a pragmatic and an epistemic value. investigated through interviews with teachers and students, and we conducted such interviews in our two cases in the UK and in France.

In the UK, we explored the use of resources at a general level, at a medium-size for all mathematical content areas (taught during the first year), focusing on the general organisa Beside conducting student interviews (with selected students), we surveyed all students of that particular mathematics course; we observed lectures, and interviewed lecturers, in addition to other support staff (e.g. teaching assistants at university) and teachers at secondary school. We also collected documents at both institutions (school, university, in the UK). For the work presented here, we selected a case study subject who is an ethnic minority student, Simar and his (ethnic minority) friends, all studying mathematics in the faculty of mathematics at the same university. We followed Simar over approximately two years: starting when he transited from a local upper secondary school into a university mathematics course, and into his second year at university. Over that period observations and interviews were conducted at several data points: (1) in his previous school (with his mathematics teacher/s); (2) at entry to university (with Simar and his friends; and with lecturers and support tutors); (3) towards the end of the first year at university (with Simar and his friends); and (4) during the first semester of his second year at university (as under ( 2)). For this study, interviews Other students

Teacher

Student

Object of the activity: work in mathematics

Resources: textbook, list of exercises, Internet with Simar (3); focus interviews with his friends (2); and interviews with selected lecturers (5) were analysed using our theoretical frame explained earlier.

The case in France has been conducted with a focus on a particular mathematical content, Number Theory. We collected data in a medium-size university where Number theory was taught in a first year teaching unit. During the academic year 2014-2015, we interviewed the teacher responsible of the course about the use of resources she expected from students. We also proposed an online questionnaire to the 140 students about their use of resources. We collected 85 answers (around 61 %). In 2015-2016 we complemented this study by proposing exercises as homework to a group of students, and interviewing three of them about their use of resources for solving these exercises. The collection of data at different levels, we claim, is likely to give us insights at different phenomena, which (it is hoped) complement each other.

CASE 1: STUDENT RESOURCES AND THEIR MEDIATION AT GENERAL LEVEL

From interviews with Simar (and his friends), who studied at City University, we could identify the main resources used in their first year of study: the lecture and lecture notes; the coursework; notes and tutorial notes made during tutorial and/or with his friends/study group. It was clear that these resources were quite different, in nature and quantity, from what students were used to at school: at school students had one textbook (which was portraying mathematics as something that one can , whereas at university they were expected to work with many resources: e.g. different textbooks; lecture notes; examination papers; etc.

At City University the main resources were clearly the lecture, usually in halls of up to 300 students, and the lecture notes provided by the lecturer/professor (sometimes supported by a textbook). Different lecturers had different styles of lecture notes. Some lecturers would produce hand-written notes projected onto a screen (and talk students through the content during the lecture) (e.g. calculus). Simar and his friends would copy these notes, most of the time with little or no understanding.

However, Simar talked quite enthusiastically about one of his lecturers, and he pointed to what he (and his peers) would regard as a good lecture and lecture notes.

Another lecturer would provide notes that students had to read in advance of the a way, students were not only obliged to prepare the lecture in advance (in order to be able to understand the notes in the lecture and fill in the missing text), but they In addition to lectures, the coursework (provided once a week) was to support student understanding of the lecture, through exercises. Simar and his friends/learning group were clear that unless the coursework was well aligned with the lectures, it did not help their understanding of the subject area (see Calculus as compared to geometry lectures/coursework). Indeed, in some cases students did not know what to ask in tutorial time, or in lectures, so little had they understood of the topic area. Other resources included textbooks (suggested/approved by the lecturers), but these were seen as less helpful than the lecture notes and coursework (provided by lecturers and tutors) However, the same resources (e.g. lecture notes) were often evaluated very differently by students, in terms of support for their learning, so much so that Simar (as student representative) had asked for a change in form and practice concerning .

Interestingly, in terms of lecture notes students distinguished between different types useful for understanding and coursework (and tests); (2) you want to come to At the same time institutional practices, such as lectures, and accompanying resources played a crucial role in the ways that mathematics, and what it meant to was portrayed. On the basis of video footage of selected lectures and pre-and post-video stimulated recall discussions with lecturers one could identify meanings that were attached to particular practices. Particular lectures reflected the kinds of things that developed through Geometry and Linear Algebra; --practical and context relatedness was regarded to be developed through Statistics.

However, it was clear that during lectures student would not learn how to work as a neither did students expect this from lectures and lecture notes. What studen worked examples suitable for studying for examinations.

to do your coursework, and this is how I tend to them more and during the tutorials, and I think the tutorials group interview)

In terms of mediation of resources, in particular lecture notes, it can be argued that for students they had mainly pragmatic value: Simar and his friends were content, if However, for university lecturers the mediation of (for example) lecture notes had , and lecture notes (and lectures) would show them how worked. How that could be learnt was not clear, except for alignment with what the lecture notes showed as examples. fact, at City University one lecturer realized the provide students with the language they needed to appreciate the epistemic side of the subject.

CASE 2: STUDENT RESOURCES AND THEIR MEDIATION IN NUMBER THEORY

Our investigations took place in a first year teaching of Number Theory spanned over twelve weeks, with four hours each week (two hours of lecture, and two hours of tutorial). The first half of the course concerned logic, sets and combinatorics; the second half more directly number theory, with Euclidean division, Euclid algorithm, prime numbers and congruencies.

Use of resources by students and pragmatic mediation

In 2014-2015, at the beginning of this course the students were provided with a exercises. They could also access complementary resources, on the webpage of the teacher responsible for the teaching: previous exam texts; references of books; links towards online exercises. We proposed an online questionnaire to the 140 students concerned and obtained 85 answers. Only 52% of these 85 students declared that they found the polycopie useful. They considered that the text of the lecture of their teacher was enough, and used the polycopie only before the final exam (83%). Moreover 90% would have liked to find worked examples in the polycopie; and 44% looked for additional resources on the Internet, in particular worked examples. We contend that these answers evidence that the polycopie mediation remains pragmatic for the students, similar to their use of textbooks at secondary school where number theory is limited to the application of some techniques [START_REF] Battie | Number theory in the national compulsory examination at the end of the French secondary level between organizing and operative dimensions[END_REF]. Alike students in the UK, they search for worked examples in order to reproduce techniques, whereas the teacher expects that the polycopie has a strong epistemic value, and is used to work on the course: learn definitions, understand proofs of theorems etc. (declared by the teacher responsible for this teaching).

Mathematical resources for Number Theory

In 2015-2016, following the results of the questionnaire evoked above (a report about the answers was presented to the teachers of the number theory unit), no polycopie was given to the students. A book was recommended instead, together with a website, Braise 2 proposing exercises associated with different mathematical texts: description of methods, extracts of the course, hints, partial solution etc.

The teacher proposed to the students the following exercise as homework: This exercise can be solved by searching for a and b such that a(2n+3)+b(n 2 +3n+2)=1. There is a specific difficulty in this exercise, since a had to be itself of the form un+v, where u and v were constant integers.

A student can solve this exercise without any explicit use of resources; nevertheless, since the homework was given at a stage where such tasks are not yet familiar to students, we consider that this will probably not happen. Amongst the resources provided by the institution, the students can use their course notes, in particular to find the definition of coprime. They can also use the notes taken during the tutorial, where an exercise using a similar method has been done: such that m divides both 8n+7 and 6n+5. Show that m=±1 (Exercise 8). In this exercise, the students also need to find a linear combination of 8n+7 and 6n+5, which does not depend on n. The students can also visit the Braise website; they could find on it w also be used here as a resource. It does not mean that the students will only need to reproduce the same method: in particular, the presence of n 2 in exercise 28 requires a significant adaptation of the method. We claim that, in such a case, the resources mediation has an important epistemic value.

16 students did the proposed homework; ten of them proposed a correct solution, and for six of them proposed a wrong solution. We met three of these students for an individual interview about the resources they used to solve the exercise: Brian, and Franck, who did not succeed; and Tom who found a correct solution. Tom used exercise 8 that he found in his tutorial notes, and correctly adapted the method. Franck used exercises that he found on the Internet (but not on BRAISE), which he identified as useful in terms of including the idea of Euclidean division of polynomials, and divides n 2 +3n+2 by 2n+3.This lead him to conclude that the gcd is (-). Brian searched his lecture notes, found p is prime and n an integer, then p and n apply it but realised that 2n+3 is not always prime. Then he searched grade 12 textbooks, found the linear combination method but only with constant coefficients, and thought that the coefficient cannot depend on n. We claim that, while for Tom and Franck epistemic mediations of the resources took place, for Brian the mediation was limited to a pragmatic aspect. He searched for a method that he wanted to apply without any adaptation. Franck took personal initiative searching for resources that where not proposed by the teacher. He tried to build an original method, but was not successful in controlling its correctness. It is noticeable that, while Brian came directly from secondary school to university, Tom and Franck have had previous experiences at university: Tom did a first year of law studies before deciding to study mathematics, while Franck studied two year We can assume that the influence of secondary school was less important for them, as they had previous experiences at university.

CONCLUSIONS

The changes in resources (for teachers and for students), and in the use of these resources, at the secondary-tertiary transition have been under-researched. We contend that the study of resource use is an important theme for research, likely to deepen our understandings of teaching and learning processes (initiated or supported by resources) at the beginning of university. In the two case studies presented here we observe that the institution provides the students with numerous resources, mainly mathematical texts. According to the teachers, the epistemic mediation of these resources should support students transiting from school to university mathematics, both in their ways of learning mathematics (more self-regulated work, more autonomous reading of mathematical texts, see e.g. Farah 2015), and in their mathematical work, so that it would become similar to the work of a . However, asking the students about their actual use of resources led to a different picture. The learning at university seemed to be based on listening to the teacher (in lecture), writing down notes, trying similar worked examples, reading the polycopie for exam preparation (in France) in other words an alignment based on a kind of apprenticeship learning. We contend that, based on the two cases we studied, the pragmatic mediation of resources took over the epistemic aspect, at least in the first year. Students used worked examples and lecture notes in order to produce the desired results. At the level of mathematical content, in our case number theory, we detected a potential epistemic aspect in the use of worked examples and lecture notes, when a significant adaptation of a given method was needed. associated with the change from school to university mathematics education take longer than expected (by university staff), and more awareness and didactical flexibility (from the side of university staff) might help students to bridge this gap more successfully.

Figure 2 .

 2 Figure 2. Exercise given as homework (Let n be a positive integer, prove that 2n+3 and n 2 +3n+2 are coprime).

NOTES 1 .

 1 TransMaths project, University of Manchester: http://www.transmaths.org 2. http://braise.univ-rennes1.fr/braise.cgi

Acknowledgement

The second author of this paper recognizes the contribution made by the TransMaths team in the collection of data, the design of instruments and the project, and discussions involving analyses and interpretations of the results. She would also like to acknowledge the support of the ESRC-TLRP award RES-139-25-0241, and continuing support from ESRC-TransMaths award(s) RES-139-25-0241 and RES-000-22-2890.


Frédéric Livet

Mathieu Fèvre

Guillaume Beutier

Federico Zontone

Yuriy Chushkin

Mark Sutton

Measuring the dynamical critical exponent of an ordering alloy using x-ray photon correlation spectroscopy

The dynamics of the order fluctuations in the AuAgZn 2 close to the critical point (T c ≃ 609 K) was observed by coherent x-ray scattering. With the high beam intensity of the ID10 ESRF beamline and with the new pixel detector, the dynamics was measured with a few tens of millisecond resolution. The intensity connected to the diffuse scattering corresponding to fluctuations was unambiguously distinguished from the surface pretransitional ordering occurring in this system close to T c . The variations of the fluctuation time with temperature and wave vectors were measured in this system belonging to the universality class of Ising second order transition with a non-conserved order parameter. The direct observation of the critical slowing down in the vicinity of the second-order transition led to an estimate of the dynamic exponent z ≃ 1.96(11), in rough agreement with theory (model "A" of Ref. 1).

INTRODUCTION

X-ray photon correlation spectroscopy (XPCS) has raised considerable interest to probe equilibrium and non-equilibrium fluctuations in soft and hard systems [2][3][4] at the nanoscale and on timescales increasingly smaller [5,6]. In alloys, the technique was used e.g., to study domain formation [7][8][9] or phason fluctuations [10]. In AuAgZn 2 , a second order transition occurs at 609 K [11,12]. Above this critical temperature T c , the alloy has a simple cubic B2 structure where Au and Ag atoms share the same simple cubic atomic sites. Below T c , the alloy has a face centered cubic L2 1 structure, as shown in Fig. 1. This transition ("Heusler") is equivalent to a SC antiferromagnetic Ising system [13]. In the vicinity of the critical temperature, the fluctuation length ξ is observed to increase to ξ > ∼ 500 Å, with the scaling law: ξ ≃ 1.76 × ((T -T c )/T c ) -ν , with ξ in Anströms ( Å) and ν = 0.63, the classical Ising value.

As ordering requires only atomic jumps between neighboring sites, the AuAgZn 2 alloy can be considered as a model system to investigate second order transitions with "non conserved order" parameters i.e., the "A" model of Hohenberg and Halperin [1,14]. Experiments are carried out in a small q-range in the vicinity of the 1 Outside the critical region, the driving force of the dynamics is the interface motion between ordered domains. The domain size L increases according to the diffusion law L= √ Dt, which applies for non-conserved dynamics (see e.g., Ref. 15). In a previous experiment [16],

the diffusion coefficient D in AuAgZn 2 was estimated to be 5 × 10 5 Å2 s -1 at T = 607 K.

This led to an elementary diffusion time τ 0 ≃ 20 µs corresponding to atomic jumps between two sites of the Au-Ag simple cubic sublattice (see Fig. 1b). In a quenched sample with antiphase domains, speckles are observed [16,17] and changes in the speckle pattern close to T c reflect the microscopic time evolution of the fluctuations.

In the critical region, equilibrium fluctuations have a characteristic relaxation time τ and, from the dynamic renormalization model [1], τ is related to ξ and q by the following expression:

τ (q, T ) ∝ ξ z Φ(qξ) ∝ { |(T -T c )/T c | -zν for qξ ≪ 1 q -z for qξ ≫ 1 (1)
where the wave vector q corresponds to the distance in the reciprocal lattice to the superstructure Bragg position. Dynamic renormalization group calculations [14], [18] yield z ≃ 2.02. The experimental characterization of the critical fluctuation dynamics is a real challenge since in this system, as explained in Ref. 16, heating devices providing high temperature stability (a few mK) and very good accuracy (∆q < ∼ 10 -4 Å-1 ) are requested. Pioneering experiments on the dynamics of critical fluctuations were attempted by Brauer et al. [19] in the Fe 3 Al binary alloy which exhibits a B2-DO 3 second order transition when the temperature is close to 824 K. In our previous measurements [16], only correlation functions at a fixed temperature and in a narrow q-range could be measured.

The measurement of reliable time correlation functions whose behavior can be compared to Renormalization Group calculations or numerical simulations is still missing.

In this article, some significant improvements in the characterization of critical fluctuation dynamics in an ordering alloy are outlined and the first measurement of the related dynamical critical exponents z is reported.

EXPERIMENTAL DETAILS

Measurements were carried out at the ID10 beamline of the ESRF (France). Undulators were tuned to produce 8 keV photons and a Si(111) monochromator was used to select a wavelenth of 1.55 Å with an intrinsic energy resolution ∆E/E ≃ 1.4 × 10 -4 . The incident beam was focused with Beryllium lenses to have a 3.5 µm (v) × 20 µm (h) section at the sample position and a partially coherent flux of ≃ 2.× 10 10 photons per second. The (111) oriented sample surface was tilted by about 12 • to observe the scattering intensities relative to the 1 2 1 2 1 2 superstructure reflection of the B2 structure, with a vertical diffraction plane. The gauge volume probed by the x-ray beam was thus close to 17µm ×20µm×1µm; the latter value being the x-ray penetration depth into the tilted AuAgZn 2 sample. A 2D pixel detector [20] with 516 × 516 pixels of 55µm ×55µm (2×2 Medipix2 chips) was placed at a distance of 1.73 m from the sample position. The q-resolution of one pixel was thus 1.29 × 10 -4 Å-1 . Measurements were performed in the vertical diffraction plane with a detector angle 2θ B close to 24 • . The results described in the following were obtained by recording sets of 40,000 frames with a maximum time spacing between frames of 0.02 s.

The temperature control of the sample was performed with a Lakeshore controller and a small furnace in vacuum developed to have a maximum deviation of a few mK from a fixed temperature between 553 K and 638 K. The AuAgZn 2 samples were cylinders with a diameter of 8 mm and a height of ≃ 1 -2 mm. Two 70 Å thick aluminum layers were deposited on the sample surfaces to prevent losses of zinc during heating (see detail in Ref. 16). The coherence contrast β of this setup was estimated from the observation of speckles in a quenched sample [16,21] to 3.5%.

Well above T c an isotropic diffuse intensity was recorded at the ordering reflection position. Close to T c , some pretransitional scattering (PTS) was observed, as shown in Fig. 2(a). This scattering is cross-shaped: in the q x direction (horizontal in the figure), some scattering corresponding to surface mosaicity is observed and in the q y direction, the scattering is elongated perpendicular to the surface, corresponding to the small thickness of the surface ordering.

PTS was strongly temperature-dependent even when only a few mK variations of the temperature did occur. For large scattering vectors, the intensity appeared isotropic, but the fluctuation dynamics was too fast and the intensity too low to obtain reliable results.

For all these reasons, only a part of the recorded intensities was utilized for data processing. This was done by using masks like the one represented in Fig. 2(b), where data in black were discarded. After an angular average, the intensity profile of each diffraction pattern was separated into two contributions using [16]:

I(q) =

S(q = 0) 1 + (qξ) 1.97

+ B (1 + (qL) 2 ) 2 (2)
The first term in Eq. 2 is a very good approximation of the static critical scattering S(q) and the second term roughly fits to PTS. PTS is strongly peaking and exhibits some anisotropy even after masking for the smallest q-values (q¡0.0008 Å-1 ) where the fit is poor. This second term is introduced because for intermediate q-values the scattering has a q -4 behavior, which means that PTS corresponds to small well defined ordered domains.

As shown in Fig. 3, the analysis of the temperature behavior of S(q = 0) and ξ in the vicinity of T c led to a reliable estimate of the critical temperature. T c was estimated to 609.162(15) K (or 336.012 • C), and its variations was less than 0.02 K during the 48 hours of the measurements. Our results were in a good agreement with the theoretical values γ=1.241 and ν=0.63 [22].

The dynamics was characterized from the intensity autocorrelation function g 2 (q, t) calculated by a time average and an angular average [16] over circular domains of thickness ∆q

RESULTS

The correlation functions g 2 (q, t)

Fig. 4 represents the values of the time correlation functions g 2 (q, t) obtained for a set of q-vectors at a fixed temperature T c + 0.113. The g 2 (q, t) functions calculated for a given temperature, were fitted by the following expression:

g 2 (q, t) = 1 + a 0 + b 0 e -2t/τ (q) (3)
where τ (q) is the characteristic relaxation time of the critical fluctuations, a 0 is essentially due to the anisotropy of the PTS (see Fig. 2) and, b 0 ≃ βx 2 is connected to the coherence β ≃ 0.035 and to x, the ratio of the critical to PTS intensities. Equation 3 assumes that heterodyne interferences between the two scattering in Eq. 2 were negligible. b 0 increases with increasing values of q: in Fig. 4, its value is in the 0.02-0.004 range. The results of the least-square fitting using Eq. 3 correspond to the solid lines in Figs. 4. Owing to the high experimental stability of the beamline, a reasonable accuracy was obtained for q ≥ 10 -3 Å-1 .

For lower q values, b 0 was too small (b 0 < 10 -3 ) and some experimental issues prevented from obtaining reliable results. For measurements performed at T-T c ¡0.2 K, oscillations of the g 2 (q, t) correlation functions for q < 0.0008 Å-1 were caused by the high sensiblity of the PTS to small temperature instabilities (fluctuations of ±2 mK). Long term (t > 100s)

instabilities of the experiment were also observed.

Critical fluctuation relaxation times τ (q, T ) and τ 0 (T )

The characteristric relaxation times τ (q, T ) in Eq. 3 deduced from the intensity autocorrelation calculated at different temperatures T and q-vectors are plotted in Fig. 5(a). The shortest values of τ close to 0.1 s were measured with a good accuracy well above T c (i. e.

T c + 0.7K) and for large q values (i. e. q ≃ 0.006 Å-1 ). The largest τ values are close to 10 s and were obtained from measurements at small q values (i.e. q ≃ 0.0006 Å-1 )) in the vicinity of T c , in a domain where large PTS was measured. The temperature dependent fluctuation time τ 0 (T ) (for q = 0) was determined using least-squares fits with the following equation:

τ (q, T ) = τ 0 (T )/(1 + (q/q 0 ) 2 ) (4)
Equation 4 assumes a Lorentzian q-dependence of τ (q) (the solid lines in Fig. 5

τ 0 (T ) = ((T -T c )/T c ) νz (5)
with νz=1.235(70). From ν = 0.630 [22], we estimate z ≃ 1.96 (11). This value is in agreement with the best theoretical value of 2.017 [18], but with a poor accuracy.

DISCUSSION

The AuAgZn 2 alloy is an example for the study of the "model A" dynamics of phase transitions, where the order parameter is changed without long range diffusion, analogous to a local spin flip in the Ising model. In order to carry out an XPCS study at large angle, a short x-ray penetration depth is mandatory and the strong contrast between Au and Ag atoms provides a reasonable intensity close to the ordering Bragg peak. This penetration depth (µ -1 ≃ 10µm) leads to small path length differences of the beam in the sample (L ≃ 2.µ -1 sin 2 (θ B ) ≃ 0.9µm is nevertheless significantly larger than the longitudinal coherence length of the experiment (Λ L = λ 2 /(2.δλ) ≃ 0.5µm).

Main difficulty is in dealing with PTS, which makes results in the close vicinity of T c (in Fig. 4, temperatures T < T c + 0.07K were discarded) and for very small angles (q < 1. × 10 -3 Å-1 ) unreliable. PTS is connected to the two 70 Å thick surface layers which were deposited in order to prevent zinc evaporation. While these layers are an excellent chemical barrier, they probably add some constraints that favor surface ordering. As a result, better results in the dynamical study were obtained for regions where the scattering is small. As PTS is connected to surface defects, its interferences with the critical scattering from the inside of the sample is here assumed negligible. For this reason, heterodyning has been neglected, contrary to our previous paper (see 16).

A measurement of the dynamical critical exponent z ≃ 1.96(11) of a 3D Ising system was performed using the XPCS method. The precision of the result does not allow a distinction from mean-field theory (which yields z = 2), and this makes valid the approximations used here: an exponential in Eq.3 and a Lorentzian for Φ(x) in Eq.4.

Obviously, these experiments will be significantly eased with the next ESRF upgrade, where the source brillance should be increased more than one order of magnitude. This will bring a large improvement in the measurement of the dynamic exponent z.

FIG. 1 .

 1 FIG. 1. (a) B2 and (b) L2 1 ordered structures of the AuAgZn 2 alloy near T = 609 K. The lattice parameter a of the B2 phase is close to 3.17 Å at 603 K.

FIG. 2 .

 2 FIG. 2. (a) Measured scattering at 609.35 K (T c + 0.188 K) and (b) the mask used and the typical circular rings where the dynamics was oberved. Both figures show the same 100 × 100 pixel area, each pixel is 1.29 × 10 -4 Å-1

FIG. 3 .

 3 FIG. 3. Determination of the critical temperature by fitting intensity profiles to Eq. 2. (a) from the variation of S -1/γ (q = 0) with γ= 1.241. (b) from the variation of ξ -1/ν with ν = 0.63.

FIG. 4 .

 4 FIG.4. Time correlations obtained for various q vectors at T =T c +0.113K. Curves are shifted vertically for sake of visibility. Symbols correspond to values of the g 2 (q, t) function and solid lines to a least-square fits with Eq. 3

FIG. 5 .

 5 FIG. 5. (a) Fluctuation times τ (q, T ) for various temperatures and the lorentzian fits. (b) Temperature variations of τ 0 (T ). An estimate of the exponent yielded zν = 1.235(70)

  (a)), compatible with mean-field approximations (i. e. z = γ/ν = 2). The double logarithmic plot of τ 0 (T ) shown in Fig. 5(b), roughly corresponds to a scaling law:

The authors want to thank the experimental team of the ID10 beamline of ESRF. * frederic.livet@simap.grenoble-inp.fr


C J Villarraga

J Vaunat
email: jean.vaunat@upc.edu

Muriel Gasc-Barbier

Muriel Gasc
email: muriel.gasc@cerema.fr

MODELLING THERMAL INDUCED DAMAGE IN PERMEABLE ROCKS

Keywords: GiD interface, THM coupled analysis, Geological media

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

INTRODUCTION

Rock degradation due to weather conditions can be a triggering factor for rockfalls, generally, those weathering process are related to freeze-thaw phenomena [START_REF] Hall | Evidence for freeze-thaw events and their implications for rock weathering in northern Canada: II. The temperature at which water freezes in rock[END_REF][START_REF] Matsuoka | Rockfall activity from an alpine cliff during thawing periods[END_REF] . Nevertheless, thermal cycling whit temperatures far from extreme values may lead to the generation and propagation of existent fissures, inducing damage in the rock mass and eventually producing rockfalls, that may occurred without any evident triggering factor, like earthquakes or rains periods. [START_REF] Amaral Vargas | On the effect of thermally induced stresses in failures of some rock slopes in Rio de Janeiro, Brazil[END_REF][START_REF] Vlcko | Rock displacement and thermal expansion study at historic heritage sites in Slovakia[END_REF] The thermal degradation of rocks is mostly influenced by the different thermal properties of the minerals that compose the rock. Differences in thermal expansion, thermal conductivity and principal mineral axis direction, controls the mineral to mineral thermal stress induced by the thermal gradients. Hence, additional stress can be generated leading to failure across the minerals boundaries, and lines of weakness, producing some intergranular disintegration and propagation of existent fissures. [START_REF] Hall | Thermal fatigue and thermal shock in bedrock: An attempt to unravel the geomorphic processes and products[END_REF][START_REF] Walsh | Micromechanical modeling of thermal spallation in granitic rock[END_REF] This environmental thermal induce damage can be analysed under laboratory conditions by the evolution of different physical and mechanical properties, under prescribed thermal cycles. The elastic wave propagation velocity is a common parameter used for this purpose, since, it is extremely sensitive to the presence of defects in the sample, (like the internal fissures) due to the considerable difference of this parameter between the air and the solid materials [START_REF] Brotóns | Temperature influence on the physical and mechanical properties of a porous rock: San Julian's calcarenite[END_REF][START_REF] Martínez-Martínez | Spatial attenuation: The most sensitive ultrasonic parameter for detecting petrographic features and decay processes in carbonate rocks[END_REF] .

Additionally, the elastic wave propagation velocity can easily be related to the dynamic Young and Poisson modulus, allowing the comparison with the results obtained with mechanical test [START_REF] Kendrick | The influence of thermal and cyclic stressing on the strength of rocks from Mount St. Helens, Washington[END_REF] The response of this phenomena is commonly represented with microscopic models, were the mineralogical interaction is analysed [START_REF] Wanne | Bonded-particle modeling of thermally fractured granite[END_REF][START_REF] Shushakova | Microcracking in calcite and dolomite marble: Microstructural influences and effects on properties[END_REF] . For this reason, a numerical model is proposed with the aim of reproducing the macroscale phenomena observed whit an experimental program focused in the thermal fatigue of a Limestone coming from a rock cliff prone to rockfalls.

PROBLEM DESCRIPTION.

La Roque Gageac is a small town located at the south-west of France. It is placed in the border of La Dordogne river at the toe of an 80m high limestone cliff. As presented in Figure 1. This cliff has a troglodyte cavern that is an important tourist attraction to the region.. This cliff is mainly composed by limestones conformed under marine conditions from the intermediate and the upper coniacian. Due to its calcareous nature, the most permeable layer placed in the upper coniacian formation is affected by the karst phenomena. Nevertheless, this cliff is composed by layer of different ages: the upper coniacian and the lower coniacian (an impermeable layer) in the inferior zone. As consequence, the karst phenomena cannot be completely developed in vertical direction.

The town of La Roque Gageac has suffered several rock falls, with a periodicity of more or less 50 years. The last event occurred in January 2010 when a part of the roof of the troglodyte caver collapsed. After this event, several instrumentation systems were installed within and around the troglodyte cavern. Description and first results of this instrumentation program are presented in Ruiz (2013) [START_REF] Ruiz | Thermo-mechanical analysis of the stability of a rock-cliff under climatic actions[END_REF] and Gasc-Barbier (2015) [START_REF] Gasc-Babier | Termal fatigue in rocks-La Roque-Gageac's case study[END_REF] . The response observed with the devices installed show an important relation between deformations in the cliff and natural thermal cycles. As result, a laboratory program has been developed to evaluate the thermal-mechanical rock properties.

EXPERIMENTAL PROGRAM

An experimental program is developed to evaluate the possible damage induced natural thermal variations. For this, samples from La Roque Gageac limestone were submitted to thermal cycles. Since, the aim of this laboratory study is to evaluate the effect of natural thermal variations, the thermal gradient considered is close to the one observed in the cliff, 40ºC. Each cycle last 3 hours and temperature varies from 10ºC to 50ºC.

The degradation of the samples was evaluated through different mechanical parameters, like the elastic wave propagation velocity. The reduction in this parameter can be related to the rock damage, even if, the sample don't show any superficial crack [START_REF] Yavuz | Estimating the index properties of deteriorated carbonate rocks due to freeze-thaw and thermal shock weathering[END_REF] . Moreover, this parameter is obtained whit a non-destructive test, hence, is possible to evaluate its evolution in the same sample through the complete laboratory study.

Samples description

For this experimental program two types of samples were considered: samples drilled from rock blocks dropped within the cavern during the 2010 event and samples drilled directly from the cliff face, nearby the troglodyte cavern.

Even if, both groups of samples were obtained from the same cliff, they show some mineralogical differences associated to the heterogeneity of this geological formation. From mineralogical analysis, it is observed that samples obtained from the blocks are mainly composed by calcite (90% calcite) but, samples form the cliff face are composed by quartz (50%) and calcite (45%)

This mineralogical composition difference may be the responsible of differences observed in the initial elastic wave propagation velocity Samples from the rock blocks have a mean compression elastic wave propagation velocity (VP) of 4400 m/s while samples from the cliff face 2900m/s

Elastic wave propagation velocities results.

The elastic compression (Vp) and shear (Vs) wave propagation velocities are obtained by the ultrasonic pulse technique following the French Standard NF P 94-411. From the results of elastic compression (Vp) and shear (Vs) wave propagation velocities, the dynamic Young (E) and Poisson () modulus can be calculated [START_REF] Kendrick | The influence of thermal and cyclic stressing on the strength of rocks from Mount St. Helens, Washington[END_REF] , following the next equation:

        ) 2 1 ( 3  E K Where,               ) 1 ( ) 2 1 )( 1 ( 2     p V E and                                      2 1 2 5 . 0 p V s V p V s V   (1)
After 500 cycles of temperature a reduction in Vp and Vs values is observed. This reduction represents a stiffness reduction around 6% for the calcite mainly composed samples and 18% for quartz-calcite composed samples, as example Figure 2 shows the evolution of Vp and Stiffness (K) for sample the 5 obtained from the rock blocks and the sample 27 taken from the face cliff.

A)

B) From de mineralogical point of view the thermal induce damage in rocks can be considered as a microstructural re-accommodation associated to the different thermo-mechanical response of the minerals that compose the rock. Nevertheless, the aim of this study is to be available of modelling this phenomenon at a macroscopic scale, not at mineral scale.

For this, the composite material model purposed by Vaunat and Gens (2003) [START_REF] Vaunat | Bond degradation and irreversible strains in soft argillaceous rock Degradación y deformación irreversible de una roca blanda arcillosa[END_REF] for argillaceous rocks is considered, this original model was modify in order to fulfilled the characteristics of the phenomena under study.

The model considers that the rock is composed by two different components: a strong matrix and bonds. As illustrated in Figure 3. For the case of La Roque Gageac, the matrix represents the sandy minerals, and the bonds, the calcite minerals. The mechanical response of the matrix may be defined by different failure criteria. In the other hand, the bond behaviour is defined by the scalar damage model proposed by Carol et al (2001) [START_REF] Carol | On the formulation of anisotropic elastic degradation. I. Theory based on a pseudo-logarithmic damage tensor rate[END_REF] , Its response is elastic, but with moduli that progressively decrease during loading as microcracks developed inside the material.

Figure 3 Schematic rock composition

From the original composite model, the principal difference introduced is that the two components of the material are considered to work in parallel; therefore, the new strain partition is defined as:

(2)

Where is the volumetric strain increment applied to the medium, the volumetric strain used to update soil porosity the volumetric strain inside bonds and the volumetric strain inside the matrix, Cb is the bond content and CM the matrix content. Two proportional coefficients between , and are defined, these coefficients are model parameters:

(3)

The stress partition is obtained by stating that the work input into the medium under any compatible external strain increment is equal to sum of the work input into the matrix and the bonds. Finally, the stress partition is defined by the equation 4. The composite material model previously described was implemented in the finite element code CODE_BRIGHT [START_REF] Olivella | Numerical formulation for a simulator (CODE_BRIGHT) for the coupled analysis of saline media[END_REF] . This model is used to reproduce the experimental results obtained with samples from La Roque Gageac subjected to thermal cycles between 10ºC to 50ºC. For instance, the responses of samples 5 and 27 (presented in section 3), are modelled at Gausspoint level. Table 1 presents the damage parameter considered, the mechanical behaviour of the matrix is defined by a Mohr-Coulomb failure criteria. Note that the bonding content (Cb) consider for the modelling is different for each sample, this is related to the different mineralogical composition of the samples considered. Sample 5 is mainly composed by calcite. Therefore, the predominant component is the bond, nevertheless, sample 27 is composed in a 50-50 relation by quartz and calcite, then its bonding content is lower that for the 90% calcite sample.

Parameter

Figure 4 presents the results obtained in the numerical modelling. In both cases the model is capable of reproducing the stiffness reduction observed whit the laboratory tests. Note also, that sample 5 has a slower damage rate that the one observed in sample 27, as is experienced whit the laboratory analysis. The thermal degradation of the rocks, induce the generation of micro-fissures but also the propagation of existent ones. This last phenomenon could be evaluated whit the numerical model presented. If we consider that the pre-existing fissure is a weaker zone whit lower bond damage locus stress. The response of this zone may increment the stress of the surrounding material inducing more damage.

A) B)

Figure 1

 1 Figure 1 La Roque Gageac localization

Figure 2

 2 Figure 2 Laboratory results for samples 5 and 27 A) Vp evolution B) Stiffness evolution

Figure 4

 4 Figure 4 Numerical model results A) sample 5 B) sample 27

Table 1 Damage model parameters

 1

		Sample 5	Sample 27
	Matrix stiffness Ku (Mpa)	15900	15900
	Initial bonding stiffness KB0 (MPa)	18000	15200
	Porosity (n)	0.15	0.12
	Matrix thermal expansion M (1/ºC)	6e-6	6e-6
	Bond thermal expansion b (1/ºC)	6e-6	6e-6
	Bond damage locus rb0 (MPa)	1e-5	1e-5
	Damage evolution rate r1 (MPa)	0.005	0.005
	CB	0.7	0.5
	XM	3.5	2.5
	XB	5.5	5.5