Alexandra Rolland

Renaud Toussaint

Patrick Baud

Jean Schmittbuhl

Nathalie Conil

Daniel Koehn

François Renard

Jean-Pierre Gratier

Modeling the growth of stylolites in sedimentary rocks

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

Introduction

Stylolites are undulated surfaces resulting from localized stress-driven dissolution of some minerals of the rock. Insoluble minerals as clay particles, oxides and organic matters are concentrated in the interface and make stylolites visible. [START_REF] Bathurst | Diagenetically enhanced bedding in argillaceous platform limestones: stratified cementation and selective compaction[END_REF] describes stylolites as serrated interfaces with an amplitude greater than the diameter of the transected grains giving them a sutured appearance. He makes a difference with dissolution seams or 'flaser' which are smooth, undulating, lacking in sutures and fitting around grains instead of cutting through them. Stylolites are most often found in carbonates [START_REF] Stockdale | Stylolites: their nature and origin[END_REF][START_REF] Stockdale | The stratigraphic significance of solution in rocks[END_REF][START_REF] Stockdale | Rare stylolites[END_REF][START_REF] Stockdale | Stylolites: primary or secondary?[END_REF][START_REF] Dunnington | Stylolite development post-dates rock induration[END_REF][START_REF] Bushinskiy | Stylolites[END_REF][START_REF] Park | Stylolites: their nature and origin[END_REF][START_REF] Bathurst | Carbonate sediments and their diagenesis[END_REF][START_REF] Buxton | Pressure solution features in a shallow buried limestone[END_REF][START_REF] Railsback | Lithologic controls on morphology of pressure-dissolution surfaces (stylolites and dissolution seams) in Paleozoic carbonate rocks from the mideastern United States[END_REF] but also in sandstones [START_REF] Young | Stylolitic solution in Wit-watersrand quartzites[END_REF][START_REF] Heald | Stylolites in sandstones[END_REF], shales [START_REF] Wright | Pressure dissolution and cleavage in Martingsburg shale[END_REF][START_REF] Rutter | Pressure solution in nature, theory and experiment[END_REF], cherts [START_REF] Bushinskiy | Stylolites[END_REF][START_REF] Iijima | Nature and origin of the Paleogene cherts in the Setogawa Terrain, Shizuoka, central Japan[END_REF][START_REF] Cox | Stylolites in the Caballos Novaculite[END_REF]] and some-times in coal [START_REF] Stutzer | Geology of coal[END_REF]. Stylolites are divided in two groups according (i) to their orientation with respect to the bedding of the surrounding rock or (ii) to the orientation of their 'tooth' with respect to the mean plane of the stylolite. The first group shows two types of orientation: stylolites parallel to the bedding plane, designated as sedimentary, and formed under the lithostatic pressure and stylolites oblique or even perpendicular to the bedding, designated as tectonic, and depending on the maximum tectonic stress. The tooth orientation is in both cases an indicator of the direction of the incremental displacement which is parallel to the major principal stress in co-axial deformation. The stylolites of the second group are called 'slickolites' [Ebner et al., 2010a]. They develop when there is a preferential plane for their growth (bedding or fracture). In this case, the stress is not perpendicular to the mean plane of the stylolite [START_REF] Stockdale | Stylolites: their nature and origin[END_REF], but the edges of the tooth are subparallel to the maximum principal stress axis. Various studies [START_REF] Park | Stylolites: their nature and origin[END_REF][START_REF] Renard | Pressure solution in sandstones: influence of clays and dependence on temperature and stress[END_REF][START_REF] Renard | Enhanced pressure solution creep rates induced by clay particles: experimental evidence in salt aggregates[END_REF][START_REF] André | Caractérisation des déformations mésocénozoïques et des circulations de fluides dans l'Est du Bassin de Paris[END_REF][START_REF] Aharonov | Interaction between pressure solution and clays in stylolite development: insights from modeling[END_REF] suggest that many parameters play an important role in the stylolite growth such as confining pressure, deviatoric stress, fluid pressure, temperature, shape and assemblage of grains, anisotropy of minerals, rates of dissolution and presence of clay (acting potentially as catalyst for the dissolution).

Only few papers report experiments about stylolites development. Indeed, they are inherently difficult to reproduce as the kinetics of pressure-solution processes is very slow [START_REF] Rutter | The kinetics of rock deformation by pressure solution[END_REF]. Experiments were conducted either on aggregates [START_REF] Cox | Experimental dissolutionprecipitation creep in quartz aggregates at high temperatures[END_REF][START_REF] Brok | The effect of elastic strain on the microstructure of free surfaces of stressed minerals in contact with an aqueous solution[END_REF][START_REF] Renard | Enhanced pressure solution creep rates induced by clay particles: experimental evidence in salt aggregates[END_REF][START_REF] Gratier | Experimental microstylolites in quartz and modeled application to natural stylolitic structures[END_REF] or with indenter techniques [START_REF] Gratier | Experimental pressure solution-deposition on quartz grains: the crucial effect of the nature of the fluid[END_REF][START_REF] Gratier | Experimental pressure solution of halite by an indenter technique[END_REF][START_REF] Gratier | Indenter studies of the swelling, creep and pressure solution of Bure argillite[END_REF][START_REF] Dysthe | Universal scaling in transient creep[END_REF][START_REF] Dysthe | High resolution measurements of pressure solution creep[END_REF][START_REF] Karcz | Deformation by dissolution and plastic flow of a single crystal sodium chloride indenter: an experimental study under the confocal microscope[END_REF]. [START_REF] Dysthe | Universal scaling in transient creep[END_REF][START_REF] Dysthe | High resolution measurements of pressure solution creep[END_REF] used an indenter technique where a sodium chloride crystal was kept in contact with a piston at given pressure and temperature for several months. A fluid at compositional equilibrium with the crystal is trapped between the sample and the indenter. The contact evolved due to pressure-solution during the indentation. A power law time dependence with an exponent value of 1/3 as in Andrade creep law was shown to control the indentation rate. The observed microstructures in the contact seem to be different from stylolites. [START_REF] Karcz | Deformation by dissolution and plastic flow of a single crystal sodium chloride indenter: an experimental study under the confocal microscope[END_REF] loaded a halite cone-shaped indenter against a flat silicate surface immersed in an undersaturated brine. Using confocal microscopy techniques, they observed that the evolution of the system is dictated by an interaction between two deformation mechanisms: undercutting dissolution reducing the area of the contact and plastic flow increasing it. Recently, similar experiments were carried out with a brine at chemical equilibrium with the crystal [START_REF] Ben-Itzhak | Pressure solution and stylolites in carbonate rocks[END_REF]Laronne Ben-Itzhak et al., 2012a]. Emerging evolving islands and channels were observed at the contact. Such islands and channels structures were previously observed at the contacts during experiments on aggregates [START_REF] Schutjens | Intergranular pressure solution in NaCl: Grain-to-grain contact experiments under the optical microscope[END_REF][START_REF] Brok | Experimental evidence for water weakening of quartzite by microcracking plus solution precipitation creep[END_REF]. Other experiments on aggregates were performed by [START_REF] Gratier | Experimental microstylolites in quartz and modeled application to natural stylolitic structures[END_REF]. They loaded layers of fine quartz sand grains. The experiments lasted several months at 350 • C, under 50 MPa of differential stress and in presence of an aqueous silica solution. Microstylolites were created for the first time in the laboratory at the stressed contacts between the quartz grains. An interesting observation is that the stylolites peaks are always located in front of dislocation pits. Consequently, stylolites appear to be localized by the heterogeneities of the mineral. Den [START_REF] Brok | The effect of elastic strain on the microstructure of free surfaces of stressed minerals in contact with an aqueous solution[END_REF] loaded elastically K-alum crystals at a controled temperature and in a saturated K-alum solution. A hole was drilled in the middle of the crystals to provide an elastic strain gradient. They observed macroscopic etch grooves on the originally smooth free surfaces of the soluble crystals which disappear when removing the stress. [START_REF] Koehn | Transient dissolution patterns on stressed crystal surfaces[END_REF] stressed crystals of N aClO3 in a N aClO3 solution at room temperature. Parallel dissolution grooves developed on their free surface in a 1D geometry to a 2D geometry with the coarsening of the pattern. The pressure-solution process slowed down or stopped progressively with the increasing concentration of the solution during the experiments. [START_REF] Gratier | Indenter studies of the swelling, creep and pressure solution of Bure argillite[END_REF] used a similar technique in which a sample of Bure claystone was kept in contact with a piston, with a saturated brine in the contact, at an imposed pressure and temperature for several months. No evidence of localized pressure-solution (dissolution seam) was observed in this case, grain to grain sliding being more efficient in presence of clay. [START_REF] Renard | Enhanced pressure solution creep rates induced by clay particles: experimental evidence in salt aggregates[END_REF], studied chemical compaction of aggregates of halite (salt) mixed with clay. They showed that clay particles enhance pressure-solution. Moreover, [START_REF] Renard | Pressure solution in sandstones: influence of clays and dependence on temperature and stress[END_REF] studied the effect of clay on clay-rich sandstones. They suggested that pressure-solution is enhanced by clay because a thick film of water is preserved between clay particles. They also concluded that the depth determines the limiting factor for the process: at great depth, the water film between grains should be thinner and diffusion limits the process. Conversely, at low depth water films are bigger, transport is easier and the reaction kinetics is the limiting factor.

The clay particles effect on pressure-solution was recently simulated in numerical modeling. [START_REF] Aharonov | Interaction between pressure solution and clays in stylolite development: insights from modeling[END_REF] used the two-dimensional Spring Network Model to study the stylolites growth in a medium with a uniform clay distribution. They showed that clay plays a role of enhancing pressure-solution and that stylolites propagation is possible only when both pressure-solution and clay-enhanced dissolution operate together. [START_REF] Koehn | Growth of stylolite teeth patterns depending on normal stress and finite compaction[END_REF] developed a new discrete simulation technique that reproduces successfully the roughening of stylolites from a preferential existing surface with no clay. This model is based on molecular dynamics, with a dissolution speed depending on the local free energy that includes stress dependent terms and surface energy terms. Two different spatial regimes arise from this modeling: a small-scale regime where surface energy is dominant with significant fluctuations of the roughness and a large-scale regime where elastic energy dominates. The dependence on the cross-over scale between both regime on the imposed stress has been recently investigated numerically [START_REF] Koehn | Modelling of stylolite geometries and stress scaling[END_REF]. This model shows that the growth of the stylolite tooth follows the main compressive stress direction. The nature and structure of the small scale disorder for the dissolution properties of grains were systematically analyzed [Ebner et al., 2009a]. Moreover, Ebner et al. [2010b] performed detailed microstructural analysis to investigate the interplay between this disorder and the compositional nature of the grains surrounding a stylolite.

Stylolites are localised features for which deformation is purely compactant as for compaction bands [START_REF] Mollema | Compaction bands: a structural analog for anti-mode I crack in aeolian sandstone[END_REF][START_REF] Baud | Compaction localization in porous sandstones: spatial evolution of damage and acoustic emission activity[END_REF]Katsman et al., 2006b;[START_REF] Tembe | Stress conditions for the propagation of discrete compaction bands in porous sandstone[END_REF]. Stylolites and compaction bands development was modeled as anticracks or anti-mode I fracture [START_REF] Fletcher | Anticrack model for pressure solution surfaces[END_REF][START_REF] Rispoli | Stress field about strike-slip faults inferred from stylolites and tension gashes[END_REF][START_REF] Mollema | Compaction bands: a structural analog for anti-mode I crack in aeolian sandstone[END_REF]. [START_REF] Fletcher | Anticrack model for pressure solution surfaces[END_REF] assume that the rate of pressure-solution is only a function of the normal stress. They observed an elliptic dissolution pattern i.e. more dissolution in the central part of stylolites than at the tips. With these observations they proposed an analogy between propagation of stylolites and propagation of mode I fractures. They observed that the relative displacement between the sides of a stylolite should have the opposite sign than that of a crack, and thus termed their model an anticrack. Note however that cracks can bear zero surface traction, contrary to stylolites. This distinction between crack solutions and stylolites was introduced, and it was shown by Katsman et al. [2006a] that, as compaction bands, stylolites are Localized Volume Reduction zones (LVR). The shape of the displacement along stylolites, and how the stress perturbation can be determined from the concept of LVR, is discussed in details by [START_REF] Katsman | Extensional veins induced by self-similar dissolution at stylolites: analytical modeling[END_REF].

In LVR where the dissolution amount is constant across the surface of the LVR, as for a compaction band, the stress enhancement was shown to be that of a dislocation [Katsman et al., 2006a]. In later models [START_REF] Katsman | Extensional veins induced by self-similar dissolution at stylolites: analytical modeling[END_REF], it was shown that if more dissolution is allowed in the center of a stylolite, another type of stress enhancement, with a dependence on the distance to the tips analogous to the one for a crack (rather than to a dislocation), can be observed it is given by the Eshelby inclusion problem. Such an increase of the dissolution in the center of a stylolite, where the dissolution does not stop in the already dissolved zone in the middle of the stylolite, can be observed in models with a positive feedback do the dissolution, as for example the one that can be modeled from a clay concentration mechanism [START_REF] Aharonov | Interaction between pressure solution and clays in stylolite development: insights from modeling[END_REF].

In general, in stylolites, the stress concentrates at the tips and the largest stress is perpendicular to the stylolites. Recent models [START_REF] Koehn | Growth of stylolite teeth patterns depending on normal stress and finite compaction[END_REF]Ebner et al., 2009b;[START_REF] Zhou | Mechanics of pressure solution seam growth and evolution[END_REF] suggest that a higher stress concentration at the top of the tooth should be responsible of localized high rates of dissolution. [START_REF] Benedicto | Stylolites in limestones: magnitude of contractional strain accommodated and scaling relationships[END_REF] investigated the topography of stylolites (along-strike trace length, maximum and average amplitudes) from the damaged zone of the Gubbio normal fault zone in central Italy. They showed that the amount of contractional strain accommodated by stylolites as well as their length and their number increase according to the topography parameters. Analyses of cores from boreholes reveal also an increase in X -3 stylolite abundance with depth [START_REF] Lind | Stylolites in chalk from leg 130, Ontong Java Plateau[END_REF]. [START_REF] Fabricius | Stylolites, porosity, depositional texture, and silicates in chalk facies sediments. Ontong Java Plateau -Gorm and Tyra fields[END_REF] compared formations of chalk from boreholes on the Ontong Java Plateau and in the central North Sea. They showed that the burial stress and the temperature play distinct roles in the burial diagenesis and porosity development of chalk. Pressure-solution and physical compaction are controled by the burial stress while the temperature controls recrystallization and cementation. Moreover, [START_REF] Lind | Stylolites in chalk from leg 130, Ontong Java Plateau[END_REF] suggests that mineralogical anomaly is an initializing factor in stylolite formation such as burrows, shale clasts or flaser structures. Many studies were conducted on the morphology of sedimentary stylolites [START_REF] Renard | Three-dimensional roughness of stylolites in limestones[END_REF][START_REF] Brouste | Variety of stylolites' morphologies and statistical characterization of the amount of heterogeneities in the rock[END_REF]Ebner et al., 2009b]. Morphology analyses can be done on 1D profiles or 2D opened surfaces. They consist on studying a stylolitic profile or surface height variations (standard deviation, height differences, power spectrum, average wavelet coefficient spectrum, etc.) over different scales [START_REF] Schmittbuhl | Scaling invariance of crack surfaces[END_REF][START_REF] Schmittbuhl | Roughness of stylolites: implications of 3D high resolution topography measurements[END_REF][START_REF] Renard | Three-dimensional roughness of stylolites in limestones[END_REF]. These analyses reveal two distinct scaling regimes that could be described by power laws. The power laws are function of a roughness exponent also called Hurst exponent inferred to be 1 and 0.5 for small and large scale respectively [START_REF] Renard | Three-dimensional roughness of stylolites in limestones[END_REF][START_REF] Schmittbuhl | Roughness of stylolites: implications of 3D high resolution topography measurements[END_REF][START_REF] Brouste | Variety of stylolites' morphologies and statistical characterization of the amount of heterogeneities in the rock[END_REF]Ebner et al., 2009b]. The two regimes are separated by a cross-over length typically around 1 mm [START_REF] Renard | Three-dimensional roughness of stylolites in limestones[END_REF][START_REF] Schmittbuhl | Roughness of stylolites: implications of 3D high resolution topography measurements[END_REF]. For sedimentary stylolites, the two dimensional (2D) analysis of their surface does not show any significant inplane anisotropy reflecting the fact that horizontal stresses are isotropic. Ebner et al. [2010a] observed that the profiles of tectonic stylolites show the same geometric attributes as sedimentary ones. Two different regimes are also observed with Hurst exponent around 1 and 0.5 for small and large scale respectively. However, for tectonic stylolites, the 2D analysis revealed an anisotropy of the cross-over length which varies with the direction in the plane of stylolites. Ebner et al. [2010a] argue that this anisotropy develops because the stylolite roughens in an anisotropic inplane stress field. The vertical and inplane horizontal stresses are significantly differents. In recent papers, stylolites are presented as fossilized signatures of the stress field [START_REF] Renard | Three-dimensional roughness of stylolites in limestones[END_REF][START_REF] Schmittbuhl | Roughness of stylolites: implications of 3D high resolution topography measurements[END_REF]Ebner et al., 2009bEbner et al., , 2010a]]. The existence of two scaling regimes for sedimentary stylolites was shown in [START_REF] Schmittbuhl | Roughness of stylolites: implications of 3D high resolution topography measurements[END_REF] where a brief theoretical derivation was performed. It was shown that the cross-over length between both scaling regimes is expected to be dependent on the stress acting on the stylolite during its growth. Their conclusion was that stylolite morphology can be used as a paleostress magnitude indicator. This conclusion was later probed independently on two types of approaches: first, on field data sampled from the same formation at different heights, Ebner et al. [2009b] showed that the measured cross-over length in the morphology followed the expected scaling with the burial stress, evaluated from the position in the formation. Next, discrete numerical simulations were carried out at different stress magnitudes, allowing for the dissolution of grains along the fluid/rock interface, with free energy depending on interfacial tension and local stress. It was shown that the two expected scaling regimes were observed [START_REF] Koehn | Growth of stylolite teeth patterns depending on normal stress and finite compaction[END_REF][START_REF] Koehn | Modelling of stylolite geometries and stress scaling[END_REF], and that the crossover length followed the predicted dependence on the farfield stress amplitude [START_REF] Koehn | Modelling of stylolite geometries and stress scaling[END_REF].

Interfaces between solids and fluids are related to models of stylolitization. In the case where a solid in contact with a fluid is stressed, an instability due to pressure-solution was shown theoretically to exist and is called the Asaro-Tiller-Grinfeld (ATG) instability [START_REF] Renard | Three-dimensional roughness of stylolites in limestones[END_REF]. In models of dissolving surfaces with a stress imposed to a solid in contact with a fluid at chemical equilibrium, this instability leads to the growth of initial large scale modulations of the surface with a wavelength selection obtained through a fastest growing mode. The basic equation depends on the particular boundary conditions e.g. when two solids with different elastic properties are in contact and submitted to a stress, the interface can undergo a fingering instability led by the contrast between the free energies applied to both solids [START_REF] Angheluta | Stress-driven phase transformation and the roughening of solid-solid interfaces[END_REF][START_REF] Angheluta | Thermodynamics and roughening of solid-solid interfaces[END_REF][START_REF] Angheluta | Morphological instabilities of stressed and reactive geological interfaces[END_REF]. The stability analysis can be performed theoretically from expressions for the kinetics using local free energy criteria for the reaction rate [START_REF] Renard | Three-dimensional roughness of stylolites in limestones[END_REF][START_REF] Schmittbuhl | Roughness of stylolites: implications of 3D high resolution topography measurements[END_REF], or global ones [START_REF] Bonnetier | Does roughening of rock-fluid-rock interfaces emerge from a stress-induced instability?[END_REF][START_REF] Angheluta | Stress-driven phase transformation and the roughening of solid-solid interfaces[END_REF]. Depending on the boundary conditions, this situation is also found to be unstable for perturbations exceeding a certain wavelength, leading to fingering (as e.g. with large stress tangential to a fluid interface, or a stress normal to fluid interfaces and lateral periodic boundary conditions [START_REF] Bonnetier | Does roughening of rock-fluid-rock interfaces emerge from a stress-induced instability?[END_REF]). With other boundary conditions, the surface energy and elastic interactions are found to stabilize the interfaces, which are only destabilized by material noise due to heterogeneities [START_REF] Schmittbuhl | Roughness of stylolites: implications of 3D high resolution topography measurements[END_REF][START_REF] Koehn | Growth of stylolite teeth patterns depending on normal stress and finite compaction[END_REF]. We will argue in details in the discussion section about the different possibilities applied to the geometry of stylolites, and the fact that stylolites displaying self-affine scaling laws for their height at large scale are compatible with the stabilizing character of elastic forces at large scale. This manuscript provides the technical development and details that lead to the final result that was previously published without derivation, in a condensed form [START_REF] Schmittbuhl | Roughness of stylolites: implications of 3D high resolution topography measurements[END_REF]. It also compares the result of the analytical development to a direct numerical simulation.

In this paper we concentrate on the following questions: (i) Is the elastic energy stabilizing or destabilizing? (ii) What is the significance of the obtained paleostress values? To answer to these questions, (i) we derive the details of the computation leading to the link between the paleostress magnitude and the cross-over length between the two scaling regimes. This is performed by a perturbative analysis of the elastic energy around an interface slightly wavy and unaligned with one of the principal stresses. Then we show in details that the mechanics and chemistry allow to relate the small and large scale behavior of stylolites to known models, with Hurst exponents corresponding to the observed ones.

(ii) We finally present and discuss an application in relation with the geological context. This is made on a stylolite from the Bure carbonates and it shows how the predicted scaling regimes can be found, and how to determinate the paleostress from the extracted cross-over length.

Analytical approach : Continuous elastostatic model for stylolite propagation

The rough morphology of stylolites arises from the disorder present in a rock and its impact on the pressure-solution process. This disorder is spatially linked to the grains constituting the rock. To understand the impact of this disorder on the chemico-mechanical coupling, we will consider the following simplified geometry: the initial stage of the stylolite is modeled as an elongated fluid pocket enclosed between two contactless rough surfaces of infinite extent. The contacts between these two surfaces can in principle modify the geometry of the resulting dissolution surface. However, they are assumed to be sufficiently loose in a stylolite and thus the main morphological results are not affected. This assumption simplifies the problem since the dissolution process, happening on both sides of the stylolite (Figure 1), can be described as the dissolution of a solid half-plane in contact with a fluid. With this geometry, the small and large scale self-affine behaviors of the dissolution surface and the associated roughness exponents (or Hurst exponents) are well reproduced. The model leads to the characteristic exponents typically observed in previous studies [START_REF] Renard | Three-dimensional roughness of stylolites in limestones[END_REF][START_REF] Schmittbuhl | Roughness of stylolites: implications of 3D high resolution topography measurements[END_REF][START_REF] Brouste | Variety of stylolites' morphologies and statistical characterization of the amount of heterogeneities in the rock[END_REF]Ebner et al., 2009b].

The average stylolitic plane is defined along the x and y-axis (Figure 2). To have better statistics on the morphology of the studied surfaces, the model is assumed to be invariant by translation along the y-axis. It allows to us to describe a larger range of scales at the same numerical cost and to numerically solve the self-affine behavior of the resulting pressure-solution surface over a larger number of orders of length scales. The same approach can be considered using invariance by translation along the x-axis. In the model we assume a mechanical equilibrium throughout the system and express the dissolution rate as a function of the stress tensor and of the area of interface per unit volume.

Force perturbation related to the mechanical equilibrium along the fluid-solid interface

First, we express the mechanical equilibrium at the solidfluid interface (Figure 2). The convention adopted is that compressive stresses and compactive strains are negative [START_REF] Landau | Theory of elasticity[END_REF]. The far-field stress applied to the host rock is denoted by σ0 . The largest principal stress axis, perpendicular to the average plane of the stylolite, is defined along the z-axis. The fluid pocket transmits all the load through itself (The boundary condition of the fluid pocket is approximated as undrained for that respect: if there is any flow, from or into the fluid pocket, it happens slowly, via the lateral ends. If there is any contact between the opposite walls perpendicular to the main fluid direction, the load transmitted through this contact is neglected). The fluid pressure is thus homogeneous and equal to the largest principal stress applied to the host rock, considering the integral of the local stress field σ along an elongated rectangular boundary (dashed line in Figure 1):

p = -σ 0 zz (1)
Locally, the local stress σ is split between the far-field asymptotic value σ0 and a perturbation generated by the irregular nature of the interface σ1 :

σ(x) = σ0 + σ1 (x)
(2)

The far-field stress unit vectors x and ẑ along the x and z-axis are assumed to be the principal directions i.e.

σ0 = σ 0 xx xx + σ 0 zz ẑ ẑ (3)
Here, the notations xx and ẑẑ correspond to unit matrixes composed from the unit vectors, as e.g. are ŷ ŷ, xŷ, or ẑ x. This canonical basis for the matrixes is composed from the doublets of unit vectors x, ŷ and ẑ. For example, xẑ represents the unit matrix with all components equal to zero, apart from a unit in the lign corresponding to the x coordinate, and the column corresponding to the z one, so that for a pair of vectors u, v applied to the left and right of this matrix, u

• (xẑ) • v = (u • x)(ẑ • v) = uxvz.
In other terms, with cartesian components along directions of indexes i and j, and the help of the Kronecker symbol δ, the components of the matrix xẑ, for example, are: (xẑ)ij = δixδjz. This convention to define the canonical basis of matrix space (nine elementary second order dyadic products like xẑ) from the three basic unitary vectors of the vectorial space, x, ŷ and ẑ is, for example, defined by [START_REF] Gonzalez | A first course in continuum mechanics[END_REF].

For a stylolite, the largest compressive stress axis is normal to its average plane and thus to the average fluid pocket direction:

|σ 0 zz | > |σ 0 xx | (4)
This relation has strong implications on the stability of the surface pattern emerging from the dissolution process. The far-field deviatoric stress is defined as:

σ 0 s = (|σ 0 zz | -|σ 0 xx |) = (σ 0 xx -σ 0 zz ) (5)
To express the force perturbation related to the curved nature of the interface, we define the unit vector n normal to the surface pointing towards the fluid. This vector is assumed to be close to the principal stress axis. In the following, we will consider small-angle deviations from a straight surface, and the results will therefore be valid for small surface slopes only. The model presented below aims to describe the onset of the stylolite propagation from a flat surface, and it will also describe the evolution of large wavelength modes, if the aspect ratio of such modes (ratio of the amplitude over the wavelength) stays small, corresponding to small effective slopes at large wavelength. The interface is described as a single-valued function z(x) and the slopes are assumed to be of the order ǫ i.e. that

X -5

|∂x(z)| ∈ O(ǫ) ≪ 1. Since the normal n to the interface of slope ∂xz can be expressed by the conditions of normality to the interface, n • (1, ∂xz) T = 0 (at any order or ǫ), and by its unitary norm n2 = 1, it is in general n = (-(∂xz)x + ẑ)/ 1 + (∂xz) 2 . Using the above limit of small slopes, developing in ǫ, we obtain to leading order

n = ẑ -(∂xz)x + O(ǫ 2 ) (6) (The order O(ǫ) is absent from n).
The local mechanical equilibrium at the solid-fluid interface is expressed as:

σ • n = -pn (7)
And with equations (1-7) the force perturbation (illustrated in Figure 3) becomes:

δf (x) = σ 1 (x) • n = -pn -σ 0 • n = σ 0 zz [ẑ -(∂xz)x] -σ 0 zz ẑ + σ 0 xx (∂xz)x = (σ 0 xx -σ 0 zz )(∂xz)x = σ 0 s (∂xz)x (8)

Chemico-mechanical coupling

Next, we express the chemico-mechanical coupling. The dissolution speed normal to the solid/fluid interface (in mol.m -2 .s -1 ), is to the first order proportional to the chemical potential ∆µ of the chemical product dissolving [START_REF] Kassner | Phase-field modeling of stress-induced instabilities[END_REF][START_REF] Misbah | Dynamics of a dissolution front for solids under stress[END_REF][START_REF] Schmittbuhl | Roughness of stylolites: implications of 3D high resolution topography measurements[END_REF][START_REF] Koehn | Growth of stylolite teeth patterns depending on normal stress and finite compaction[END_REF]:

v = m∆µ (9)
where

m = k0Ω/RT (10)
is the mobility of the dissolving species, R = 8.31 J.mol -1 .K -1 is the universal gas constant, T is the temperature in Kelvin, k0 is a dissolution rate which can be measured experimentally, and Ω is a molar volume. For calcite, Ω ≃ 4 • 10 -5 m 3 .mol -1 and k0 ≃ 10 -4 mol.m -2 .s -1 for dissolution in water at atmospheric pressure and 298 • K [START_REF] De Giudici | Surface control vs. diffusion control during calcite dissolution: Dependence of step-edge velocity upon solution pH[END_REF][START_REF] Schmittbuhl | Roughness of stylolites: implications of 3D high resolution topography measurements[END_REF]. The difference in chemical potential from the solid state to the fluid state is [START_REF] Kassner | Phase-field modeling of stress-induced instabilities[END_REF][START_REF] Misbah | Dynamics of a dissolution front for solids under stress[END_REF][START_REF] Koehn | Growth of stylolite teeth patterns depending on normal stress and finite compaction[END_REF]:

∆µ = ∆Ψs + Ω∆Pn + Ωγκ (11)
Considering a solid state at given pressure and elastic free energy in chemical equilibrium with the fluid, ∆Ψs and ∆Pn are defined respectively as the change in Helmoltz free energy per mole and the change in stress normal to the interface. The last term corresponds to the surface energy with κ = ∂xxz, the surface curvature (the inverse of the radius curvature) and γ the surface tension between the solid and the fluid phase. In a particular case, neglecting temperature variation effects and assuming that the fluid composition is in chemical equilibrium with a solid flat surface at normal pressure p and stress σ ref , Eq.( 11) reduces to :

∆µ = 0 (12) κ = 0 (13)
More generally, by definition ( [START_REF] Kassner | Phase-field modeling of stress-induced instabilities[END_REF]:

∆Ψs + Ω∆Pn = Ω∆ue, (14)
= [(1 + ν)σij σij -νσ kk σ ll ]/4E (16)
is the elastic free energy per unit volume with E the Young's modulus and ν the Poisson's ratio of the elastic solid [START_REF] Kassner | Phase-field modeling of stress-induced instabilities[END_REF][START_REF] Landau | Theory of elasticity[END_REF].

To take into account the dissolution speed variations associated to the morphology of the stylolite, we develop the dissolution speed to the leading order as:

v = v 0 + v 1 (17)
With equations (9-16):

v 0 = k0Ω 2 RT [(1 + ν)σ 0 ij σ 0 ij -νσ 0 kk σ 0 ll ] 4E -u ref e = k0Ω 2 RT E (αp 2 0 -α ref p 2 ref ) (18)
The geometrical factor α is computed assuming σ 0 xx = σ 0 yy = -p0 + σs/3 and σ 0 zz = -p0 -2σs/3:

α = 9(1 -2ν) + 2(1 + ν)σ 2 s /p 2 0 12 ( 19
)
α ref is expressed with a similar expression and characterizes the chemical equilibrium with the fluid at the referential state as a function of the pressure p ref and the shear stress

σ ref .
Typically, for a limestone with a Young's modulus E = 80 GPa stressed at p0 ≃ 10 MPa (which corresponds to a few hundred of meters deep in sedimentary rocks) and for a fluid with a chemical composition in equilibrium with the solid, the dissolution speed at the solid-fluid interface in a limestone is of the order of: v 0 n ≃ 10 -6 to 10 -5 m.year -1

Consequences for the stability of the dissolution process

From the local mechanical equilibrium and the nature of the chemico-mechanical coupling, some important considerations can be inferred about the morphological stability of the dissolution surfaces. This behavior depends on the orientation of the surfaces with respect to the far-field stress.

Previously we have shown how to express the force perturbation arising from the mismatch between the solid-fluid interface orientation and the principal axis of the far-field stress tensor x (equation 8).

This relationship holds independently of the relative magnitudes of the principal stresses σxx and σzz. If the largest principal stress is tangential to the interface, which is not the case for stylolites, σ 0 s < 0 and the sign of δf (x) • x is opposite to the slope of the interface ∂xz. Such tangential force perturbation is concentrated at the points lying ahead of the average dissolution front (Figure 4). The elastic forces concentrate stress at the valleys of the dissolution front where the free energy is thus higher. This leads to an increased dissolution speed at the points lying ahead of the averaged front. The dissolution propagates downwards. The points at the crests, i.e. located behind the averaged dissolution front, show a reduced rate of dissolution thus pushing them further from the average front. The points lying out of the average dissolution plane tend therefore to depart further from the average position. The elastic force is in this situation a destabilizing force. On the contrary, the surface tension tends to stabilize the process by decreasing the surface area by flattening the interface.

The competition between the elastic long-range destabilizing forces and the surface tension short-range stabilizing forces leads to the ATG interface instability. The fastest growing wavelength is determined by the balance between these long-range destabilizing and short-range stabilizing effects. Such instability arising in stressed solids was studied theoretically [START_REF] Asaro | Interface morphology development during stress-corrosion cracking. 1. Via surface diffusion[END_REF][START_REF] Grinfeld | Instability of the interface between a nonhydrostatically stressed elastic body and a melt[END_REF][START_REF] Misbah | Dynamics of a dissolution front for solids under stress[END_REF] and observed experimentally in stressed soluble crystals immersed in a saturated fluid [START_REF] Brok | The effect of elastic strain on the microstructure of free surfaces of stressed minerals in contact with an aqueous solution[END_REF][START_REF] Koehn | Transient dissolution patterns on stressed crystal surfaces[END_REF].

If the largest principal stress lies perpendicular to the interface, as for stylolites, σ 0 s > 0 and the sign of δf (x) • x is the same as the slope of the interface ∂xz. Such tangential force perturbation is concentrated at the points lying behind of the average dissolution front (Figure 5). The elastic forces concentrate stress at the crests of the dissolution front where the free energy is thus higher. This leads to an increased dissolution speed for the points lying behind the averaged front. The dissolution propagates downwards. The points at the valleys, i.e. located ahead of the averaged dissolution front, tend to come back to the average position. The elastic force is a stabilizing force in this situation. Here, the surface tension is again a stabilizing process.

Since the long-range elastic force and the short-range surface tension force are stabilizing forces, if the modeled solid properties are purely homogeneous (i.e. homogeneous elastic solid with homogeneous dissolution rate properties), the model predicts the flattening of any initial non-plane surface with time.

Consequently, to model the morphogenesis of stylolites, which are rough surfaces, we will take here into account the disorder linked to the material properties.

Consequence on initial evolution of trapped fluid pocket

In summary, the above arguments show that an elementary bump of a flat surface disappears for σs > 0, or grows for σs < 0. Qualitatively, if the argument on the stability of surfaces depending on their orientation on the principal stress axis extends for more local orientations along trapped fluid pockets, one should observe the following: for the sides of a fluid pocket lying tangentially to the largest stress, these should develop instable grooves penetrating into the solid, similarly to the ATG instability case. On the contrary, the sides normal to the largest stress direction should remain relatively flat, apart from the fluctuations due to the disorder. These small variations along the surfaces normal to the principal stress axis, and the penetrations of grooves of characteristic wavelength in the rock along the direction of the weakest stress, should lead to the development of elongated structures, and merge initially separated fluid pockets (or clay-enriched pockets). This qualitative mechanism is illustrated on Figure 6. This expectation of qualitative evolution is indeed compatible with the mechanism of development of anti-cracks numerically obtained by [START_REF] Koehn | Instabilities in stress corrosion and the transition to brittle failure[END_REF]. The experimental grooves observed along the free surface on the sides of a fluid-filled cylindrical pocket by Den [START_REF] Brok | The effect of elastic strain on the microstructure of free surfaces of stressed minerals in contact with an aqueous solution[END_REF] also displayed this trend.

Expression of the dissolution speed perturbation as a function of the interface shape

To model the disorder in the solid we assume that the material properties (related to the solid grains) vary in a random and spatially uncorrelated way. This disorder can originate from the diversity of grain composition, grain size or orientation, i.e. it represents the small scale heterogeneities present in the rock For example, the dissolution rate k can be expressed as an averaged term k0 plus some spatial variations of zero average η(x, z) • k0:

k = k0(1 + η(x, z(x))) (20)
The random variable η is a quenched disorder with no spatial correlations and is characterized by its mean < η >= 0 and its variance < η 2 > assumed to be small enough to keep small local slopes. The dynamics of the dissolving interface z(x, t) can be expressed from equations (9-17) as:

v = -∂tz = kΩ 2 RT (∆ue + γ∂xxz), = k0Ω 2 RT (1 + η){(1 + ν)[(σ 0 ij + σ 1 ij )(σ 0 ij + σ 1 ij ) -ν(σ 0 kk + σ 1 kk ) = k0Ω 2 RT {(1 + ν)[σ 0 ij σ 0 ij -ν(σ 0 kk ) 2 ]/4E -u e ref } + k0Ω 2 RT η{(1 + ν)[σ 0 ij σ 0 ij -ν(σ 0 kk ) 2 ]/4E -u e ref } + k0Ω 2 RT {(1 + ν)[2σ 0 ij σ 1 ij -2ν(σ 0 kk σ 1 kk )]/4E + γ∂xxz}
i.e., using Eq.( 18) for the expression of σ 0 ij σ 0 ij -ν(σ 0 kk ) 2 , a dissolution speed separated between an average homogeneous speed v0 and a leading order of the perturbations v1, Figure 6. Expected stability or instability of the dissolution front around a trapped fluid pocket.

X -7

first order in ǫ as ∂tz(x, t) = -v0 -v1(x, t)

with v0 the dissolution speed given by equation ( 18) and v1 a deviation of the dissolution speed with respect to the average dissolution speed v0 expressed as:

v 1 = k0Ω 2 RT (αp 2 0 -α ref p 2 ref ) E η(x, z(x)) + k0Ω 2 RT γ∂xxz(x) + k0Ω 2 RT [(1 + ν)σ 0 ij σ 1 ij -νσ 0 kk σ 1 ll ] 2E ( 23
)
σ 1 is the stress perturbation mentioned previously in equation (2). It is generated by the surface distribution of the tangential force perturbation δf (x) due to the irregular nature of the interface.

The first term is a quenched disorder term leading to the roughening of the interface. The second one is a stabilizing quadratic short-range term arising from the surface tension. The last term can be expressed via a non-local kernel from the shape of the interface z(x) by integrating the elastostatic equations in the solid half-plane.

Detailed form of the elastic long-range interaction kernel

The stress perturbation induced by the force perturbation δf (x) (equation 8) exerted on the surface can be determined via the Green function method. Following [START_REF] Landau | Theory of elasticity[END_REF], the displacement induced by an elementary force x applied at the origin (0, 0, 0) on a semi-infinite solid is:

ax (x, y, z) = 1 + ν 2πE 2(1 -ν)r + z r(r + z) + (2r(νr + z) + z 2 r 3 (r + z) 2 x 2 ay (x, y, z) = 1 + ν 2πE 2r(νr + z) + z 2 r 3 (r + z) 2 xy az (x, y, z) = 1 + ν 2πE (1 -2ν)x r(r + z) + zx r 3 ( 24
)
where r is the distance relatively to the force application point at (0, 0, 0), i.e. r 2 = x 2 + y 2 + z 2 . The associated strain applied on the solid is:

ǫ e ij = 1 2 (∂iaj + ∂jai) (25)
and the associated stress is:

fij (x, y, z) = E 1 + ν (ǫ e ij + ν 1 -2ν ǫ e kk δij ) (26)
The stress associated to the point force x applied on the surface of normal ẑ at the origin is equal at the origin itself to xẑ + ẑ x. Since the model treated here is invariant by translation along y, the force perturbation δf (u) = σ 0 s (∂uz)(u)x is exerted at any v ∈] -∞, ∞[ and the resulting displacement field at (x, y, z), is solely dependent on (x, z) and can be expressed, by linearity of the elastostatics equations, (similarly to the elastostatic Green function method detailed in Eq.(8.14) by [START_REF] Landau | Theory of elasticity[END_REF]), as a displacement field w of components

wi(x, y = 0, z) = ∞ u=-∞ ∞ v=-∞ ai(x -u, -v, z)dudv δf (u).x (27)
The associated strain perturbation is

ǫ p ij = 1 2 (∂iwj + ∂jwi), (28)
and the associated stress,

σ 1 ij (x) = E 1 + ν (ǫ p ij + ν 1 -2ν ǫ p kk δij ) + δf (x)(δixδjz + δizδjx)δ(z), (29)
where the first term represents the stress induced by the elastic deformation, and the second one the direct application of the force perturbation on the surface. In the above, the spatial derivative of Eq.( 28) can be exchanged with the integration in Eq.( 27), to obtain

ǫ p ij (x, y = 0, z) = ∞ u=-∞ ∞ v=-∞ ǫ e ij (x -u, -v, z)dudv δf (u).x. ( 30
)
Recalling the expression of the force perturbation, Eq.( 8), from Eq.( 29), the stress perturbation along the surface, at z = 0, is thus:

σ 1 ij (x) = σ 0 s • p.p.[ ∞ u=-∞ du (∂uz)(u) * ∞ v=-∞ fij (x -u, -v, 0)dv] +σ 0 s (∂xz)(x)(δixδjz + δizδjx) (31)
where p.p. refers to the principal part of the integral. Taking the derivatives of the displacement field (equation 25), we can calculate the associated stress. Integrating this result along the y-axis gives:

∞ v=-∞ fij (x, -v, 0)dv = - 2ν πx (δixδjx + δiyδjy) (32)
and thus,

σ 1 (x) = - 2νσ 0 s π • p.p.[ ∞ x ′ =-∞ du (∂uz)(u) x -u ](xx + ŷ ŷ) +σ 0 s • (∂xz)(x)(xẑ + ẑx) (33)
The elastic energy perturbation associated to the interface deformation can be computed using equation ( 33) and the relation:

σ 0 = -(p0 -σ 0 s /3)(xx + ŷ ŷ) -(p0 + 2σ 0 s /3)ẑ ẑ (34)
It results in:

u 1 e = [(1 + ν)σ 0 ij σ 1 ij -νσ 0 kk σ 1 ll ] 2E = 2ν[(1 -2ν)p0] πE * σ 0 s • p.p.[ ∞ u=-∞ du (∂uz)(u) x -u ] (35)

Dynamic equation for the dissolution interface

The equation ( 22) rules the dynamics of the interface dissolution. When computed with equation ( 18), it gives:

RT k0Ω 2 v 1 = (αp 2 0 -α ref p 2 ref ) E η(x, z(x)) -γ∂xxz(x) +β p0σ 0 s E • p.p.[ ∞ u=-∞ du (∂uz)(u) x -u ] ( 36
)
where β is a geometrical factor:

β = [2ν(1 -2ν)]/π (37)
Equation ( 23) can be expressed in a dimensionless form by using length and time units as:

L * = γE/(βp0σs) (38) τ = (L * ) 2 RT /(γk0Ω 2 ) (39)
We define the dimensionless variables in the reference frame moving at the average velocity -v0 as:

z ′ = [z + (v0t)]/L * (40) x ′ = x/L * (41) t ′ = t/τ (42)
and the reduced quenched noise as:

η ′ (x ′ , z ′ (x, t) -v0t/L * ) = [(αp 2 0 -α ref p 2 ref )/(βp0σs)]η(x, z(x, t)) (43)
The dimensionless stochastic equation for the stylolite growth process is then:

∂ t ′ z ′ (x ′ , t ′ ) = η ′ (x ′ , z ′ (x ′ , t ′ ) -v0τ t ′ /L * ) + ∂ x ′ x ′ z ′ -p.p.[ ∞ u=-∞ du (∂uz ′ )(u) x ′ -u ] (44)
At large average dissolution speed, the term v0τ t ′ /L * takes over z quickly and the noise is annealed, becoming mostly time-dependent. On the contrary, for sufficiently slow processes such as the extend of the surface roughness over several grains, the noise can be considered as quenched. This is the case here as the changes in η ′ arising from z(x, t) are significantly larger than the changes due to some variations of the average dissolution front position v0τ t ′ /L * . To the first order, the noise dependence is mainly η ′ (x ′ , z ′ (x ′ , t ′ )) and the noise will therefore be considered here as quenched.

The dynamic equation then becomes:

∂ t ′ z ′ (x ′ , t ′ ) = η ′ (x ′ , z ′ (x ′ , t ′ )) + ∂ x ′ x ′ z ′ -p.p.[ ∞ u=-∞ du (∂uz ′ )(u) x ′ -u ] (45)
Alternatively, in some arbitrary spatial unit ℓ, this can also be written:

∂tz(x, t) = η ′′ (x ′ , z ′ (x ′ , t ′ )) + ∂xxz - ℓ L * dy ∂yz x -y ( 46
)
with L * = γE/(βp0σs) and τ = ℓ 2 RT /(γk0Ω 2 ), the time unit.

Small and large scale behavior of the model

Elastic interactions can be neglected in equation ( 46) for small scales such as ℓ ≪ L * (the lower limit corresponds to the resolution of the analyzed signal) reducing the model to a Laplacian description:

∂tz ′ (x, t) = ∂xxz ′ + η(x, z ′ (x)) ( 47
)
This equation is known as the Edwards Wilkinson model [START_REF] Edwards | The surface statistics of a granular aggregate[END_REF] modified with a quenched random noise. It has been studied in the literature and leads to the growth of self-affine surfaces of roughness ζ ∼ 1.2 [START_REF] Roux | Interface roughening and pinning[END_REF], in agreement with existing data on stylolites where ζ ∼ 1.1 [START_REF] Schmittbuhl | Roughness of stylolites: implications of 3D high resolution topography measurements[END_REF].

Conversely, for large scales ℓ ≫ L * (the upper limit corresponds to the system size), surface tension can be neglected reducing equation ( 46) to a mechanical regime:

∂tz ′ (x, t) = - ℓ L * dy ∂yz x -y + η(x, z ′ (x)) (48)
In this case, the model is similar to known models describing the propagation of an elastic line on a disordered pinning landscape or the propagation of a mode I fracture front in a disordered solid. It leads to the growth of self-affine surfaces of roughness ζ ≃ 0.5 [START_REF] Tanguy | From individual to collective pinning: Effect of long-range elastic interactions[END_REF]. In summary, the model derived above predicts the growth of dissolution surfaces with different self-affine characteristics at small scale (ζ1 ∼ 1.2) and large scale (ζ2 ∼ 0.5). The transition between these regimes is expected to occur at a certain cross-over length L * .

Numerical approach : Dynamic evolution of the interface

From a purely analytical point of view and via the similarity of asymptotic form of the dynamic equation with known models for large and small scales, we have shown that two different scaling laws are expected for small and large scales, and that the cross-over length should depend on the far-field stress magnitude. Independently from this general analytical analysis, we will now show how to solve the problem numerically, i.e. implement the dynamic evolution of the interface with all the large and small scale terms and random variables to represent the disorder and analyse the resulting morphogenesis.

Practical implementation of the model

We simulate the dissolution process for a calcite-water interface. This is done in an event-driven discrete lattice code, with algorithms corresponding to a discrete Langevin equation leading to grains getting dissolved one at a time: for each grain along the interface, a time to dissolution is computed from the above Langevin equation, and the grain with the shortest dissolution time is removed. After what, the times are recomputed for all grains along the interface, and the next grain with shortest dissolution time is removed, and so on (see [START_REF] Renard | Three-dimensional roughness of stylolites in limestones[END_REF] for details of the practical implementation). The selected constants correspond to a calcite-water system, γ = 0.27J • m -2 , Ω = 4 • 10 -5 m 3 • mol -1 , ν = 0.25, E = 80 GPa and k0 = 10 -4 mol • m -2 • s -1 [START_REF] Renard | Three-dimensional roughness of stylolites in limestones[END_REF]. The chosen physical conditions are T = 420 K, < p >= 10 MPa and < σs >= 40 MPa. The amount of quenched noise is associated to the natural variations of grain properties. The typical scale associated to the quenched disorder (or typical grain size) is considered here to be around ℓ = 10µm, with no correlation above this scale. This quenched disorder has a standard deviation < η 2 > = [αℓp0/(βL * σs)] • [(δE/E) + (δk/k0) + (δα/α)] corresponding to some relative variations of the dissolution rate of around 10% (i.e. δk/k0 ∼ 0.1).

The dimensionless surface dynamic equation without disorder is:

∂tz(x, t) = v0 + ∂xxz - ℓ L * dy ∂yz x -y (49)
where L * = γE/(βp0σs), ℓ is the unit length, and τ = ℓ 2 RT /(γkΩ 2 ) is the time unit. We assume a small disorder in the implied quantities (e.g. Young's modulus), that are quenched in the material properties of the rock heterogeneity associated with micrometric grains, typically ℓ = 10µm. The interface is supposed to be normal to the largest stress direction (stabilizing elastic interactions).

Considering a perturbation to the first order, in the referential frame of the homogeneously moving average front, z ′ = z -v0t, the equation ruling the surface growth becomes:

∂tz ′ (x, t) = ∂xxz - ℓ L * dy ∂yz x -y + η(x, z(x)) (50) with a quenched random term η(x, z ′ (x)) = [αℓp0/(βL * σs)]• [(δE/E) + (δk/k) -(δα/α)]
The first and second terms are stabilizing terms. The third term refering to the quenched disorder destabilizes the interface. We perform the simulation of this dynamic equation with both stabilizing terms and quenched noise.

The prefactors in equation 50 depend on the rock type and on the applied stress. In addition to these mappings, the characteristic units are known as function of the rock properties. The cross-over scale L * = γE/(βp0σs) is function of the pressure during the growth, through p0 and σs.

Determining the cross-over length L * for natural samples allows to determine such stress value during the growth, and consequently the depth of the rock during the stylolite propagation. Assuming as an order of magnitude p0 ∼ σs and typical values for the limestone elastic properties and the water calcite reaction rates, L * ∼ 1mm leads to a typical depth of 1 km. Stylolites can thus be considered as fossils of the stress magnitude.

We solved the dynamic equation ( 46) with an eventdriven algorithm where the fastest dissolving grain is removed at each step.

The problem is considered as L-periodic and the long-range elastic kernel p.p. dy ∂y z

x-y = -p.p. dy z(x)-z(y) (x-y) 2 is replaced by its finite-size form -p.p.

L 0 dy z(x)-z(y) sin 2 (π(x-y)/L) π 2
L 2 . This standard form can be obtained by solving the elastostatic equations in the Fourier space and performing an inverse Fourier integration. When a new grain is reached, the realization of its quenched disorder η is evaluated using a Gaussian distribution. For the dissolution surface simulated which is 4096ℓ long, 8000000 grains were dissolved.

Analysis of the small-scale and large-scale roughness of the saturated interface

The simulation of the calcite-water system leads to the growth of a dissolution interface. Starting from a flat interface and after a certain transient time, the Fourier modes saturate to a characteristic amplitude. A snapshot of the developed stylolitic interface is shown in Figure 7.

This interface fluctuates around the average progressing flat dissolution front. The Fourier power spectrum P (k, t) = z(k, t) 2 of each front z(x, t) is extracted, and the ensemble average of this power spectrum P (k) = z(k) 2 is obtained for developed interfaces, by averaging over all fronts after 80000 grains have been dissolved. The expected small and large scale self-affine characteristics correspond to the theoretical predictions, as shown in Figure 8. Indeed, the power-spectrum is a power-law of scale, with two different exponents at large and small scale, and a cross-over length around the scale L * : For k > 2π/L * i.e at small scale, we have P (k) ∼ k -1-2ζ with ζS = 1.2, and for the large scales, the roughness exponent is found to be around ζL = 0.35. The straight lines in the bilogarithmic axes (Figure 8) correspond to these power law behaviors, determined by linear regression over the two domains k > 2π/L * and k < 2π/L * . The ensemble used for the roughness estimate is the following: it corresponds roughly to 100000 grain being dissolved after the first 80000 first ones, which are discarded. We thus compute the average power spectral density profiles over all these states, representative of a saturated situation with fluctations of the Fourier mode amplitude around some characteristic magnitude for each wavelength. The linear regression have been performed in bilogarithmic space on the ranges 0 < log10(k) < 1.5 and 2 < log10(k) < 3, with k-unit of 2π/L, with L = 4096ℓ and a grain size ℓ = 10µm. The standard error bar provided by the linear regression over this two ranges is around ±0.2 in slope (i.e. ±0.1 for the Hurst exponents ζ).

Thus, we find that the scaling of saturated surfaces in this model is compatible with observations made on natural surfaces, and with the previous analytical predictions.

In addition, the dynamic behavior of these models (Edwards Wilkinson in a quenched noise [START_REF] Roux | Interface roughening and pinning[END_REF], or elastic string in a disordered landscape [START_REF] Tanguy | From individual to collective pinning: Effect of long-range elastic interactions[END_REF]) is known. The prefactor (characteristic time) associated with the dynamics can be evaluated through the previous computations from the rock material properties. The time to saturation at an observation scale of a few centimeters is estimated to be around a few thousands of years. The stylolite roughness is hence always in a saturation state for a geologist at small observation scale.

However, for longer systems, e.g. decametric ones, much longer times would be required for saturation. Such long stylolites are sometimes observed but rarely analyzed in terms of scaling of the height. To our knowledge, the only analysis performed on decametric size stylolites [START_REF] Ben-Itzhak | Pressure solution and stylolites in carbonate rocks[END_REF]Laronne Ben-Itzhak et al., 2012b] showed that these large scale structures were not saturated. This means that the time during which the stylolitization was active on such very long stylolites was only enough to lead the small scales to saturated amplitude, but not the large ones (above a few tenth of centimeters). Figure 10. Functions obtained from profiles 1, 2, 3 and 4. A grey-level threshold was imposed on the pictures to isolate the stylolites. The functions were then obtained by selecting the mean limit of the pixels.

Example: Application of the model to natural data

The model is applied to a sedimentary stylolite collected in a core at the Andra (French national radioactive waste management agency) Underground Research Laboratory (URL) at Bure in Eastern France. The selected sample comes from the borehole EST433 at a depth of 720 m. The host rock is a fine-grained, homogeneous grainstone from the Dogger age. The core was cut in three parts thus giving four profiles for analysis (Figure 9).

Profiles 1 and 2 and profiles 3 and 4 are spaced by 3 mm (thickness of the drilling saw) and profiles 2 and 3 are spaced by 30 mm. Each profile has a length around 90 mm. The stylolites were photographed at a resolution of 30 µm. A systematic method was used to extract profiles from the photographs. It consists on isolating the black pixels constituting the clay particles in the stylolite from photographs converted in grey level pictures. The profiles will be used as functions in the spectral analysis (integral transforms) and thus are required to be single-valued. Stylolites show a self-affinity geometry [START_REF] Schmittbuhl | Scaling invariance of crack surfaces[END_REF][START_REF] Barabási | Fractal concepts in surface growth[END_REF] meaning that they are statistically invariant under an affine transformation. Thus, for ∆x and ∆y the horizontal direction amplitude and ∆z the vertical direction amplitude: ∆x → λ∆x, ∆y → λ∆y and ∆z → λ ζ ∆z, where λ can take any value and ζ is the Hurst exponent which describes the scaling invariance [START_REF] Schmittbuhl | Roughness of stylolites: implications of 3D high resolution topography measurements[END_REF][START_REF] Renard | Three-dimensional roughness of stylolites in limestones[END_REF]. As in Ebner et al. [2009b] we used both the Fourier power spectrum [START_REF] Schmittbuhl | Scaling invariance of crack surfaces[END_REF] and the averaged wavelet coefficient [START_REF] Simonsen | Determination of the Hurst exponent by use of wavelet transforms[END_REF] signal processing methods to analyse the profiles (Figure 10). We used two different methods to check the repeatibility of the results. First we calculated the Fourier power spectrum P (k), which is the square of the modulus of the Fourier transform, as a function of the wave-number k (k = 2π/L, where L is the wavelength). The power spectrum expressed as a function of the length for a selfaffine profile behaves as P (L) ≃ L 2ζ+1 . We calculated also the averaged wavelet coefficient spectrum as a function of the scale a with Daubechies 4 wavelets which behaves as W (a) ≃ a 1/2+ζ . The results show the two scaling regimes predicted by the theory presented above, described by two different power laws. Figure 11 shows the Fourier power spectrum for the profile 1 as a function of the length L. The raw data are more concentrated at small scale. The lower limit for the length corresponds to the Nyquist length which is the resolution multiplied by 2. As the profiles have a finite-size the upper limit for the analysis (corresponding to small wavenumber) is given by the size of the profile. To analyze the data, we apply a logarithmic binning so that the weight on each point is equal. To estimate the cross-over length, we used a linear-by-part fit with a cross-over function changing the scaling law from small to large scale as explained in Ebner et al. [2009b]. The averaged wavelet coefficient spectrum (Figure 12) does not require a binning. The same kind of fitting was used to appraise the cross-over length L * .

The intersection between both regimes (whose slopes are imposed by ζS = 1 and ζL = 0.5 for small and large scale respectively) gives the cross-over length L * . We summarize the estimated cross-over length for all the analyzed profiles with both methods in Table (1). The uncertainties on the cross-over length (68% and 44% for Fourier power spectrum and averaged wavelet coefficient respectively) are due to the spatial variability of the intersection between the small and large regimes.

Discussion

Interpretations of the estimated paleostress

We discuss three theories about the meaning of the estimated paleostress for the studied sedimentary stylolite:

• Present day stress: if the conditions for pressuresolution (lithostatic pressure in competition with the presence of a fluid at an appropriate state of equilibrium) are present, a stylolite should show the present day state of morphology and is still evolving. This means that we should measure the current applied stress and see the last evolution of the morphology. This can be compared with recent studies where vertical and horizontal stresses were measured in boreholes at Bure [START_REF] Wileveau | Complete in situ stress determination in an argillite sedimentary formation[END_REF][START_REF] Gunzburger | Rheological characterization of a sedimentary formation from a stress profile inversion[END_REF] to assess if the estimated stress corresponds to the measured ones.

• Evolution stopped: this can occur if the lithostatic stress becomes too small to encourage the process (change in the magnitude due to a tectonic phase for example). It can also be associated with the closing of the porosity by recrystallization. Indeed, if the pore size decreases because of recrystallization at the pore surface, the surface tension increases preventing more recrystallization. Thus, the water is getting more charged in dissolved materials and the chemistry of the water changes and can stop the stylolite evolution. Moreover, the decrease of the pore size can limit or stop the fluid flow and close the system.

• Reactivation: both previous theories can act on the history of a stylolite. After its initiation, a stylolite can see its growth stopped by the kind of process we developped just before. If in the geological history an event as emerged soil and/or erosion allows to change the applied stress or to meteoritic fluids to flow in the soil, the system can have its properties changed and pressure-solution process can start again until it is stopped or it can still evolve.

These three theories will be discussed with regard to the paleostress results.

Estimation of the paleostresses 5.2.1. Hypothesis on the basin evolution

To estimate the paleostress from the model developped in this study, we use the average of the cross-over lengths determined for the four profiles. The cross-over length is related to the stresses by equation ( 38). However, this equation can be simplified by making assumptions on the surrounding rock formation. We use the same assumptions as in Ebner et al. [2009b] for the stylolites from Cirque de Navacelle (Cévennes, France) considering that the initiation of stylolites occurs at the early stage of a basin. The major principal stress is vertical (σzz) as we analyzed a sedimentary stylolite. The principal horizontal stresses are isotropic (σxx = σyy). Thus the mean stress p0 and the shear stress σS are:

p0 = -(2σxx + σzz)/3 (51) σS = σxx -σzz (52)
As stylolites are known to develop in the early stage of sedimentation of basins, the strain is assumed to be uniaxial:

σxx = σyy = ν 1 -ν σzz (53)
Using equations (51-53), equation (38) becomes:

σ 2 zz = γE αβL * (54)
where

α = 1 3 (1 + ν) (1 -ν) (1 -2ν) (1 -ν) (55)
is a dimensionless geometrical factor. The geometrical factor β (equation 37) is β = ν(1 -2ν)/π. Using the average cross-over length L * in equation ( 54), we can estimate the main principal paleostress σzz. The Poisson's ratio ν of the host rock was determined by measuring the P and S elastic Figure 12. Averaged wavelet coefficient spectrum of the profile 1. A linear-by-part fitting were run on the data [Ebner et al., 2009b]. Two different scaling regimes are observed at small and large scale with Hurst exponent around 1 and 0.5 respectively. The fit reveals a cross-over length L * around 1.95 mm. This is in good aggreement with the length inferred using the Fourier power spectrum method.

wave velocities (ν = 0.5(V P /V S ) 2 -1 ((V P /V S ) 2 -1) ). The relative errors for the measurements of VP and VS are 1 and 2% respectively [START_REF] Benson | Pore fabric shape anisotropy in porous sandstones and its relation to elastic wave velocity and permeability anisotropy under hydrostatic pressure[END_REF]. The relative error for the Poisson's ratio is thus equal to 12%. Consequently, α and β have error bars equal to 2% and 22%, respectively. The last constant to be determined is the Young's modulus. The next paragraph details our choices for this matter.

Uncertainties on the Young's modulus E

In their paper, Ebner et al. [2009b] determined E assuming the vertical stress is equal to the lithostatic stress as in equation ( 56) where z is the current depth of their samples. They plotted the determined stress as a function of L -1/2 . The slope of the curve is proportional to E 1/2 (see equation ( 54)). They found E = 15 GPa which is the lowest acceptable limit for limestones [START_REF] Clark | Handbook of physical constants[END_REF]. Based on uniaxial loading made in our laboratory, we determined E = 36.2 GPa for the rock surrounding the analyzed stylolite. Considering that the limestones from Bure replaced in the geological context of the Paris basin cannot be excessively harder than what we observe today, the value determined in the laboratory is taken as the upper limit for E. Thus we can estimate the paleostress in a small range of E. The values used for the calculation of the paleostress are summarized in Table (2).

To calculate the paleostress σzz, we take into account the error bars for each parameter. The computed error for the calculation of σzz is 66% for the Fourier power spectrum method and 54% for the averaged wavelet coefficient method. The results are summarized in Table (3).

Geological context

By doing some assumptions on the sedimentary overburden, the depth of development of the stylolite can be assessed. The lithostatic pressure σzz can be expressed as:

σzz = ρgh ( 56
)
where ρ is the density in g • m -3 , g is the Earth's gravity (g = 9.81 m • s -2 ) and h is the depth in m. We make the assumption that at the initiation of the stylolite, linked to the early stage of formation of the sedimentary basin, the overburden was made of limestones only. Thus, we consider the density of limestones ρ = 2710 g • m -3 . The estimated depths of development of the stylolites are summarized in Figure ( 4). The error bars on h are of the same order as for the paleostress. Now we can wonder what is the interpretation of the estimated paleostress with regard to the three theories exposed previously:

1. [START_REF] Wileveau | Complete in situ stress determination in an argillite sedimentary formation[END_REF] and [START_REF] Gunzburger | Rheological characterization of a sedimentary formation from a stress profile inversion[END_REF] measured the vertical stress at Bure which is equivalent to the lithostatic pressure as in equation ( 56). Our results show that the calculated depth corresponds to the depth where we cored the analyzed stylolite. Thus the studied stylolite is more likely to be still active and to show the present day stress.

2. [START_REF] André | Chronology of fracture sealing under a meteoritic fluid environment: microtectonic and isotopic evidence of major Cainozoic events in the eastern Paris Basin (France)[END_REF] discussed about a reactivation of the stylolitization during the Tertiary age (end of Cretaceous more precisely) by the change in the stress orientation or by the emergence of the Cretaceous sediments which Table 3. Results for the calculation of the paleostress σzz using the averaged cross-over length for the Fourier power spectrum (FPS) and average wavelet coefficient (AWC) methods. We calculated the paleostress taking into account the variability of the Young's modulus E where E low is the lower limit for the Young's modulus for limestones and Eup is the determined Young's modulus for the studied sample. were eroded and permitted to meteoritic fluid to spread in the sediments. This reactivation process could have acted on the growth of the studied stylolite until today. But still it seems that the theory of the present day stress is more applicable on that example.

Conclusions

Analysing the local boundary conditions due to the fact that the inside of a stylolite does not sustain shear stress and an elastic surrounding, we derived the dependence of the free energy along a stylolite surface on the shape of the stylolite. Adding up a surface energy term we derived a dynamic surface evolution model for a stylolitic interface. This model, in the situation where a stylolite is perpendicular to the largest principal stress axis -as in most case -includes terms that lead to the stabilization of the surface dynamics, i.e. to the vanishing of initial perturbations towards a flattening surface. Hence, the presence of disorder linked to the heterogeneities of the material properties is required to explain the rough nature of stylolites. Introducing such non correlated quenched disorder, the model predicts the occurence of two scaling laws. At small scale, a destabilizing disorder competing with a stabilizing surface energy term give a model similar to the Edwards Wilkinson model in a quenched noise leasing to a saturated surface with a Hurst exponent around 1. At large scale, the competition between destabilizing disorder and stabilizing elastic interactions is similar to models of evolution of an elastic interface in quenched disorder leading to a Hurst exponent of 0.5.

The cross-over scale between these two scaling regimes was shown to be directly linked to the stress magnitude. Hence, the determination of this cross-over and other physical rock properties allows to use stylolites as markers of the paleostress magnitude.

Both scaling laws and the dependence of this cross-over scale on the stress magnitude were derived in two ways: by purely analytical derivation and similarity to known models in section 2 and by numerical integration in section 3.

Importantly, it should be noted that the elastic forces, depending on the boundary conditions, can be stabilizing, as here, or destabilizing. The existence of several models and techniques of global or local calculation of the free energy can raise the question of a stabilizing or destabilizing nature of the elastic forces in the context of a stylolite. Independently from the derivation carried out in details in this paper, we note the following argument that can distinguish between stabilizing and destabilizing terms. The only difference between models with stabilizing or destabilizing elastic kernel is the sign of the prefactor in front of the elastic operator in the dynamic equation. However, when this sign is reverted, all large scale wavelength Fourier modes become unstable (with a selection of fastest growing mode, as e.g. shown in [START_REF] Misbah | Dynamics of a dissolution front for solids under stress[END_REF] or [START_REF] Bonnetier | Does roughening of rock-fluid-rock interfaces emerge from a stress-induced instability?[END_REF]). Numerical simulations similar to the ones shown above, with a destabilizing mode, do not lead to any saturation of the amplitude of the large modes at long times, and the Fourier power spectrum at a given time does not display any scaling law at fixed time for the large scales. Thus, the scaling laws observed in field stylolites are compatible with a model where elastic forces are stabilizing: we take this as a good sign of validity of the proposed approximations to take the boundary conditions into account in the proposed model.

The results from both analytical and numerical independent resolutions presented in this study are also consistent with three other independent observations:

• The existence of two Hurst exponents at small and large scales, as observed in [START_REF] Schmittbuhl | Roughness of stylolites: implications of 3D high resolution topography measurements[END_REF], in the stylolites from the log cores of Bure (section 4).

• The results of recent molecular dynamic models of dissolution with pressure reliance and surface energy terms in the free energy displaying similar scaling laws and an identical law for the dependence of the cross-over length over the applied stress [START_REF] Koehn | Modelling of stylolite geometries and stress scaling[END_REF].

• The model was applied in a previous study to stylolites found at various depths in a limestone formation at Cirques de Navacelles (Cévennes, France). The inferred formation stresses were compatible with the derived weight of overburden at the time of formation [Ebner et al., 2009b].

We show finally on the example of sedimentary stylolites in Bure, how the confinement stress can be derived from morphological studies of stylolites. The ubiquitous nature of these pressure-solution features makes them a versatile marker for paleostress magnitude that can give access to the stress during the growth of stylolites. This easily available paleostress marker opens the way for systematic studies of paleostress in large rock formations for different stylolite families. However, it must be used carefully as the error bars are not minor. An important number of measurement is required to constrain the results. Together with dating indications for the time of occurence of such stylolites (as e.g. times of tectonic events) and current stress assessment methods it opens the way for the determination of stress evolution in large basins, which is a key to understand their evolution.

Figure 2 .

 2 Figure 1. Initial stage of a stylolite: Trapped elongated fluid pocket.

Figure 3 .

 3 Figure

Figure 4 .

 4 Figure 4. Surface tangential to the largest stress (σxx) axis: unstable case, Azaro-Tiller-Grinsfeld instability.

Figure 5 .

 5 Figure 5. Surface normal to the largest stress (σzz) axis: stable case.

Figure 7 .

 7 Figure 7. Snapshot of the pressure-solution profile.

Figure 8 .

 8 Figure 8. Average power spectrum of simulated stylolitic fronts, in bilogarithmic representation. The k-unit is 2π/L, with L = 4096ℓ and a grain size ℓ = 10µm. The vertical unit is arbitrary. The crossover is obtained at 2π/L * .

Figure 9 .

 9 Figure 9. Profiles 1, 2, 3 and 4 from right to left. A core from the Dogger formation (EST433 well) was cut in three parts to obtain four profiles. Each profile was photographed at high resolution. The picture at the bottom shows the profile number 2.

Figure 11 .

 11 Figure11. Fourier power spectrum of the profile 1. The raw data were binned logarithmically to run a linear-bypart fitting on the data[Ebner et al., 2009b]. Two different scaling regimes are observed at small and large scale with Hurst exponent around 1 and 0,5 respectively. The fit reveals a cross-over length L * around 1.14 mm.

Table 1 .

 1 Summary of the cross-over length found for the four profiles analyzed by Fourier power spectrum (FPS) and average wavelet coefficient (AWC).

	Profile number	1	2	3	4	Average, L *
	L * F P S (mm) L * AW C (mm) 1.95 1.52 0.72 1.60 1.14 0.37 0.37 1.13	0.75±0.51 1.45±0.64

Table 2 .

 2 Summary of the estimated paleostress for the stylolites from the cross-over length.

	γ(J • m -2 ) Eup(GP a)	ν	α	β
	0.27	36.2±0.4 0.37±0.04 0.32±0.01 0.033±0.007

Table 4 .

 4 Results for the calculation of the depth for the Fourier power spectrum (FPS) and average wavelet coefficient (AWC) methods. We consider an early stage of formation of a sedimentary basin with an overburden made of limestones only.

	Paleostress (MPa) L * , E low	L * , Eup
	σ F P S	22.6±14.9 35.1±23.2
	σ AW C	16.3±8.8 25.3±13.7
	Depth (m)	L * , E low	L * , Eup
	h F P S	850.1±561.1 1320.3±871.4
	h AW C	613.1±331.1 951.7±513.9

Acknowledgments

We would like to thank C. Aurière, J.D. Bernard, T. Reuschlé, A. Steyer, E. Aharonov, Y. Bernabé and J.L. Alves for technical help and discussions. This work was partially funded by the French national radioactive waste management agency (Andra) and supported by a FORPRO grant. The authors would also like to thank the CAPES-COFECUB program for stimulating fruitful discussions.