Abdallah Said
email: abdallah.said@hds.utc.fr

Reine Talj
email: reine.talj@hds.utc.fr

Clovis Francis

Hassan Shraim
email: hassan.shraim@ul.edu.lb

Tuning and Costs Analysis for a Trajectory Planning Algorithm for Autonomous Vehicles

Keywords: Autonomous vehicle, trajectory planning, cost analysis

  L'archive ouverte pluridisciplinaire

INTRODUCTION

An autonomous vehicle needs a local trajectory based on real-time vehicle status and dynamic environment perception data (e.g., nearby cars, road conditions) to guarantee safe maneuvers while following the global trajectory. Local trajectory planning is defined as the planning of a vehicle's transition from one possible state to the next one while satisfying the vehicle's kinematic restrictions based on the vehicle dynamics and constrained by the passenger's comfort, lane boundaries, and traffic rules while avoiding static and dynamic obstacles. Different trajectory planning approaches have been developed for the navigation of autonomous vehicles [START_REF] Dixit | Trajectory planning and tracking for autonomous overtaking: State-of-the-art and future prospects[END_REF], [START_REF] Katrakazas | Real-time motion planning methods for autonomous on-road driving: State-of-the-art and future research directions[END_REF]. They differ in how they deal with the environment and the vehicle dynamics limitations. According to the literature, there are four well-known approaches to trajectory planning methods: grids generalization [START_REF] Pivtoraiko | Efficient constrained path planning via search in state lattices[END_REF], samplingbased planning [START_REF] Karaman | Sampling-based algorithms for optimal motion planning[END_REF], numerical optimization methods [START_REF] Lim | Hierarchical trajectory planning of an autonomous car based on the integration of a sampling and an optimization method[END_REF] and curve interpolation methods [START_REF] Alia | Local trajectory planning and tracking of autonomous vehicles, using clothoid tentacles method[END_REF]. The last approach is adopted, and it aims to generate trajectories on a given horizon with a specified geometric shape that responds to one or more constraints such as vehicle dynamics and kinematics, comfort, road shape, and curvature continuity. Each candidate trajectory is assessed using a cost function that takes into account several factors such as traveled distance and execution time costs, acceleration, collision verification, and other performance criteria. Then, the cost function is minimized to find the best trajectory among several candidate trajectories generated by the planning algorithm. Each one of these costs is multiplied by a weighting coefficient to rank its significance. The major weakness of this approach is that, for some trajectory shapes, the result of the planning algorithm is highly dependent on the tuning of the cost weights . Moreover, function of the shape of geometric curves, which represent candidate trajectories, the algorithm is more or less dependent on tuning weights. It is important to evaluate the sensitivity of the planning algorithm to the tuning parameters in order to evaluate its ability to deal with the different driving situations without the need to adapt the parameters. In [START_REF] Zhang | Toward a more complete, flexible, and safer speed planning for autonomous driving via convex optimization[END_REF], the planner's flexibility is defined by how many different types of scenarios it can manage by just altering parameters rather than modifying the fundamental issue formulation or problem structures. In [START_REF] Mouhagir | Trajectory planning for autonomous vehicle in uncertain environment using evidential grid[END_REF], trajectory planning based on clothoid tentacles was adopted. Different combinations for cost function weights were proposed and the results show that the proposed method was very sensitive to parameter variation and some combinations were not suitable for real-time application. A study is presented in [START_REF] Arnay | Safe and reliable path planning for the autonomous vehicle verdino[END_REF] on the impact of the weights of a cost function with five criteria on the final behavior of the vehicle. This planner deals with only pedestrians as obstacles and low velocities. Several configurations were chosen and ranked based on two behaviors.

In order to choose the best combination of weights, the role of each cost must be investigated and its influence in choosing the best candidate trajectory must be understood. In this paper, impact analysis and tuning of weighting coefficients are done on the planning algorithm presented in [START_REF] Said | Local trajectory planning for autonomous vehicle with static and dynamic obstacles avoidance[END_REF]. This paper is organized as follows: in Section 2, we present the trajectory planning algorithm. The global cost function and its components are detailed in Section 2.2. Section 3 presents the proposed method and reports the simulation results with some analysis on the cost weights, while the final Section 4 concludes the paper.

TRAJECTORY PLANNING

Presentation of the trajectory planning method

The trajectory planning algorithm must provide the best trajectory from a set of candidate trajectories that helps the vehicle to track a reference trajectory while avoiding static and mobile obstacles and ensuring safety and passenger comfort. ing the vehicle on the reference trajectory, a set of candidate trajectories are generated. they cover either the host lane or the entire width of the road depending on their navigability. Each of them consists of two phases: The transient phase, which is modeled by a 4 th order polynomial curve to provide a smooth change ensuring continuity of the curvature, starts at the actual position of the vehicle up to a defined lateral offset from the reference trajectory. Then, the candidate trajectory continues with a permanent phase parallel to the reference trajectory where the lateral offset becomes constant (see Fig. 2). Secondly, an obstacle detection procedure is carried out. A classification area is formed along the candidate trajectory by the footprint of the vehicle. The collision distance d obs , which is the free distance traveled on the trajectory to reach the first obstacle, is then detected. Based on the collision distance and security distance, we classify the candidate trajectories into three classes: non, partially, or fully navigable trajectory. Note that the security distance is the distance the vehicle must maintain between it and the encountered obstacle. It depends on the obstacle state (static or dynamic). It is calculated based on the safe stop distance, the needed distance to stop the vehicle from its actual speed, with a defined comfortable deceleration. If there are no navigable trajectories, the algorithm selects the one with the longest distance to the obstacle to stop the vehicle with high deceleration (safe stop scenario). Thirdly, the navigable trajectories are evaluated according to various criteria, including smoothness, safety, consistency with the previously selected trajectory of the vehicle, and tracking of the reference center lane. All these costs are detailed in Section 2.2. Once these criteria are costed and merged into a weighted global cost function, the chosen trajectory is the best one with the lowest cost. Finally, a set of points, defined by the curvilinear abscissa, x and y coordinates, velocity, and curvature, depicts the best trajectory. For more details on the planning algorithm and its implementation, please refer to [START_REF] Said | Local trajectory planning for autonomous vehicle with static and dynamic obstacles avoidance[END_REF].

Cost function definition

The cost function is a weighted combination of costs that should be minimized in order to find the best trajectory among the many candidate trajectories generated by the planning algorithm. These costs are:

1. Smoothness cost C ρ [i]
This cost seeks to prefer lower-curvature trajectories. On the other hand, soft smoothness causes increased lateral acceleration, which affects passenger comfort. As a result, the integration of the curvature squared along a trajectory is chosen as a smoothness requirement for this trajectory in order to reduce its slackness:

C ρ [i] = ρ i M 2 ds i M (1)
Figure 2: Trajectory planning: simulation environment where ρ M and ds i M are the curvature and the curvilinear abscissa of a point M along the candidate trajectory i, respectively.

Cost of tracking the reference trajectory C r [i]:

In order to position the ego vehicle in the center of the intended lane, the vehicle must follow the global reference trajectory. The cost of following the reference trajectory is related to the offset since the candidate trajectories are created around this lane center with a given lateral offset (Fig. 2).

The reference cost in multi-lane structured environment is determined by the vehicle's position in relation to the lanes. So, the reference cost is given by:

C r [i] = (q i f -q re f c ) 2 , ( 2
)
where q i f is the lateral shift of the permanent phase of the candidate trajectory i to the host lane, q re f c is the lateral offset of the center of the reference lane. By default, the reference lane is the one in which the vehicle is located: it is either the host lane (q re f c = 0) or adjacent lane (q re f c = q L2 ). However, in case of returning to the host lane in an obstacle overtaking scenario, the host lane becomes the reference if the vehicle passes the obstacle (dotted lane in Fig. 3) and the host lane is navigable. Fig. 3 shows the different cases for the definition of the reference lane.

Consistency cost C c [i]:

A quick shift in the selected trajectory between two iterations needs a considerable control effort and can affect the vehicle's stability. As a result, a measure of similarity between the prior trajectory and the current candidate trajectory is required to avoid generating a trajectory that is drastically different from that of the previous step. Our technique reduces it to the difference between the current trajectory index i and the chosen one from a prior planning iteration.

C c [i] = |i -i t-1 des | (3) 4. Safety cost C s [i]:
We define two safety costs: a longitudinal cost C sl and a lateral cost C sL . The first one assesses the trajectory's navigability and the distance to the first obstacle. This cost increases when the first detected obstacle is longitudinally close and decreases when the obstacle is far (Fig. 4-left). The second cost considers the navigability of nearby trajectories and attempts to position the vehicle laterally in the middle of the navigable zone. The collision distance d obs (distance to the first obstacle) is transformed into a safety cost as follows:

C sl [i] = 0 obstacle-free trajectory 2 - 2 1+e -c 1 d i obs otherwise (4) g[k] = 1 √ 2πσ exp -(k ∆q) 2 2σ 2 (5) C sL [i] = (∑ j∈Γ i C sl [ j] g[i -j]) /n Γ i (6)
where c 1 is a cost setting parameter, ∆q is the desired lateral sampling resolution, σ is the standard deviation of the discrete Gaussian convolution for the risk of collision and is calculated by considering that g[ L v ∆q ] = 0.5, where L v is the vehicle width (Fig. 4-right). Γ i is the set of indexes of generated trajectories without i, n Γ i is their number.

Total cost function C T [i]:

We use the normalization technique to make the different cost criteria dimensionless because they are assessed in different units. As a result, we use the following formula:

C . [i] = C . [i] -min(C . ) max(C . ) -min(C . ) (7)
Then, the total cost function C T [i] of a given trajectory is defined as the weighted sum of the normalized cost functions defined above:

C T [i] =W ρ C ρ [i] +W c C c [i] +W r C r [i] +W sL C sl [i] +W sl C sL [i] (8)
with W ρ , W c , W r , W sL and W sl the weighting coefficients of smoothness, consistency, reference tracking, longitudinal safety and lateral safety costs respectively.

3 Impact analysis and tuning of weighting coefficients

Proposed method

Our goal is to understand better each cost influence on choosing the best trajectory and to find a generic weighting coefficients range for the cost function of the planning algorithm. First, we propose scenarios with or without static or dynamic obstacles. Second, we generate n random combinations of weighting coefficients, each weight is bounded between two values (0 ≤ l b ≤ W . ≤ u b ≤ 1). The random method was used due to the high number of combinations in a classic grid search method with five weights. Then, we run the planning algorithm n times. The n vehicle's trajectories obtained are classified into two categories:

• T c complete vehicle's trajectory set: the algorithm manages to reach the end of the proposed scenario without stopping or colliding with any obstacle.

• T i incomplete vehicle's trajectory set: the algorithm fails to reach the end of the proposed scenario due to hitting an obstacle or lane's borders.

The correlation between the obtained results and the corresponding combinations is then analyzed according to various factors such as the distance traveled, the longitudinal and the lateral acceleration, the reference tracking, the safety distance, and the overall behavior (conservative or aggressive). Based on this analysis, the wide range chosen first is contracted. this procedure is repeated until we obtain the desired generic range.

Case study

To study the proposed method, we apply the planning algorithm on a portion of a realistic trajectory track taken from SCANeR T M studio simulator, with different scenarios, for 50 random combinations. This trajectory and the corresponding limit speed and curvature are shown in Fig. 5. Initial speed may vary depending on the scenario. The trajectory planning module is implemented in an autonomous driving system composed of main modules for navigation: localization, perception, trajectory planning, control, and vehicle model. The two last modules are running at 50 Hz while the others are running at 10 Hz. The Renault-Zoe robotic vehicle model validated on SCANeR T M studio is used.

The presented scenarios are simulated using the combinations of weighting coefficients presented in Fig. 6. We start with simple scenarios and large ranges for weighting coefficients. The smoothness, Figure 5: Trajectory, longitudinal speed and curvature (same legend as Fig. 2) reference, and consistency coefficients ranges are chosen to be between 5% and 40% as they are used for performance purposes,while the coefficients ranges of safety costs are greater: The longitudinal safety one is between 5% and 65%. The lateral safety one is between 5% and 50%. These ranges are chosen based on evident knowledge about the algorithm and preliminary simulations. Safety is prior to performance while deciding the trajectory for autonomous navigation. Outside of these ranges, the planning algorithm does not give reasonable trajectories that can be executed. It will either give an extremely smooth with high smoothness weight or aggressive trajectory with high reference weight that will try to follow the lane's center. In the following, for each scenario, we illustrate the generated trajectories for all combinations and the final reached points for the same time horizon in order to analyze the results.

Scenario 1: static obstacle overtaking

The scenario presented in Fig. 7 consists of an overtaking of two static obstacles, each one on a lane sep-arated by a distance of 15m. This distance is considered small, taking into account the lane's desired speed (65km/h), and for this, the algorithm pushes the vehicle to brake when it slaloms between obstacles. In some combinations, where the sum of smoothness and consistency weights is significant (Set 2: W ρ > 12% or W c > 13% ), the vehicle crosses the trajectory with a high lateral acceleration value (> 4m/s 2 ) as shown in set 2 in Fig. 8, where set 1 represents trajectories with the remaining combinations. The smoothness cost leads the vehicle to move Figure 8: Scenario 1: lateral acceleration straight towards the obstacle, especially if the longitudinal safety weighting coefficient is too small. In contrast, the consistency cost delays the lane change until the last possible moment before reaching the safe stop distance, and this results in strong lateral acceleration. Then the algorithm proceeds to additional braking to lateral We note here that the developed algorithm is able to execute the trajectory safely, even with small safety weighting coefficients. As we can see in Fig. 7, all trajectories can reach the end of the scenario, and with almost a 10m difference between their ends due to additional braking on some trajectories (set 2). However, we can note that the trajectories with high lateral safety (W sl > 40%) and reference (W r > 20%) weighting coefficients execute the longest distance and hence offers the highest longitudinal speed.

Scenario 2: Mobile obstacle overtaking

The second scenario is a simple overtaking of a mobile obstacle at 40km/h where the vehicle speed during overtaking reaches 50 to 60km/h. As we can see in Fig. 9, all trajectories can reach the end of the scenario.

Figure 9: Scenario 2: trajectories in a mobile obstacle overtaking (same legend as Fig. 2)

The difference between trajectories can be seen in Fig. 10 when overtaking and returning to the host lane. The security inter-distance is respected in all combinations but leads to more or less conservative behavior depending on the reference weighting coefficient. For all combinations, the vehicle starts overtaking before 35m, which is the security distance corresponding to its velocity. However, the overtaking starts earlier with a high lateral safety weight. However, when returning to the lane, the inter-distance of the first trajectory that reaches the host lane's center is equal to 22m, which is equal to 2 seconds of obstacle's speed. Set 4 in Fig. 10 represents trajectories with smoothness weight W ρ higher than > 30%. In this set, the three trajectories with a late return to the lane have in addition a low reference weight W r < 7%, while the rest of the trajectories try hard to follow the center of the track. Furthermore, the reference cost acts on the vehicle stability: a high reference weight with low smoothness or consistency weight (set5: W r > 15% and (W ρ < 6% or W c < 6%) ) leads to a higher lateral acceleration compared to other combinations (see Fig. 11) and may result in a passenger's discomfort, due to late change lane and rapid return to the host lane. Note that set 3 represents trajectories with the remaining combinations. As we can see in Fig. 12, 92% of the generated trajectories can complete the scenario. The incomplete trajectories fail at the tight passage: the classification area around the trajectory, which is a little larger than the vehicle's imprint, detects hitting one of the obstacles. The four failed combinations have low lateral safety weights ( W sl < 10%) value and high reference weight (W r > 30%). This result is expected as the reference cost forces the vehicle to track the reference lane while the lateral safety cost tries to pass it in the middle of the navigable zone. represents the generated trajectories in the Frenet frame. We can see a significant difference between the trajectories. After t = 10s, we notice that a high reference and lateral safety weights (set 7: W r > 18% and W sl > 30%) lead to an early return to the host lane, while those with low reference weight and high smoothness weights (set 8: W r < 8% and W ρ > 25%) were late to return and overtake. In addition, the last one to return was the one with a high value of smoothness weight (W ρ = 39%) and caused over-steering when trying to follow the reference (s = 170m). At the overtaking of the mobile obstacle (t = 20s), the high lateral safety weight leads to earlier overtaking and then to an earlier return to the lane. The longest trajectories were those with the high lateral safety weight (W sl = 45% of the longest one). Note that set 6 represents trajectories with the remaining combinations. From these results, we can deduce that the smoothness, reference, and consistency weighting coefficients should be limited to relatively small values. Indeed, a high reference weight coefficient can be aggressive in a dual lane change scenario. Furthermore, high smoothness and consistency weights may cause a delay in response to the change in curvature. Consequently, bounded combinations are applied in the following simulation scenario.

Scenario 3.2: multi-maneuver scenario with limited combinations

After many tests and simulations, we decided to limit the ranges of weighting coefficients around the most favorable values. The new range proposed for combinations of weighting coefficients is between 5% and 10% for smoothness and consistency cost's weights, between 5% and 15% for reference cost's weight, and between 15% and 65% for longitudinal and lateral safety cost's weights. These new ranges are applied to the same scenario of 3.2.3. As we can see in Fig. 14, all trajectories can reach the end of the scenario, and the difference at the arrival point is less than 5m for different combinations of trajectories. In addition, it shows that all trajectories can pass smoothly and follow the reference without any accident. Note that these combinations are simulated for all the previous scenarios and succeed in reaching their

Scenario 3.3: longitudinal and lateral safety costs

To better understand the impact of the two components of the safety cost of the algorithm behavior, we simulate the third scenario with a fixed value for smoothness, reference, and consistency weights as 8%, 14% and 8% respectively, while the range of the two safety costs is between 15% and 65%. All combinations succeed in finishing the scenario smoothly with a lateral acceleration inside the

Discussion about influence correlation

results presented above give us a better understanding of the different costs that influence choosing the best trajectory. These simulations helped us to improve the planning algorithm. A good combination weight coefficients the vehicle to reference lane and to overtake smoothly and safely in the presence of obstacles, when overtaking is possible, without an abrupt lane change. Moreover, the planning algorithm respects the safety distances with all static and mobile obstacles and stops the vehicle if needed when no navigable trajectory is detected.

The smoothness cost, as defined, tries to move the vehicle straight forward as much as possible. Unfortunately, it delays the lane return with low reference cost or high consistency cost. The consistency cost acts as a damper of the vehicle's motion by delaying the variation of the vehicle's orientation. Concerning the reference cost, it plays an essential role, as presented before in all the scenarios. When encountering an obstacle, a high reference weight causes a high lateral acceleration, especially for the combination where the consistency or smoothness weight is low or the lateral safety weight is low. The safety cost guides the vehicle very well in the navigable zone and respects the inter-distances. In some cases, the costs identified above may correlate with each other. For example, in the case of lane following, the lateral safety cost correlates with the reference tracking cost by maximizing the lateral safety cost of candidate trajectories near the lane borders and thus minimizing those towards the lane's center. On the other hand, they behave differently when the lateral safety cost tries to position the vehicle in the middle of the navigable zone, away from the reference lane. The higher is the reference cost; the earlier is the host lane tracking. Finally, choosing the different weights of the cost function components must obey to a compromise between the other considered criteria. It is crucial to prioritize safety on performance. The present analysis shows that the proposed planning algorithm is robust against the cost weights variation. Moreover, this study allows us to identify the more suitable ranges for the different weights to arrive to a robust combination. Hence, the objective is to arrive to a planning algorithm tuning as generic as possible to the variation of driving conditions in a dynamic environment. The combination 8%, 14%, 8%, 40%, and 30% for smoothness, reference, consistency, longitudinal and lateral safety weights, respectively, is adopted as being a good compromise between safety and performance.

Conclusion

In this paper, after briefly the local trajectory planning method developed, we have introduced the method used to determine, analyze and fine tune the best ranges of the different weights of the trajectory planning cost function. The planning algorithm is tested using a variety of scenarios and ranges of weights combinations. This study shows the role of each cost in determining the best overall trajectory and leads us to select range for cost weight. We that the planning algorithm is robust to the variation of the cost weighting, which represents a significant advantage for encountering various driving situations and conditions in a dynamic environment, without the need for and tuning the algorithm.

Figure 1 :

 1 Figure 1: Local trajectory planning module

Figure 3 :

 3 Figure 3: Reference cost cases

Figure 4 :

 4 Figure 4: Safety cost

Figure 6 :

 6 Figure 6: Box plot of 50 combinations of weighting coefficients

Figure 7 :

 7 Figure 7: Scenario 1: trajectories in presence of 2 static obstacles (same legend as Fig. 2)

Figure 10 :

 10 Figure 10: Scenario 2: Frenet frame

Figure 12 :Figure 13 :

 1213 Figure 12: Scenario 3.1: trajectories (same legend as Fig. 2)

  smoothness weight leads to a late lane change, high a y and over/under steer W r high value of reference weight leads to an early lane change, high a y , prioritize reference tracking over safety W c high consistency weight leads to high a y (> 0.4g) W sL high lateral safety weight leads to longest path, early lane change;W sl longitudinal safety weight has a similar effect to lateral safety ;

Figure 14 :Figure 16 :

 1416 Figure 14: Scenario 3.2: trajectories (same legend as Fig. 2)

Table 1 :

 1 Cost analysis results


Arthur P C Lavenu

Juliette Lamarche

What controls diffuse fractures in platform carbonates? Insights from Provence (France) and Apulia (Italy)

Fractures are widespread in rocks and regional opening-mode arrays are commonly ascribed to major tectonic events. However, fractures occur in otherwise undeformed rocks. Some of these are earlydeveloped features independent of tectonics and forming a background network at regional scale.

To overcome this lack of understanding, two hydrocarbon reservoir analogues from platform carbonates have been targeted: the Provence (SE France), and the Apulian platform (SE Italy). In both areas, an early fracturing stage has been observed, made of high-angle-to-bedding opening-mode fractures, and bed-parallel stylolites. These features developed synchronously during the first burial stages and prior to major tectonic events. The fracture sets are not genetically related to the present-day layering. Contrarily, fractures developed in a brittle media where facies transitions were not sharp and did not act as mechanical discontinuities.

Carbonate facies distribution and early diagenetic imprint constrained the mechanical stratigraphy when fractures occurred. In addition, we observed that fractures related to late tectonic inversion were partly inhibited. Indeed, rock mechanical properties change through time. Characterizing the temporal evolution of carbonate rocks has revealed that diagenesis and sedimentary facies are the prime actors for brittleness and mechanical layering in carbonates.

Introduction

For more than a century, opening-mode fracture processes, mechanics, and origins have been investigated. In the first part of the 20th century, principles of fracture formation and propagation were established [START_REF] Griffith | The phenomena of rupture and flow in solids[END_REF][START_REF] Anderson | The Dynamics of Faulting[END_REF]. Since the early 1960s, extensive academic research reported in theses and books, outlined attributes and refined understanding of pattern development [START_REF] Price | Fault and Joint Development in Brittle and Semi-brittle Rock[END_REF][START_REF] Ramsay | Folding and Fracturing of Rocks[END_REF][START_REF] Suppe | Principles of Structural Geology[END_REF][START_REF] Engelder | Loading paths to joint propagation during a tectonic cycle: an example from the Appalachian Plateau[END_REF][START_REF] Atkinson | Fracture Mechanics of Rock[END_REF][START_REF] Ramsay | The techniques of modern structural geology[END_REF][START_REF] Pollard | Progress in understanding jointing over the past century[END_REF] and references therein; [START_REF] Price | Analysis of Geological Structures[END_REF]. Since the late 1960s, natural fractures have been increasingly important to the oil and gas industry owing to the role of fractures in reservoir behavior.

Researchers have worked at characterizing fractures in subsurface reservoirs [START_REF] Stearns | Reservoirs in fractured rocks[END_REF][START_REF] Nelson | Natural fracture systems: description and classification[END_REF][START_REF] Aguilera | Naturally Fractured Reservoirs[END_REF][START_REF] Narr | Origin of fracture porosityeExample from Altamont field, Utah[END_REF][START_REF] Nelson | Geological Analysis of Naturally Fractured Reservoirs[END_REF][START_REF] Narr | Estimating average fracture spacing in subsurface rock[END_REF][START_REF] Aydin | Fractures, faults, and hydrocarbon entrapment, migration and flow[END_REF][START_REF] Laubach | Practical approaches to identifying sealed and open fractures[END_REF][START_REF] Narr | Naturally Fractured Reservoir Characterization[END_REF][START_REF] Gale | Natural Fractures in shale: a review and new observations[END_REF], where fractures provide a pathway for fluid flow. Indeed, deciphering fracture origin and architecture in the subsurface would help in the understanding of reservoir behavior, which is critical in the oil and gas industry.

Fractures are widespread structures in rocks, and for practical applications it is imperative to accurately predict fracture attributes that affect fluid flow, including their size and spatial arrangement. Correlating the stress-strain relationship to fracture density in folded formations [START_REF] Murray | Quantitative fracture study; sanish pool, McKenzie county, North Dakota[END_REF][START_REF] Roberts | Curvature attributes and their application to 3D interpreted horizons[END_REF][START_REF] Keating | An experimental evaluation of the curvature-strain relation in fault-related folds[END_REF][START_REF] Hunt | On calibrating curvature data to fracture density: causes[END_REF] or faulted rocks [START_REF] Caine | Fault zone architecture and permeability structure[END_REF][START_REF] Aydin | Fractures, faults, and hydrocarbon entrapment, migration and flow[END_REF][START_REF] Billi | The damage zone-fault core transition in carbonate rocks: implications for fault growth, structure and permeability[END_REF][START_REF] Kim | Fault damage zones[END_REF][START_REF] Agosta | From fractures to flow : a field-based quantitative analysis of an outcropping carbonate reservoir[END_REF][START_REF] Matonti | Structural and petrophysical characterization of mixed conduit/seal falt zones in carbonates: example from the Castelas fault (SE France)[END_REF][START_REF] Michie | Influence of carbonate facies on fault zone architecture[END_REF]a r e among the most widely used concepts for fracture prediction. Where folds and faults are absent, regional opening-mode arrays are commonly ascribed to major tectonic events. However, fractures occur in otherwise undeformed rocks [START_REF] Laubach | Subsurface fractures and their relationship to stress history in East Texas basin sandstone[END_REF][START_REF] Frost | Controls on syndepositional fracture patterns, Devonian reef complexes, Canning Basin, Western Australia[END_REF]. Some of these are early-developed features independent of tectonics that form a background network at regional scale. Their occurrence and architecture remains poorly known and, as a consequence, largely unquantified in reservoir modeling, especially in the case of carbonate rocks. A fundamental, but frequently overlooked aspect of fracture spatial arrangement is how fracture occurrence partitions within a stratigraphic sequence, in other words, the fracture stratigraphy [START_REF] Laubach | Mechanical and fracture stratigraphy[END_REF]. Many outcrops have distinct segregations of fracture occurrence corresponding to depositional facies and interbeds, and in some cases these patterns correspond to current differences in mechanical properties (e.g. [START_REF] Corbett | Fracture development and mechanical stratigraphy of Austin Chalk, Texas[END_REF]. The prevalence of lateformed, barren joints to cluster within specific beds may mask (and in many cases, probably overprint) important patterns in which early-formed fractures respond to earlier states of mechanical stratigraphy. Here we describe fracture array spatial distribution that differs from the modern mechanical stratigraphy, having unbounded fracture height patterns likely reflecting early mechanical homogeneity. Differences in height patterns are potentially important variables in governing vertical migration of fluids.

Diffuse d or regional, background d fractures [START_REF] Guiton | Contribution de la fracturation diffuse al ad eformation d'ensemble durant le plissement de roches s edimentaires[END_REF][START_REF] Bazalgette | Relations plissement/fracturation multi echelle dans les multicouches s edimentaires du domaine elastique/fragile: Accommodation discontinue de la courbure par la fracturation de petite echelle et par les articulations[END_REF][START_REF] Vitale | Modes and timing of fracture network development in poly-deformed carbonate reservoir analogues, Mt. Chianello, southern Italy[END_REF] are opening-mode fractures that range from barren joints to fully mineral filled veins. One set of such an array in the subsurface can comprise fully cement-filled veins in one location or bed, and apparently barren (albeit usually cement lined) joints in adjacent locations or beds; moreover, small fractures of a given set typically qualify as veins whereas larger fractures may fitad e finition of joints (e.g., [START_REF] Laubach | Practical approaches to identifying sealed and open fractures[END_REF][START_REF] Gale | Predicting and characterizing fractures in dolostone reservoirs: using the link between diagenesis and fracturing[END_REF][START_REF] Lander | Insights into rates of fracture growth and sealing from a model for quartz cementation in fractured sandstones[END_REF]. Arrays are low-strain features that do not require tectonic events to form [START_REF] Olson | Natural fracture characterization in tight gas sandstones: integrating mechanics and diagenesis[END_REF]. Diffuse arrays can arise throughout cycles of burial and uplift by processes such as regionally elevated pore fluid pressure (for example, by gas generation: [START_REF] Fall | Natural hydraulic fracturing of tight-gas sandstone reservoirs, Piceance Basin, Colorado[END_REF] or thermoelastic contraction during uplift [START_REF] English | Thermomechanical origin of regional fracture systems[END_REF]. In carbonate rocks, an important but little appreciated type of diffuse fracture array may form in the first stages of the host rock brittle history. Owing to early and near-surface diagenesis, and consolidation and early embrittlement of some carbonate rocks [START_REF] Mcneill | Early load-induced fracturing in a prograding carbonate margin[END_REF][START_REF] Frost | Controls on syndepositional fracture patterns, Devonian reef complexes, Canning Basin, Western Australia[END_REF], diagenesis can play an essential role in fracture formation.

This paper provides insights on processes governing the occurrence, spacing, and height distribution of diffuse fractures in carbonate rocks. Synthesizing and comparing results from several recent case studies [START_REF] Lamarche | Relationships between fracture patterns, geodynamics and mechanical stratigraphy in Carbonates (South-East Basin, France)[END_REF][START_REF] Lavenu | Tectonic versus diagenetic origin of fractures in naturally fractured carbonate reservoir analogue: the Nerthe anticline (southeastern France)[END_REF][START_REF] Lavenu | Relating background fractures to diagenesis and rock physical properties in a platformeslope transect. Example of the Maiella Mountain (central Italy)[END_REF]2015), we identify commonalities that govern the close spacing and unbounded height distributions that characterize diffuse fracture arrays in some early fractured carbonate rocks.

Methods

Outcrops have been selected away from fold hinges and fault zones, in order to restrict our findings to diffuse fractures only. For each outcrop, a six-step workflow was applied to fully characterize the fracture network and the host rock. The first steps consist of gathering hard data: fracture parameters (strike, dip, size, spacing, cross-cutting relationships, cement fill …) measured using scanlines [START_REF] Priest | Estimation of discontinuity spacing and trace length using scanline surveys[END_REF][START_REF] Priest | Discontinuity Analysis for Rock Engineering[END_REF][START_REF] Zeeb | Evaluation of sampling methods for fracture network characterization using outcrops[END_REF]o r scan area method d or window sampling (e.g. [START_REF] Wu | An experimental study of the relationship between joint spacing and layer thickness[END_REF][START_REF] Hardebol | DigiFract: a software and data model impletmentation for felxible acquisition and processing of fracture data from outcrops[END_REF][START_REF] Lavenu | Background fractures in carbonates: inference on control of sedimentary facies, diagenesis and petrophysics on rock mechanical behavior. Example of the Murge Plateau (southern Italy)[END_REF]. The difference between the two techniques is the representativity of the sampling: the scanline method measures fracture attributes along a tape, in 1D, and the scan area method measures fractures inside a defined area (in this case, a rectangle shape), in 2D. Fractures have been analyzed in terms of geometry, chronology between the fracture sets and with stylolites, to determine cross cutting relationships. Rock samples were collected from each outcrop for further laboratory investigations. Petrophysics measurements were made to capture elasticity and porosity for each collected sample. V p was measured in free air using a Pundit7 apparatus (©CNS Farnell ltd.) that calculates the Pewave travel time through opposite sides of a sample. Porosity is measured on 2.54 Â 2.54 cm (1 Â 1 in.) dry plugs using a helium porosimeter. Petrographic analyses were used to unravel the sequence of diagenetic phases in the host rock, using cathodoluminescence techniques [START_REF] Amieux | La cathodoluminescence dans les roches s edimentaires[END_REF][START_REF] Boggs | Application of Cathodoluminescence Imaging to the Study of Sedimentary Rocks[END_REF] and diagenetic stratigraphy concepts [START_REF] Meyers | Carbonate cement stratigraphy of the Lake Valley formation (Mississippian) Sacramento mountains, New Mexico[END_REF][START_REF] Loreau | Diagenetic stratigraphy of discontinuity surfaces; an application to paleo-environments and sequence stratigraphy[END_REF][START_REF] Boggs | Petrology of Sedimentary Rocks[END_REF]. To perform the diagenetic sequencing, we used a Technosyn Cold Cathode Luminescence 8200 Mk II coupled to an Olympus™ BH2 microscope and to a Zeiss™ MRc5 digital camera. The last step of the workflow is creating a subsidence curve for specific areas taking account of regional geology and back-stripping [START_REF] Allen | Basin Analysis. Principles and Applications[END_REF], to define the burialeuplift history of the host rocks. For this study, the datasets used consists of about 3500 fracture measurements, 230 samples and 280 thin-sections, collected from 14 outcrops.

We account for all available geologic attributes, both internal to the rock, such as lithology, porosity, fracture orientation, density, diagenetic phases, and external, for example structural position, regional stress, geodynamic history, to determine critical parameters that rule fracturing and mechanical property evolution.

Geological settings

As mentioned above, each area has been described in detail elsewhere. Here we briefly summarize the geology to provide context for our synthesis of commonalities amongst our examples.

The Provence

The Urgonian Formation is a 400 m-thick platform carbonate succession, deposited in Valanginianeearly Aptian times [START_REF] Masse | Les calcaires urgoniens de Provence (Valanginien-Aptien Inf erieur) -Stratigraphie, pal eontologie, pal eoenvironnements et leur evolution[END_REF]. At that time, Provence was entirely within a tropical inner platform environment and sediments deposition evolved little throughout the region and the time interval such that the original depositional facies are essentially the same in all outcrops. We examined six outcrops of Urgonian carbonates distributed from the north to the south across the region (Fig. 2). Four outcrops are located far from folds and major faults to avoid tectonic-related fractures. Two outcrops are on the flanks of the Nerthe anticline, far from the hinges, to test the fracture/fold relative timing. One outcrop (B) is presented in Fig. 3 as an illustration.

In A, B, and C (Fig. 2), the exposure is composed of highly compacted packstones-to-grainstones containing numerous rudists, miliolids, bioclasts, and rare stromatolitic laminae but with many bed-parallel stylolites (8 per meter, Table 1) within layers and at layer interfaces. Beds thicknesses range from 0.5 to 10 m. In D, E, and F, the exposed outcrops are made of porous to microporous grainstones with rudist clasts and miliolids, and with few bedparallel stylolites (about 2 per meter, see Table 1). Beds range in thickness between 0.5 and 6 m (see [START_REF] Lamarche | Relationships between fracture patterns, geodynamics and mechanical stratigraphy in Carbonates (South-East Basin, France)[END_REF][START_REF] Lavenu | Tectonic versus diagenetic origin of fractures in naturally fractured carbonate reservoir analogue: the Nerthe anticline (southeastern France)[END_REF]. The studied outcrops are characterized by gradual facies transitions in thick beds. Each individual bed interface is enhanced by stylolites.

Each outcrop exhibits a pattern made of two non-stratabound vertical-to-bedding fracture sets, striking about NeS and EeW, ±30 . Fractures in both sets range from 0.1 to 10 m in height, and both sets are perpendicular to bed-parallel stylolites.

The Maiella Mountain

Maiella Mountain belonged formerly to the Apulian platform realm but is now a thrust-related anticline, NeS-oriented and eastward verging. The sequence exhibits an uninterrupted platform-to-basin profile of Late Jurassic to Late Cretaceous age [START_REF] Accarie | Dynamique s edimentaire et structurale au passage plate-forme/ bassin[END_REF][START_REF] Vecsei | Aggradation und progradation eines karbonatplattform-randes: kreides bis Mittleres Terti€ ar der Montagna della Maiella, Abruzzen. Mitteilungen aus dem Geologischen Institu der Eidgen€ osseischen Technischen Hochschule un der Universit€ at[END_REF]. Three outcrops were selected in the platform domain (Fig. 4). The investigated carbonate formation, socalled Lower Cretaceous Morrone di Pacentro Formation [START_REF] Bruni | Uppermost Jurassic-Lower Cretaceous carbonate deposits from Fara San Martino (Maiella, Italy): biostratigraphic remarks[END_REF], comprises 1500 m-thick inner platform deposits characterized by transitional facies variations. The observed texture gradually shifts from mudstone-wackestone with fenestrae to pelletoidal wackestone-packstone, characteristic of an inter-tosupratidal environment. Beds range in thickness between 0.5 m and 2 m. Some high-energy facies such as oolitic grainstones have been observed albeit very sparsely. The layers dip 35 eastward (Fig. 5). Fractures are vertical-to-bedding but no fracture abutment on present-day bed is evident. They are clustered into four sets: one major set N070 -oriented, 3 minor sets N030 ,N 1 2 0 and N160 (see [START_REF] Lavenu | Relating background fractures to diagenesis and rock physical properties in a platformeslope transect. Example of the Maiella Mountain (central Italy)[END_REF]. Fractures range from 0.1 to 10 m apart. Bed-parallel burial stylolites are ubiquitous (average of 12 per meter, Table 1), and are located either within layers or at layer interface (Fig. 5). Fractures and bed-parallel stylolites are mutually intersecting, the latter emphasizing the present-day bedding.

The Murge Plateau

The Murge area is a NWeSE-oriented continental plateau along the southeastern shore of Italy. Like the Maiella Mountain, the Murge Plateau is part of the former Apulian platform but is in the undeformed Apennine foreland. It is made of two slightly deformed superimposed carbonate platform successions: the Lower Cretaceous-to-Cenomanian Calcare di Bari Formation [START_REF] Borgomano | The upper cretaceous carbonates of the gargano-murge region, southern Italy: a model of Platform-to-basin transition[END_REF][START_REF] Spalluto | Stratigraphy of the mid-creataceous shallow-water limestones of the Apulia carbonate platform (Murge, Apulia, southern Italy)[END_REF], and the Coniacian-to-early Campanian Calcare di Altamura Formation (Luperto [START_REF] Sinni | Le Cr etac e sup erieur des Murges sudorientales (Italie m eridionale): stratigraphie et evolution des pal eoenvironnements[END_REF][START_REF] Borgomano | The upper cretaceous carbonates of the gargano-murge region, southern Italy: a model of Platform-to-basin transition[END_REF]. Five sites have been investigated (Fig. 6).

On outcrop, the Calcare di Bari Formation is made of thin mudstone intervals (20 cm to 1 m-thick) with levels of stromatolitic laminae, evolving to wackestoneepackstone (one to 2 m-thick) with pelletoids, ostracods and foraminifers. The Calcare di Altamura Formation is defined at outcrop scale by an alternation of massive wackestones with benthic foraminifera and laminated mudstones, in beds ranging from 5 to 70 cm thick. Some rudist shells and bioclasts are present. These facies are characteristic of a peritidalto-shallow subtidal depositional environment. In both formations, the facies changes are gradual rather than sharp.

In each outcrop, the fracture pattern is made of two vertical-tobedding sets, perpendicular to and mutually intersecting with bedparallel burial stylolites (Fig. 7). Fractures are not d or only rarely d stratabound and range between 0.2 and 5 m in height. Fracture sets are NeS and EeWt oN 0 4 5 and N135 in the Calcare di Altamura Formation, and N070e160 in the Calcare di Bari Formation [START_REF] Lavenu | Background fractures in carbonates: inference on control of sedimentary facies, diagenesis and petrophysics on rock mechanical behavior. Example of the Murge Plateau (southern Italy)[END_REF]. Bed-parallel stylolites are common with a density ranging from 5 per meter in the Calcare di Altamura Formation to about 11 per meter in the Calcare di Bari Formation (Table 1). Bedding is clearly marked, but the present day mechanical stratigraphy does not fit with the fracture stratigraphy. Looking closely, the current bedding is enhanced by pressure solution, meaning that either the early bedding has undergone dissolution or there was no bedding at all.

Results

General fracture stratigraphy

The studied outcrops display high-angle to vertical-to-bedding opening-mode fractures. These fractures are arranged in at least two perpendicular sets, and as many as four in some locations. Fracture heights do not correlate to the present-day bedding. In other words, fractures that strike parallel to each other have a wide range of heights and fractures terminate vertically at different positions within the unit, not all at the same bed boundary or interface. Thus fractures are not stratabound and either cut through bed boundaries or bypass interfaces (Fig. 8). In the terminology of [START_REF] Hooker | Fracture-aperture size-frequency, spatial distribution, and growth processes in strata-bounded and non-stratabounded fractures, Cambrian Mes on Group, NW Argentina[END_REF], the fracture stratigraphy is unbounded.

We interpret the unbounded height pattern to result from fracture development in a rock having little or no mechanical properties contrast. Although at present-day, bed interfaces are visible, which should constrain the fracture stratigraphy according to [START_REF] Shackleton | Evidence for temporally changing mechanical stratigraphy and effects on joint-network architecture[END_REF]. The lack of fracture abutment on any interface means that any mechanical stratigraphy that prevailed at the time of fracturing was insufficient to arrest fracture growth, and that the host rock mechanical contrasts then differed from those existing today.

In addition to high-angle-to-bedding fractures, bed-parallel stylolites are present both within beds and at bed interfaces. These stylolites are mutually cross cutting with the high-angle-tobedding fracture sets. According to [START_REF] Nelson | Natural fracture systems: description and classification[END_REF], fractures in rocks depict the stress state at the time of fracturing. However, the mutual cross cutting relationships between the different fracture sets does not allow clear discrimination relative to timing amongst fractures in different orientations. The Griffith theory of rupture (1921) and [START_REF] Anderson | The Dynamics of Faulting[END_REF] relate mode I failure process to a maximum vertical compressive s 1 . In the case of tilted outcrops (B and C in Fig. 2, and A, B, and C in Fig. 4), a back tilting procedure has been applied (see stereoplots in Figs. 3 and5). This procedure restores the layers in their original horizontal position, as well as the features contained within the layers. Back tilting the currently tilted layers to horizontal shows that the fractures stood at high-angle to vertical to bedding, and the bed-parallel stylolites were horizontal. The cross cutting relationships between both features combined with patterns in the restored state is evidence of pretilting/folding fracturing, when bedding laid horizontally, and under a stress pattern with a vertical s 1 .

Following [START_REF] Nelson | Significance of fracture sets associated with stylolite zones[END_REF], a vertical s 1 is consonant with the state of stress that leads to burial pressure solution. Thus, verticalto-bedding fractures and bed-parallel stylolites developed under a similar stress state, most probably synchronously, and during the burial of the host rock. The synchronicity of fractures sets with stylolites argues for a stress pattern defined by a vertical loadinduced s 1 , and poorly differentiated horizontal components (s 2 and s 3 ) with low amplitudes (e.g. [START_REF] Laubach | Characteristics and origins of coal cleat: a review[END_REF][START_REF] Lamarche | Relationships between fracture patterns, geodynamics and mechanical stratigraphy in Carbonates (South-East Basin, France)[END_REF].

Petrophysical properties

Analyzing the petrophysical properties of host rocks reveals similarities in terms of pore types and rock susceptibility to fracturing across the investigated areas. Two main pore types have been identified: intergranular porosity and microporosity. In the case of intergranular porosity, porosity values are low: an average of 2% in the Calcare di Bari Formation, close to 3% in the Urgonian of South Provence, and up to 7% in the Morrone di Pacentro Formation (Table 1). In the specific case of intergranular porosity, when values are low, P-wave velocity is high: about 6250 m s À1 in the Calcare di Bari Formation, 5900 m s À1 in the Urgonian, and 5400 m s À1 in the Morrone di Pacentro Formation, with a negative correlation trend between the 2 parameters (Fig. 9). In the case of microporosity, values range from 1 to 21%, with an average value close to 11% in the case of the Urgonian of North Provence (see Table 1 and Fig. 9). Similarly, but with a narrow scatter, porosity values of the Calcare di Altamura Formation range between 0.1 and 11% (Fig. 9). Such variation in porosity results in a wide scatter of P-wave velocities: between 2100 and 6300 m s À1 with an average value close to 4940 m s À1 in Northern Provence, and between 4760 and 5900 m s À1 in the Calcare di Altamura Formation (Table 1 and Fig. 9).

Intuitively, the fracture pattern appears to be intrinsically related to the host rock pore type. For instance, in the Urgonian carbonates of Provence, the northern part is dominated by microporosity and the southern part by intergranular porosity. This latitudinal dichotomy is reflected in the variation on fracture and stylolite densities: to the north, the host rock bears few to no stylolites and scarce fractures, while to the south the fracture frequency is higher with dense burial stylolites (see Table 1; see Fig. 3. Photograph of one of the studied outcrop in Provence, here on the southern flank of the Nerthe anticline (outcrop B on Fig. 2). The bottom schematic of the picture is showing the interpretation of the photograph and the stereoplots (lower hemisphere) are presenting the fractures and bedding measurements in the present-day and back-tilted position.

Table 1

Summary of the main host rock parameters, for each of the investigated area, with the location, the name of the formation, the depositional environment, the age of the rock, the main pore type, the porosity range, the V p range, the fracture/stylolite density and the maximum estimated burial depth. [START_REF] Lamarche | Relationships between fracture patterns, geodynamics and mechanical stratigraphy in Carbonates (South-East Basin, France)[END_REF]. However, the rock susceptibility to fracturing not only relates to the host rock pore type but also to the volume of porosity. Indeed, the northern part of Provence is made of porous carbonates (11% in average) in contrast to the southern part, made up of tight carbonates (3%, see Table 1). Looking at the Murge Plateau, where investigated formations are low porosity (or tight) d 2 and 4% respectively (see Table 1) d the highest fracture density is observed in the most porous host rock, where microporosity is one of the main pore types (i.e. Calcare di Altamura Formation).

For the case of Morrone di Pacentro Formation, collected samples have been categorized in two groups according to the dominant pore type. Most of the samples have fracture porosity, and in minor proportion, a matrix intergranular porosity. Since fracture porosity is difficult to measure, the values in Table 1 only depict the matrix porosity of the sample. Based on observation, samples affected by fracture porosity can be differentiated from intact samples. However, intact samples do not mean that the overall host rock is intact, but only that the collected sample does not contain visible fractures. The samples assigned to the "fracture-stylolite" group do have intergranular porosity, but this is occluded by calcite, which results in very low values of porosity (0.5%), high V p (6400 m s À1 ) and denser fracture pattern. Consequently, not only the pore type but also the associated cementation can drastically change the mechanical behavior of the rock.

Diagenetic history

Based on natural light microscopy and cathodoluminescence, a single diagenetic sequence has been deciphered with some variations from one investigated host carbonate to another (Fig. 10).

The major phases observed in the three areas, and the relative timing of these phases is as follows: (a) an early micritization phase affected grains and bioclasts producing either a micritic envelope or full micritization. The micritization process occurs at the watersediment interface, meaning a development in the very first stage of eogenesis [START_REF] Bathurst | Carbonate sediments and their diagenesis[END_REF][START_REF] Tucker | Carbonate Sedimentology[END_REF]. The micrite appears dusty grey to black in natural light (Fig. 10e7). In the Murge Plateau, a phase of dissolution has been identified and occurs after micritization (Fig. 10e7). (b) Then, an early calcite [EC] affected grains, developing as isopachously-distributed equigranular-tofibrous cement around grains at shallow depth during eogenesis (Fig. 10e1). Consequent to the development of the early calcite, sediments lithified and turned into a brittle media. This phase of early calcite tends to preserve grains from the effects of later compaction. In Provence, this fibrous calcite is present only in the case of porous rocks, and is preceded by another phase of early calcite resulting from a neomorphism of high-Mg calcite into low-Mg calcite characterized by irregular intercrystalline boundaries, and appears non luminescent in cathodoluminescence. Tight Urgonian carbonates suffered compaction though. According to L eonide et al. ( 2014), in Provence both cement phases (neomorphism and fibrous calcite) occurred near the surface and in a marine environment. In the investigated platform carbonates of Maiella Mountain, the fibrous Early Calcite is overlapped by pendulum prismatic calcite (Fig. 10e5). Although the first early calcite might have occurred in marine environment, the second early calcite morphology argues for a development in vadose conditions. In any case, petrographic evidences argue for a development of both cements EC1 and EC2 at very shallow burial [START_REF] Tucker | Carbonate Sedimentology[END_REF]. (c) The subsequent event observed in all the studied outcrops is a brittle event: at least one fracture phase [F/F1] developed coeval with (d) a drusy-to-mosaic blocky calcite spar [BS/BS1] which filled fractures, intergranular and intragranular pore spaces. In natural light, the BS crystals appear limpid to slightly cloudy and dull blue alternating with orange luminescent fringes in cathodoluminescence (Fig. 10e2, 4, 6 and 8). The petrographic evidence argues here for a burial-related origin of the blocky calcite [START_REF] Bathurst | Carbonate sediments and their diagenesis[END_REF][START_REF] Tucker | Carbonate Sedimentology[END_REF]. In Provence and Maiella, a second fracturing occurred: F2 fractures developed cross cutting F1 and BS1 (Fig. 10e4 and6). Concomitantly with F2 fractures, a second blocky spar [BS2] cement formed, with dull blue-toviolet colors under cathodoluminescence (Fig. 10e2 and6). This cement filled F2 fractures and residual intergranular spaces. Both phases of blocky spars are responsible for porosity reduction. (e) Bed-parallel burial-related stylolites formed synchronously with fractures and blocky spars, and affected all the other diagenetic phases (Fig. 10e3 and4). No relative chronology is clear between fractures and bed-parallel stylolites because there are an equal proportion of cross-cutting relationships. The events (c)t o( e) occurred during the mesogenesis (burial diagenesis). No diagenetic phases (cementation, dissolution or fractures) related to uplift or exhumation have been evidenced. Moreover, no fractures related to the recent tectonic events have been identified in any of the studied outcrops. Only very few of the early-formed fractures have been reactivated in shear. The present study has shown in all cases that during burial, blocky cement crystallized within fractures and in surrounding pores (see Fig. 10e2, 6, 7, and 8). Locally, vertical-tobedding stylolites related to late tectonic inversion have been observed. Although they are obvious in Maiella Mountain, they are negligible in Provence and in the Murge.

Early diagenetic processes are responsible for host rock early brittleness. As previously stated, the investigated platform carbonates are made of shallow water deposits that are characterized by gradual facies transitions without sharp mechanical limits. Early diagenetic processes converted sediments to brittle rock. Because early fractures have an unbounded fracture height pattern, we infer that in the early diagenetic state, the rock was extremely brittle, accounting for anomalously close fracture spacing, and interfaces were weak (or at least mechanically indistinct) accounting for the unbounded fracture height distribution. Subsequent burial cementation filled pores and strengthened previous porous parts of the rock. Because cementation and compaction acted differentially on originally mechanically similar depositional facies, late mechanical stratigraphy is more differentiated in terms of contrasts between layers and between layers and interfaces. This later diagenetic state may have be stronger and less fracture prone. Late diagenesis could also be responsible for strengthening and possibly sealing early fractures [START_REF] Laubach | Practical approaches to identifying sealed and open fractures[END_REF][START_REF] Olson | Natural fracture characterization in tight gas sandstones: integrating mechanics and diagenesis[END_REF][START_REF] Lavenu | Tectonic versus diagenetic origin of fractures in naturally fractured carbonate reservoir analogue: the Nerthe anticline (southeastern France)[END_REF].

Geodynamic evolution

Intuitively, with greater burial depth, stress is higher, and fracture-stylolite density is expected to be higher. Consideration of the geodynamic evolution of the rocks we studied shows that the relationship between structures and burial depth is not straightforward. In Maiella Mountain, the studied formation has been submitted to a simple burial-uplift history: a long-lasting burial from deposition time to Pliocene, with a maximum depth of 1.5 km, then a quick uplift during Pliocene (Fig. 11). The diagenetic sequence compared to the geodynamic evolution allows constraining a fracturing path. The analysis of diagenetic phases has revealed that early calcites developed at shallow burial depth, what led to early host rock brittleness, and allowed fracturing when the stress conditions were favorable. In addition, fractures developed vertical and perpendicular to bedding, and synchronously with bed-parallel stylolites. Therefore, fracturing occurred in a depth window between the seafloor and a few hundred meters' depth, synchronous with burial stylolites, during the pre-folding stage of the Maiella. In Provence, the burial-uplift trend is similar for all the outcrops but with variations in the burial depth (600 me110 0 m ) from one studied outcrop to another (Fig. 11). However, no correlation exists between burial depth and fracture-stylolite density in Provence [START_REF] Lamarche | Relationships between fracture patterns, geodynamics and mechanical stratigraphy in Carbonates (South-East Basin, France)[END_REF]. Contrarily, outcrops that have undergone the deepest burial (A and D, Fig. 2) do not have the same properties and fracture pattern. For instance, in outcrop A, the host rock is tight and fractured. In outcrop D, the host rock is porous and poorly fractured. In the same way, the outcrop that has been submitted to the shallowest burial (B) is made of tight fractured carbonates. Similarly to the Maiella, fractures initiated as a function of the local stress that prevailed at an early burial stage rather than relative to far-field tectonic events. The geodynamic evolution of the platform carbonates of the Murge plateau includes a very shallow burial (less than 500 m) and gentle uplift (Fig. 11 diagenetic history reveal that fracture and bed-parallel stylolites developed very early in burial history. Because of the synchronous development of fractures and burial stylolites, we assume that both features occur subsequent to a vertical maximum s 1 . The fact that fractures are clustered in two perpendicular mutually intersecting sets suggests a single brittle event, possibly involving biaxial extension [START_REF] Olson | Predicting fracture swarms d the influence of subcritical crack growth and the crack-tip process zone on jointspacing in rock[END_REF]. However, the formations differ in age, and the fracture pattern orientation varies from one to the other. Thus, fracturing occurred diachronously. We infer that fractures developed sequentially. Once host rock brittleness was acquired, at relatively shallow depth, vertical failure occurred following the local stress pattern that prevailed at this time. Furthermore, the diachronic fracturing in the Calcare di Altamura Formation did not involve the creation of new fractures in the underlying Calcare di Bari Formation. In fact, burial cementation (blocky spar) likely strengthened the host rock, reducing its susceptibility to fracture under mild uplift and extension.

In each case, the fracture-stylolite pattern argues for a shallow burial development of brittleness between few tens to hundreds of meter-depth, as soon as the initial sediments lithified d which can occur as fast as in a year's time [START_REF] Friedman | Rapidity of marine carbonate cementation d implications for carbonate diagenesis and sequence stratigraphy: perspective[END_REF]. The burial depth is not a stand-alone key for understanding fracture distribution. Although host rock fracturing likely initiated under the influence of burial-related vertical stress, drivers of mechanical behavior are the depositional facies, which control the initial mechanical properties and susceptibility to diagenesis, and the diagenetic processes that are responsible for the host rock embrittlement and inferred petrophysical properties.

Discussion

Modelling naturally fractured reservoirs (NFR) remains challenging. Historically, as a matter of simplification for the purposes of efficient flow simulation, fractured reservoirs have been considered as a dual porosity media: matrix þ fractures [START_REF] Warren | The behavior of naturally fractured reservoirs (SPE 426)[END_REF]. Nowadays, geological models can afford to model individual fractures in Discrete Fracture Networks (DFN), where fractures are populated in a deterministic way (Lei et al., 2017 and references therein). However, the real 3D architecture of fracture networks remains a key uncertainty. Geometry, scales, lateral and vertical variations need to be defined to fully characterize the fracture pattern in 3D. Input data to model a subsurface fracture network are limited to 1D well data and 2D/3D seismic data. Most of the time, these data are not appropriate for fracture characterization either because of a lack of spatial continuity or because of an under-size resolution of the considered object. In the case of carbonate reservoirs, fracture characteristics are strongly dependent on host rock heterogeneities that control mechanical properties and reservoir quality. Host rock lithologic heterogeneities are complex in carbonates owing to a wide range of depositional environments, facies, and diagenetic evolution. Platform carbonates in particular, are subject to early and strong diagenesis owing to rapid carbonate mineral cementation at near-surface and shallow burial conditions and the potential for subsequent differential compaction of initially highly porous, highly cemented rock [START_REF] Anselmetti | Controls on sonic velocity in carbonates[END_REF][START_REF] Fournier | Pore space evolution and elastic properties of platfrom carbonates (Urgonian limestone, Barremian-Aptian, SE France)[END_REF][START_REF] Fournier | Diagenetic patterns and pore space distribution along a platform to outer-shelf transect (Urgonian limestone, BarremianeAptian, SE France)[END_REF].

To assess the 3D geometry of fracture networks and the complexity of carbonate rocks, the best alternative remains analogue outcrops [START_REF] Di Naccio | Role of mechanical stratigraphy on fracture development in carbonate reservoirs: insights from outcropping shallow water carbonates in the Umbria-Marche Apennines, Italy[END_REF][START_REF] Bertotti | Improving predictions on mechanical stratigraphy of buried sedimentary successions: lessons and workflow from outcrop analogs[END_REF][START_REF] Bisdom | Calibrating discrete fracture network models with a carbonate three-dimensional outcrop fracture network: implications for naturally fractured reservoir modelling[END_REF], although the tectonic, diagenetic and mechanical histories may differ from subsurface reservoirs [START_REF] Laubach | Mechanical and fracture stratigraphy[END_REF][START_REF] English | Thermomechanical origin of regional fracture systems[END_REF]. In this paper, we focus on carbonate rocks of two platform realms cropping out in three different areas (see Fig. 1) to unravel the contribution of depositional facies and diagenetic processes to rock mechanical properties evolution through timeegeodynamic history, and how they control fractures development and architecture. In each case study, mutual cross cutting relations amongst sets defined by orientation and associated bed-parallel stylolites show that diffuse fracturing occurred as a single brittle event, early in the burial history, and synchronous with burial. Occurrence of diffuse fractures is not related to one parameter but to a combination of several interwoven parameters, brittleness acquisition, burial loading, and subsequent cementation and compaction serving to strengthen the rock (see Fig. 12).

Brittleness acquisition of carbonates

The diagenetic history undergone by all studied host rocks includes early calcitic cements developed during the first stages of burial. These cements caused lithification and host rock early embrittlement. In addition, the fracture geometry does not mimic the present-day mechanical stratigraphy, but rather reflects an early burial condition, different from today's mechanical stratigraphy. In these platform carbonates, diffuse fractures occurred at an early stage of the mechanical evolution of the host rock (see Fig. 12). Recently [START_REF] Friedman | Rapidity of marine carbonate cementation d implications for carbonate diagenesis and sequence stratigraphy: perspective[END_REF][START_REF] Eberli | Factors controlling elastic properties in carbonate sediments and rocks[END_REF] investigated modern platform carbonates in the Bahamas, revealing that lithification processes can occur in a very short time-span, from months to years d instantly on the geological time scale. In warm and protected shallow-water conditions, carbonate sediments have been rapidly subjected to early cementation: micritization processes modified the initial sediments at the seafloor level, and initiated the lithification.

The next cement phase occurring in all studied samples is an early calcite, developing as an isopachous fringe around micritized (1, 5, and7) andcathodoluminescence (2, 3, 4, 6, and8) of selected thin sections, to illustrate the main diagenetic phases.

elements (see Fig. 10). The development of these cements coated and bonded grains and clasts, forming a rigid framework (e.g. Bjorlykke and Hoeg, 1997). This reduced the primary porosity, increased the rock stiffness, and turned the sediments into a brittle rock. Evidence of compaction on cemented grains d i.e. concavoconvex and sutured contacts d argue for an early development of these calcites, before or in the first stages of the burial. In addition, the fibrous calcite morphology and distribution shows that diagenetic fluids were most probably of marine origin [START_REF] Tucker | Carbonate Sedimentology[END_REF], rather than of burial origin. A recent study on Urgonian carbonates of Provence (L eonide et al., 2014) using isotope geochemistry of 13 C and 18 O, has determined an Urgonian age d i.e. Late BarremianeEarly Aptian d and a marine origin of the fluids which led to isopachous fibrous calcites. These results confirm the very early origin of the calcite, which makes the rock brittle at very shallow depth, and the marine origin of the fluids in platform carbonates.

Thus, in the very first stages of the rock development and at shallow depth the host rock is prone to accommodate brittle deformation (e.g. [START_REF] Frost | Controls on syndepositional fracture patterns, Devonian reef complexes, Canning Basin, Western Australia[END_REF][START_REF] Mcneill | Early load-induced fracturing in a prograding carbonate margin[END_REF][START_REF] Vandeginste | Interplay between depositional facies, diagenesis and early fractures in the early cretaceous habshan formation, Jebel Madar, Oman[END_REF].

Diagenesis and mechanical stratigraphy

As previously stated, fractures develop as an unbounded diffuse network. They rarely terminate on any interface (i.e. bed boundaries, bed-parallel stylolites), meaning that at the time of fracturing, the rock medium was characterized by no or weak interfaces. Yet, mechanical limits were not sharp enough to arrest fractures. However, field outcrops show obvious present day bedding, commonly pressure solution enhanced. The bed boundaries are sharp and sudden, giving rise to strong mechanical breaks. Giving that fractures developed soon after rock embrittlement, the unbounded character of fracture height remains misunderstood: either bed boundaries had no control on fracture propagation, or fracturing occurred prior to present day bedding development. This question has been addressed in [START_REF] Lamarche | Relationships between fracture patterns, geodynamics and mechanical stratigraphy in Carbonates (South-East Basin, France)[END_REF] and [START_REF] Lavenu | Tectonic versus diagenetic origin of fractures in naturally fractured carbonate reservoir analogue: the Nerthe anticline (southeastern France)[END_REF], and the answer resides in the specific nature of the host rock. Indeed, while depositing, sediments are characterized by gradual facies transitions within thick beds rather than sharp facies changes. Early cements froze the initial settings, and as a consequence, the poor mechanical contrast of the original host rock (Fig. 12). Fractures initiated in this context, where facies transitions were not mechanical boundaries. Later on the apparent bedding developed as pressure-solution due to burial. Therefore, the fracture stratigraphy is strongly dependent on the host rock mechanical stratigraphy very soon after deposition [START_REF] Laubach | Mechanical and fracture stratigraphy[END_REF][START_REF] Lamarche | Relationships between fracture patterns, geodynamics and mechanical stratigraphy in Carbonates (South-East Basin, France)[END_REF], and on the pre-fracturing diagenetic events [START_REF] Di Naccio | Role of mechanical stratigraphy on fracture development in carbonate reservoirs: insights from outcropping shallow water carbonates in the Umbria-Marche Apennines, Italy[END_REF][START_REF] Gale | Predicting and characterizing fractures in dolostone reservoirs: using the link between diagenesis and fracturing[END_REF][START_REF] Barbier | Sedimentary and diagenetic controls on the multiscale fracturing pattern of a carbonate reservoir: the Madison Formation (Sheep Mountain, Wyoming, USA)[END_REF]. The fossil mechanical stratigraphy is fixed with early calcite soon after sedimentation. Actually, the fracture stratigraphy can be considered as an image of the fossil mechanical stratigraphy before burial.

Consequently, diagenetic processes are the main controlling factor for both early and late mechanical behavior of the platform carbonates. Eogenetic cementations are responsible for rock embrittlement and subsequent early fracture development at the onset of mesogenesis (Fig. 10). Later and deeper in mesogenesis, burial cementation is responsible for rock over-strengthening, which reduces the host rock susceptibility to fractures.

Inference of petrophysical properties with mechanical behavior

Although the depositional environment is common in all studied sites (platform carbonates), small differences can occur in the host rock. We noticed that the pore type and the facies distribution can influence petrophysics, mechanics, and thus, fractures.

At the very moment of deposition, the environment where sediments stacked rules the initial petrophysical properties: in deep water or in protected platform, where energy is low-tomoderate, micritic mud is the main component of the sediments and leads to low porosity/permeability, which slows diagenetic alteration [START_REF] Anselmetti | Controls on sonic velocity in carbonates[END_REF]. Where the energy moderate-to-high (in shoal or open-platform conditions for instance), sediments are mainly grain-supported. They have an initial higher porosity and permeability that allows marine and/or meteoric fluid circulation, and enhances diagenesis [START_REF] Anselmetti | Controls on sonic velocity in carbonates[END_REF]. Therefore, in platform conditions, carbonate sediments are prone to early cementation at shallow depth, which results in a pre-or early burial lithification [START_REF] Eberli | Factors controlling elastic properties in carbonate sediments and rocks[END_REF][START_REF] Mcneill | Early load-induced fracturing in a prograding carbonate margin[END_REF].

Although the sediments deposited in similar conditions in all the studied outcrops, the porosity and P-wave velocities vary as a function of the observed pore type, and so is the fracture density. For instance in the case of Urgonian platform carbonates, fracture density is higher in rock hosting intergranular porosity than microporosity (Table 1). According to L eonide et al. ( 2014), 2012), [START_REF] Lavenu | Tectonic versus diagenetic origin of fractures in naturally fractured carbonate reservoir analogue: the Nerthe anticline (southeastern France)[END_REF][START_REF] Lavenu | Relating background fractures to diagenesis and rock physical properties in a platformeslope transect. Example of the Maiella Mountain (central Italy)[END_REF][START_REF] Lavenu | Background fractures in carbonates: inference on control of sedimentary facies, diagenesis and petrophysics on rock mechanical behavior. Example of the Murge Plateau (southern Italy)[END_REF]. microporous Urgonian limestones did not experience porosity occlusion in the early stages of diagenesis. Contrarily in tight Urgonian limestones, porosity is occluded by early marine and/or meteoric cementation. Petrophysical measurements on the Urgonian samples confirm the conclusions of L eonide et al. (2014). Following [START_REF] Nelson | Significance of fracture sets associated with stylolite zones[END_REF] and [START_REF] Wennberg | Fracture intensity vs. mechanical stratigraphy in platform top carbonates: the Aquitanian of the Asmari Formation, Khaviz Anticline, Zagros, SW Iran[END_REF], tight rocks are more susceptible to brittle deformation, which is in agreement with our observations in the Urgonian carbonates of Provence. Similarly, in Maiella Mountain, the most fractured samples are the ones where initial intergranular porosity has been occluded due to early cementation. Where porosity has been preserved d and more recently filled by blocky spars d early fracture density is smaller. Hence, the initial pore type and the early cementation control the diffuse fracture architecture in platform carbonates (Fig. 12). Therefore, the fracture occurrence is not only conditioned by the pre-fracturing diagenetic events, but also by the rock properties evolution through time [START_REF] Voight | Stress history and rock stress[END_REF][START_REF] Narr | Origin of fracture porosityeExample from Altamont field, Utah[END_REF][START_REF] Shackleton | Evidence for temporally changing mechanical stratigraphy and effects on joint-network architecture[END_REF].

Stress conditions for fracturing

The studied outcrops exposed carbonates of different ages. These carbonates deposited in different paleogeographic locations, and suffered different and diachronous tectonic events. Nevertheless, in all outcrops we have observed a similar pattern made of early high-angle-to-bedding fractures. In the studied platform carbonates, the diagenetic analyzes realized on the collected samples revealed an early near-surface lithification, and subsequent embrittlement, which allows failure at the onset of mesogenesis. Fracture architecture and relationships with burial stylolites strongly argue for a compaction-related development of both features, during host-rock burial, with a maximum vertical s 1 and poorly differentiated horizontal s 2 and s 3 (e.g. [START_REF] Lavenu | Relations entre motif de fracturation et g eodynamique dans les roches carbonat ees. Importance du faci es de d epôt, de la diag en ese et des propri et es m ecaniques de la roche[END_REF]. Indeed, burial depth controls the vertical stress increase. The failure is reached when s 1 is sufficiently higher than s 2 and s 3 . Once burial was initiated and stress state reached failure, fracturing occurred despite the shallow depth, with fracture density depending on host rock petrophysical properties.

Fracture timing versus depth aids understanding of the effect of diagenetic processes on fracturing. The present study has given evidence for a sequential fracture development once the host rock brittleness is acquired, rather than a single brittle event affecting the sedimentary pile. The platform carbonate successions of the Murge Plateau, where we observe the exact same fracture-stylolite pattern as in the Urgonian carbonates and in the Maiella are in favor for this interpretation: the fracture set orientations are strictly different in the two superposed formations (i.e. Calcare di Bari Formation and Calcare di Altamura Formation). Both fractures and bed-parallel stylolites developed step by step after lithification and early burial of each layer. Hence, diffuse fracture occurrence depends upon early host rock mechanical properties, being related to early diagenesis, and to local stress variations rather than far-field tectonics. Therefore, the very first burial stages enhance the onset of fracturing, after early diagenetic processes made the rock brittle with favorable petrophysical and mechanical properties.

Conclusion

Through the presented case studies, the role of early diagenesis on host rock mechanical properties, and consequently on fracture development, timing, and spatial arrangement, has been emphasized. We have shown that most of the diffuse fractures which provide the background permeability in naturally fractured carbonate reservoirs developed in response to early burial stresses, and to diagenetic stiffening by early micritization and calcite cementation of carbonates. In the investigated platform carbonates, early diagenetic phases are enabling:

(1) Early brittleness acquisition at the onset of mesogenesis, which leads to pervasive, closely spaced arrays of openingmode fractures developed in flat-lying rocks far from folds and faults under small but sufficient vertical s 1. The diagenetic strengthening of the rock at shallow depths (0e500 m) makes failure possible for small magnitudes of s 1 , s 2 , s 3, and very young host rock (years to few million years). ( 2) Unbounded fracture height patterns and close fracture spacing unrelated to current mechanical limits. Immediately subsequent to deposition, rocks were largely undifferentiated mechanically, and interfaces and various beds (facies) that later were differentiated as a consequence of progressive burial diagenesis had little or no mechanical contrast. Later and deeper burial changes the apparent bedding by pressure-solution and compaction. Implicitly in fractured platform carbonate reservoirs, the diffuse fracture pattern cannot be predicted using usual methods based on (late) mechanical stratigraphy. (3) Diffuse fracture architecture, depending upon pore types and local s 2 /s 3 orientation. Diffuse fracture occurrence in platform carbonates and distribution is neither related to largescale deformation and tectonic inversion, nor to deep burial, nor the present-day host-rock mechanical stratigraphy.

Fig. 1 .

 1 Fig. 1. Schematic map of southeastern France and Italy with the location of the studied areas.

Fig. 2 .

 2 Fig. 2. Schematic structural map of the Provence with outcrop locations, modified from Lamarche et al. (2010).

Fig. 4 .

 4 Fig. 4. Schematic structural map of the Maiella Mountain with outcrop locations, modified from Lavenu et al. (2014).

  and Lavenu et al., 2015 for more details). Both fracture geometry and

Fig. 5 .

 5 Fig. 5. Photograph of one of the studied outcrop in Maiella Mountain (outcrop B on Fig. 4). The bottom schematic of the picture is showing the structural interpretation of the photograph and the stereoplots (lower hemisphere) are presenting the fractures and bedding measurements in the present-day and back-tilted position.

Fig. 6 .

 6 Fig. 6. Schematic structural map of the Murge Plateau with outcrop locations, modified from Lavenu et al. (2015).

Fig. 7 .

 7 Fig. 7. Photograph of one of the studied outcrop in Murge Plateau (outcrop E on Fig. 6). The bottom schematic of the picture is showing the structural interpretation of the photograph and the stereoplots (lower hemisphere) are presenting the fractures and bedding measurements in the present-day and back-tilted position.

Fig. 8 .

 8 Fig. 8. Schematic of the fracture architecture at outcrop scale, with a zoom showing the relationship between fractures, bed-parallel stylolites and facies variations.

Fig. 9 .

 9 Fig. 9. Porosity versus V P distribution for all the collected samples plotted area-wise for Provence, Maiella and Murge.

Fig. 10 .

 10 Fig. 10. Diagenetic sequences (left) for all the investigated outcrops, region-wise: Provence, Maiella Mountain and Murge. Photomicrographs (right) in natural light(1, 5, and 7) and cathodoluminescence (2, 3, 4, 6, and 8) of selected thin sections, to illustrate the main diagenetic phases.

Fig. 11 .

 11 Fig. 11. Burial-uplift curves for each studied outcrops, plotted per area, with superposition of the diagenetic environments and the diagenetic events: (0) Sediment deposition, (1) early cementation, lithification and rock embrittlement, (2) fracturing, development of blocky spars and burial stylolites, and (3) tilting and development of tectonic stylolites. After Lamarche et al. (2012),[START_REF] Lavenu | Tectonic versus diagenetic origin of fractures in naturally fractured carbonate reservoir analogue: the Nerthe anticline (southeastern France)[END_REF][START_REF] Lavenu | Relating background fractures to diagenesis and rock physical properties in a platformeslope transect. Example of the Maiella Mountain (central Italy)[END_REF][START_REF] Lavenu | Background fractures in carbonates: inference on control of sedimentary facies, diagenesis and petrophysics on rock mechanical behavior. Example of the Murge Plateau (southern Italy)[END_REF].

Fig. 12 .

 12 Fig.12. Schematic of the evolution of carbonate rock through time. The middle bar compares the subsidence evolution to the diagenetic events. In platform carbonates, made of gradual facies transition, the early cementation is the key for early rock embrittlement and mechanical contrast fossilization, and the degree of cementation the main factor controlling fractures and stylolites density.

Acknowledgments

The authors want to thank Dr. Bertrand D. M. Gauthier for supporting the study, and TOTAL for funding the R&D "Carbonate" project, in which was included the first author's PhD project. Authors also gratefully acknowledged the laboratory technician from Aix-Marseille Universit e, Lionel Mari e, for the amount of work he has done for making plugs and thin sections. Dr. Christophe Matonti, Dr. Arnaud Gallois, Dr. Roland Salardon, Lisa Texier and Hanane Boukbida are thanked for their help on the field and in the lab. The authors also thanks Dr. Stephen Laubach for his valuable input and corrections on the manuscript as well as Dr. David Levitt, and to the two reviewers: Ole Petter Wennberg and Veerle Vandeginste.