Bryan Cochelin
email: bryan.cochelin.get@gmail.com

Charles Gumiaux

Dominique Chardon

Yoann Denèle

Benjamin Le Bayon

Multi-scale strainfield analysis using geostatistics: Investigating the rheological behavior of the hot Variscan crust of the Pyrenees (Axial Zone)

Keywords: Strain partitioning, variogram analysis, kriging interpolation, transpression, hot orogen

Multi-scale strainfield analysis using geostatistics :

investigating the rheological behavior of the hot Variscan crust of the Pyrenees (Axial Zone)

Introduction

Determination of principal strain and stress directions in the lithosphere has always been a key task for understanding and constraining plate tectonics. Mapping and analysis of strain/stress field patterns bring constraints on the mechanical behavior of crust/lithosphere. Such an approach allows evaluating (i) the meaning of spatial shear gradients and finite strain values in terms of displacement(s), (ii) the way deformation localizes, (iii) the role of rheological contrasts between adjacent lithologies during deformation and (iv) the rheological contrasts between structural levels and their degree of coupling. In active tectonic domains, approximate stress fields can be obtained using measured GPS velocity fields (e.g. [START_REF] Welsch | Finite element analysis of strain patterns from geodetic observations across a plate margin[END_REF][START_REF] Kahle | The GPS strain rate field in the Aegean Sea and western Anatolia[END_REF][START_REF] Grenerczy | Present crustal movement and strain distribution in Central Europe inferred from GPS measurements[END_REF], Zang et al., 2004;[START_REF] Titus | Geologic versus geodetic deformation adjacent to the San Andreas fault, central California[END_REF] or focal mechanisms from earthquakes (e.g. [START_REF] Mckenzie | Active Tectonics of the Mediterranean Region[END_REF][START_REF] Calais | Current strain regime in the Western Alps from continuous Global Positioning System measurements, 1996-2001[END_REF][START_REF] Engdahl | Global seismicity: 1900-1999[END_REF]. Structural analyses of fractures and faults further allow retrieving paleo-stress fields (e.g. [START_REF] Angelier | Tectonic analysis of fault slip data sets[END_REF][START_REF] Blès | Successive post-Variscan stress fields in the French Massif Central and its borders (Western European plate): comparison with geodynamic data[END_REF]. For ductile deformation, foliation and stretching lineations can be considered as good proxies of principal strain axes orientations [START_REF] Ramsay | Folding and Fracturing of Rocks[END_REF][START_REF] Siddans | Slaty cleavage -a review of research since 1815[END_REF][START_REF] Wood | Current Views of the Development of Slaty Cleavage[END_REF].

Foliation and lineation trajectories allow (i) investigating how the direction of principal strain axes change through space and (ii) deciphering homogeneous from heterogeneous strains by considering gradients and trajectory curvatures (i.e., convergence of trajectories) (see [START_REF] Cobbold | The geometrical significance of strain trajectory curvature[END_REF]. Following [START_REF] Cobbold | The geometrical significance of strain trajectory curvature[END_REF], the amount of strain is mostly difficult to quantify for natural case, therefore foliation and lineations trajectories -once coupled to a kinematic analysis of fabrics -are used to build paleo strainfields (or strain patterns). These trajectory maps are mostly extrapolated manually on the basis of structural measurements data (e.g. [START_REF] Brun | Strain patterns of pluton emplacement in a crust undergoing non-coaxial deformation, Sierra Morena, Southern Spain[END_REF][START_REF] Brun | L'origine des dômes gneissiques: modèles et tests[END_REF][START_REF] Chown | Tectonic evolution of the Northern Volcanic Zone, Abitibi belt, Quebec[END_REF][START_REF] Gautier | Crustal-scale geometry and kinematics of late-orogenic extension in the central Aegean (Cyclades and Ewia Island)[END_REF][START_REF] Bouchez | Two-stage deformation of the Mont-Louis-Andorra granite pluton (Variscan Pyrenees) inferred from magnetic susceptibility anisotropy[END_REF]. Moreover, extracting regional trends from individual local structural M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT 4 measurements may not be straightforward and would hardly be reproducible in case of low data density at a considered scale. The resulting trajectory maps may thus be intrinsically qualitative.

Considering the limitations of manual strainfields' interpolation, geostatistics can be used to quantitatively analyze spatial variations in the orientation of structural measurements and to obtain the best fitted interpolated maps, whatever the scale considered. The general approach has first been introduced by [START_REF] Matheron | Traité de Géostatistique Appliquée[END_REF]. Such geostatistical analysis has already been used for brittle deformation to retrieve stress fields using fracture or fault data set [START_REF] Escuder Viruete | 3-D stochastic modeling and simulation of fault zones in the Albalá granitic pluton, SW Iberian Variscan Massif[END_REF][START_REF] Rafiee | Application of geostatistical characteristics of rock mass fracture systems in 3D model generation[END_REF][START_REF] Koike | 3D geostatistical modeling of fracture system in a granitic massif to characterize hydraulic properties and fracture distribution[END_REF][START_REF] Hanke | Directional semivariogram analysis to identify and rank controls on the spatial variability of fracture networks[END_REF][START_REF] Riller | Fault-slip inversions: Their importance in terms of strain, heterogeneity, and kinematics of brittle deformation[END_REF] and 3D modeling of geological objects at local scale (e.g. folds, see [START_REF] Hillier | 3D form line construction by structural field interpolation (SFI) of geologic strike and dip observations[END_REF][START_REF] Grose | Structural data constraints for implicit modeling of folds[END_REF]. It has been also developed and successfully applied to the analysis of ductile fabrics directional data [START_REF] Gumiaux | Geostatistics applied to best-fit interpolation of orientation data[END_REF]. Based on a preliminary statistical analysis (the variogram analysis or variography), spatial interpolations can be computed, leading to automatic cleavage or foliation trajectory maps [START_REF] Gumiaux | Geostatistics applied to best-fit interpolation of orientation data[END_REF][START_REF] Gumiaux | Strain removal within the Hercynian Shear Belt of Central Brittany (western France): methodology and tectonic implications[END_REF].

In this study, we apply such geostatistical approach in order to obtain quantitative insights on the way strain has been spatially (i.e. horizontally and vertically) partitioned during convergence. Indeed, in orogenic contexts, geotherms control -at first order -vertical rheological layering and mechanical coupling in the lithosphere [START_REF] Davy | Modélisation thermo-mécanique de la collision continentale[END_REF]. Contrary to cold collisional belts such as the Alps and Taiwan [START_REF] Brun | Deformation of the continental lithosphere: Insights from brittle-ductile models[END_REF][START_REF] Lacombe | Basement-involved shortening and deep detachment tectonics in forelands of orogens: Insights from recent collision belts (Taiwan, Western Alps, Pyrenees)[END_REF], old accretionary or mature orogens have undergone critical weakening on a lithospheric scale, inducing specific vertical strain partitioning pattern and mass redistribution between the upper crust and a particularly weak lower crust [START_REF] Chardon | Flow of ultra-hot orogens: A view from the Precambrian, clues for the Phanerozoic[END_REF]. In order to investigate such space-time strain partitioning patterns, the present study focuses on a segment of the Variscides exposed in the Pyrenees, which underwent both thickening upper crust and the laterally flow of its mid-lower crust (Cochelin et , 2017). In a first step, processing of the overall structural dataset allows producing a strain trajectory map corresponding to the total finite strainfield of this crustal segment. Strainfield is then successively deconvoluted to constrain strain interference effects induced by crustal domains or units of contrasted mechanical behaviors (i.e. the upper crustal metasediments, the mid-lower crust and plutons). The analysis of strainfields based on strain trajectory modeling is performed from regional to local scale, with some focuses made on the Canigou gneiss dome and the Mt-Louis pluton.

Geological setting

The Pyrenees form an east-west trending intra-continental mountain range resulting from the collision between the Iberian and the European plates [START_REF] Choukroune | Tectonic evolution of the Pyrenees[END_REF]. During upper Cretaceous to Early Tertiary convergence, shortening was responsible for the exhumation and erosion of the Mesozoic cover, allowing exposure of the Paleozoic basement, which preserves Variscan (i.e., Late Paleozoic) orogenic deformation patterns. This portion of the Variscan orogen corresponds to the foreland of the belt and occupied the core of the Iberian-Armorican syntax during the latest stages of the orogeny (310-290 Ma) [START_REF] Cochelin | Vertical strain partitioning in hot Variscan crust: Syn-convergence escape of the Pyrenees in the Iberian-Armorican syntax[END_REF]. The largest uplifted unit of pre-Mesozoic rocks outcropping in the Pyrenees is called Axial Zone and forms the core and highest part of the Cenozoic belt (Fig. 1a). Differential denudation of the Axial Zone provides an access to a continuous vertical section of the Variscan crust from the shallowest parts in the west (i.e. Carboniferous sedimentary basins) to mid-crustal levels in the east. Despite the "recent" (i.e., Alpine) Pyrenean orogeny, the structural continuity of the Variscan crust is maintained in the Axial Zone; which forms a consistent tectonic domain affected by only local and moderate perturbations due to post-Paleozoic deformation [START_REF] Zwart | The variscan geology of the Pyrenees[END_REF][START_REF] Carreras | Tectonique Hercynienne[END_REF] [START_REF] Carreras | Framing the tectonic regime of the NE Iberian Variscan segment[END_REF][START_REF] Cochelin | Vertical strain partitioning in hot Variscan crust: Syn-convergence escape of the Pyrenees in the Iberian-Armorican syntax[END_REF][START_REF] Cochelin | Structural inheritance in the Central Pyrenees: the Variscan to Alpine tectonometamorphic evolution of the Axial Zone[END_REF]. The present study focuses on the pre-Mesozoic Axial Zone.

, M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT 6 1996,

Figure 1 here

The Axial Zone is mainly made of Ediacaran to Upper Carboniferous sedimentary rocks having undergone High Temperature-Low Pressure (HT/LP) metamorphism that induced regional partial melting (Guitard et al., 1996 and references therein). These metasedimentary rocks were intruded from ca. 310 to 290 Ma by large calc-alkaline plutons (Denèle et al., 2014 and reference therein).

Two contrasted structural levels are distinguished (de Sitter and Zwart, 1962;[START_REF] Zwart | The geology of the Central Pyrenees[END_REF][START_REF] Carreras | Structures and Tectonics at Different Lithospheric Levels Tectonic levels in the Palaeozoic basement of the Pyrenees: a review and a new interpretation[END_REF]. The low-grade metasedimentary rocks belonging to the upper crust are characterized by open to tight EW trending folds and steep N95-110°E striking axial plane cleavage (Figs. 1b &c) bearing steeply plunging stretching lineations [START_REF] Cochelin | Vertical strain partitioning in hot Variscan crust: Syn-convergence escape of the Pyrenees in the Iberian-Armorican syntax[END_REF][START_REF] Cochelin | Structural inheritance in the Central Pyrenees: the Variscan to Alpine tectonometamorphic evolution of the Axial Zone[END_REF]. The average strike of cleavage around N100°E associated to an average pitch of lineations around 83°W is compatible with transpression due to ~NS shortening [START_REF] Cochelin | Vertical strain partitioning in hot Variscan crust: Syn-convergence escape of the Pyrenees in the Iberian-Armorican syntax[END_REF]. The lower structural level comprises the partially molten mid-lower crust that appears in gneiss domes elongate parallel to the regional structural trend (Figs. 1a &c) and mostly bearing longitudinal stretching lineations [START_REF] Carreras | Tectonique Hercynienne[END_REF][START_REF] Mezger | Polymetamorphism and ductile deformation of staurolite-cordierite schist of the Bossòst dome: indication for Variscan extension in the Axial Zone of the central Pyrenees[END_REF][START_REF] Denèle | The Hospitalet gneiss dome (Pyrenees) revisited: lateral flow during Variscan transpression in the middle crust[END_REF][START_REF] Mezger | Transpressional tectonic setting during the main Variscan deformation: evidence from four structural levels in the Bossòst and Aston-Hospitalet mantled gneiss domes, central Axial Zone, M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT 37 Pyrenees[END_REF][START_REF] Cochelin | Vertical strain partitioning in hot Variscan crust: Syn-convergence escape of the Pyrenees in the Iberian-Armorican syntax[END_REF]. The chronological and structural relationship between steep and shallowly dipping structures has been subjected to debates since the 60's (e.g. de Sitter and [START_REF] De Sitter | Tectonic development in supra and infra-structures of a mountain chain[END_REF][START_REF] Zwart | The geology of the Central Pyrenees[END_REF][START_REF] Soula | A dynamic model of the structural evolution of the Hercynian Pyrenees[END_REF][START_REF] Carreras | Structures and Tectonics at Different Lithospheric Levels Tectonic levels in the Palaeozoic basement of the Pyrenees: a review and a new interpretation[END_REF][START_REF] Carreras | Tectonique Hercynienne[END_REF] but recent structural studies, especially within calc-alkaline plutons and gneiss domes, showed that they were formed during a single, protracted period of deformation, under dextral transpressional regime [START_REF] Bouchez | Two-stage deformation of the Mont-Louis-Andorra granite pluton (Variscan Pyrenees) inferred from magnetic susceptibility anisotropy[END_REF][START_REF] Gleizes | The main phase of the Hercynian orogeny in the Pyrenees is a dextral transpression[END_REF][START_REF] Mezger | Polymetamorphism and ductile deformation of staurolite-cordierite schist of the Bossòst dome: indication for Variscan extension in the Axial Zone of the central Pyrenees[END_REF][START_REF] Aurejac | Le complexe granitique de Quérigut (Pyrénées, France) ré-examiné par la technique de l'ASM : un pluton syntectonique de la transpression dextre hercynienne[END_REF][START_REF] Vilà | LP-HT metamorphism in a late orogenic transpressional setting, Albera Massif, NE Iberia: implications for the geodynamic evolution of the Variscan Pyrenees[END_REF][START_REF] Denèle | The Hospitalet gneiss dome (Pyrenees) revisited: lateral flow during Variscan transpression in the middle crust[END_REF]2008;[START_REF] Mezger | Transpressional tectonic setting during the main Variscan deformation: evidence from four structural levels in the Bossòst and Aston-Hospitalet mantled gneiss domes, central Axial Zone, M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT 37 Pyrenees[END_REF][START_REF] Mezger | Transpressional tectonic setting during the main Variscan deformation: evidence from four structural levels in the Bossòst and Aston-Hospitalet mantled gneiss domes, central Axial Zone, M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT 37 Pyrenees[END_REF]; Cochelin et al., 2017 among others). In spite of the N-S convergent context, HT-rock of the midlower crust were exhumed in gneiss domes [START_REF] Mezger | Polymetamorphism and ductile deformation of staurolite-cordierite schist of the Bossòst dome: indication for Variscan extension in the Axial Zone of the central Pyrenees[END_REF][START_REF] Vilà | LP-HT metamorphism in a late orogenic transpressional setting, Albera Massif, NE Iberia: implications for the geodynamic evolution of the Variscan Pyrenees[END_REF][START_REF] Aguilar | P-T-t-d evolution of orogenic middle crust of the Roc de Frausa Massif (Eastern Pyrenees): a result of horizontal crustal flow and Carboniferous doming?[END_REF]. This has been favored by strain localization near the biotite-andalusite transition, forming retrogressive extensional shear zones (Van den Eeckout, 1986;Van den Eeckout and Zwart, 1988;Mezger and Paschier, 2003;[START_REF] Cochelin | Vertical strain partitioning in hot Variscan crust: Syn-convergence escape of the Pyrenees in the Iberian-Armorican syntax[END_REF] while the core of the domes was affected by horizontal longitudinal flow [START_REF] Denèle | The Hospitalet gneiss dome (Pyrenees) revisited: lateral flow during Variscan transpression in the middle crust[END_REF][START_REF] Mezger | Transpressional tectonic setting during the main Variscan deformation: evidence from four structural levels in the Bossòst and Aston-Hospitalet mantled gneiss domes, central Axial Zone, M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT 37 Pyrenees[END_REF][START_REF] Mezger | Transpressional tectonic setting during the main Variscan deformation: evidence from four structural levels in the Bossòst and Aston-Hospitalet mantled gneiss domes, central Axial Zone, M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT 37 Pyrenees[END_REF][START_REF] Cochelin | Vertical strain partitioning in hot Variscan crust: Syn-convergence escape of the Pyrenees in the Iberian-Armorican syntax[END_REF].

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT 7
Figure 2 here

The transition between the thickening upper crust and the laterally flowing lower crust is well exposed in the Central and Eastern Pyrenees by the succession of ~EW trending regional scale synforms and gneiss domes (e.g. the Villefranche syncline and the Canigou gneiss dome, see

Figs. 1c & 2). The Canigou gneiss dome is characterized by a pervasive foliation reworking

Ordovician plutons turned into orthogneisses and their country-rocks (Fig. 2, [START_REF] Casas | New insights into the Late Ordovician magmatism in the Eastern Pyrenees: U-Pb SHRIMP zircon data from the Canigó massif[END_REF]Martinez et al., 2011). Stretching lineations are mainly NE trending in the core of the dome and are associated with dominant top-to-the SW shearing [START_REF] Guitard | Linéations, schistosité et phases de plissement durant l'orogenèse hercynienne dans les terrains anciens des Pyrénées orientales, leurs relations avec le métamorphisme et la granitisation [with discussion[END_REF][START_REF] Soliva | Shear deformation and emplacement of the gneissic Canigou thrust nappe (Eastern Pyrenees)[END_REF][START_REF] Carreras | Tectonique Hercynienne[END_REF][START_REF] Cochelin | Vertical strain partitioning in hot Variscan crust: Syn-convergence escape of the Pyrenees in the Iberian-Armorican syntax[END_REF]. This massif was lately cut by post-Paleozoic faults like the Alpine Canigou Thrust and the Neogene Têt-Cerdagne and Py-Fillols faults (Fig. 2).

Immediately to the north west of the Canigou gneiss dome crops the Mt-Louis pluton, which was emplaced at ca. 301-303 Ma (Figs. 1a & 2) [START_REF] Denèle | Timing of granite emplacement, crustal flow and gneiss dome formation in the Variscan segment of the Pyrenees[END_REF]. It shows a complex internal structure and can be divided into three subintrusions [START_REF] Bouchez | Two-stage deformation of the Mont-Louis-Andorra granite pluton (Variscan Pyrenees) inferred from magnetic susceptibility anisotropy[END_REF] and, like most plutons in the Axial Zone, this pluton was deformed under sub-magmatic to solid-state conditions [START_REF] Bouchez | Two-stage deformation of the Mont-Louis-Andorra granite pluton (Variscan Pyrenees) inferred from magnetic susceptibility anisotropy[END_REF][START_REF] Debon | Magmatisme hercynien[END_REF][START_REF] Denèle | Timing of granite emplacement, crustal flow and gneiss dome formation in the Variscan segment of the Pyrenees[END_REF]

Methodology and data set

The computation of cleavage trajectory maps results from an interpolation process of the striking value between individual structural measurements. This technique allows imaging the continuous spatial variations of the finite strain field that can be deduced from discrete natural observations. Thereby, some difficulties may arise when measurements (i) are heterogeneously distributed through the area and/or (ii) display significant local-scale variations in orientation. In particular, it is often not straightforward when one has to a priori choose whether to highlight local tendencies or regional ones through a hand drawing process, using curvy or smooth trajectories shape, respectively. Here, we employ a geostatistical approach to automatically map both cleavage and lineation trajectories and cleavage dip angle variations from a given field structural dataset.

Ideally, interpolations should be made from foliation-lineation pairs in order to integrate the full information of the finite strain ellipsoid's orientation (see the critical review of [START_REF] Davis | Modern methods of analysis for three-dimensional orientational data[END_REF]. However, such approach requires a very homogeneous dataset in which each lineation trend/plunge measurement can be associated to a foliation strike/dip one and vice versa. In the case of regional-scale studies (as the one exposed here), structural data are often compiled from various types and ages (e.g. digitized maps from PhD theses, articles), and such a prerequisite is unfortunately hardly tenable. Alternatively, structural data sets can be decomposed into four onedimensional numbers (i.e. strike, dip, trend and plunge) and analyzed separately. While dip value is interpolated using simple kriging on a scalar variable (for a 0-90° range in that case), kriging interpolation of the strike value (which is of circular data type) is performed following the approach developed by [START_REF] Gumiaux | Geostatistics applied to best-fit interpolation of orientation data[END_REF]. Based on a theoretical and unscaled case study and, in particular, some keys for the pre-interpolation variogram analysis (3.2). This synthetic cleavage and stretching lineation dataset simulates interactions between a roughly E-W regional trending structure and termination of an elliptical body, to the west, which may correspond to an intrusion. It results in a triple point as visible in the center of the map (Fig. 3.1).

Computation of automatic trajectory maps

The computation process has been developed in a GIS environment and all the structural measurements are firstly compiled into a spatial database (37 cleavage orientation measurements in the synthetic case study, with 20 stretching lineations; Fig. 3.1). Secondly, spatial variations of the dataset values can be quantified by the computation of experimental variograms (Figs 3.2 & 4.2) [START_REF] Matheron | Traité de Géostatistique Appliquée[END_REF]. Initially developed for mining purposes and applied to scalar values [START_REF] Matheron | Traité de Géostatistique Appliquée[END_REF], this method has been adapted by [START_REF] Gumiaux | Geostatistics applied to best-fit interpolation of orientation data[END_REF] for the analysis of circular angular variables such as structural striking data. The goal is to assess the degree of correlation between values of the studied variable (e.g. azimuth or dip of cleavage data) measured at different locations and according to the interdistance separating them. If this correlation degree changes with increasing interdistance, the variable is thus categorized as regionalized. The variogram analysis allows (i) quantifying the maximum distance to which distant data values show significant mutual dependence, (ii) constraining their dependence function to the distance, (iii) deciphering any potential anisotropy in the spatial variations of the regionalized variable and (iv) given the scale of the study area, filtering "regional" and "local" variation trends from the dataset (see details in the next section). In a third stage, kriging interpolation can be achieved on both the dip values and azimuth directions following the method developed by [START_REF] Gumiaux | Geostatistics applied to best-fit interpolation of orientation data[END_REF] proportional nugget effect and sill values; see part 3.2 below). In the example, interpolations allow obtaining cleavage direction map, cleavage dip map (Fig. 3.3) or lineation trend maps (Fig. 4.3), as continuous grids covering the overall study area, and where corresponding data is available. During the two-fold computation/interpolation of the directional map, the local direction variance can be estimated using the parameter R, which reflects the degree of dispersion of data [START_REF] Upton | Spatial Data Analysis by Example[END_REF]. In [0,1], low R values reflect high variation in measured direction values for short distances, whereas a R of a theoretical maximum value of 1 reflects strictly parallel cleavage directions in the neighborhood (see Gumiaux et al., 2003 and Figs 3.3c & 4.3b). In the example, low R values computed from the cleavage data point to significant directional variations induced by the eastern termination of the elliptical body and its associated cleavage triple point (forming a singular point; see Fig. 3.3c). In contrast, the trend of stretching lineations shows only minor variations for short distances with R values above 0.5 for most of the study area (Fig. 4.3b). Trajectories are then computed, giving a good picture of variation of cleavage directions at the scale of the studied area and reproducing at high resolution cleavage singular points (including neutral or triple points), both within and outside the simulated intrusion, and small scale variations within the elliptical geological body (Fig. 3.4). Besides, the continuous and radial distribution of stretching lineations, as implied by dataset in Figure 3.1, is also well-modeled and imaged following this fully automatic process (Fig. 4.4). In order to better constrain the strain field geometry, trajectory maps can be superposed over the interpolated dip angle map. In the example case study, shallow dipping cleavage areas correspond to the elliptical shape structure to the west and to more scattered cleavage direction zones, as highlighted by rather high R values (compare maps in Figs. 3.4 & 4.4). Note that the low R values obtained in

Variogram analysis

As a general principle, two close data measurements may have more similar values than two distant ones. Yet, superimposed phenomenon at different scales can control variations in the regionalized variable. For instance, cleavage orientation can result from the combination of (i) a tectonic regional-scale maximum shortening direction, (ii) influence of the expansion of a neighboring synkinematic intrusion (medium-scale) and (iii) very local cleavage refraction due to lithological contrast in the rock unit considered. Thus, in order to determine and evaluate the mutual dependence between data values, experimental variograms are calculated as a squared root difference between individual angle values as a function of their separating distance (Figs.

& 4.2); increasing values with distance show decreasing correlation between data values with

a rate depending on the steepness. Thus, the shape of the resulting experimental variogram plot gives key information on the spatial structure of the strain field. This is estimated by fitting a theoretical function to the experimental variogram (see continuous grey curves in Figs. 3.2 and 4.2) and extracting several parameter values from the latter, as detailed below (see Matheron,

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

12 -An increasing variogram can reach a constant sill value for a maximum correlation range distance value (Fig. 3.2 left) and is then defined as stationary. This shows that the kinematical and/or rheological phenomena controlling the spatial strain variations have a distance of influence shorter than the maximum width of the study area. In the synthetic example, the stationary cleavage directions variogram highlights deformation process with effects over 9 to 10 distance units (Fig. 3.2 left). This feature is also consistent with the computed trajectory map which displays the overall structure of the junction area (Fig. 3.4). On the contrary, the non-stationary variograms computed for the lineation directions or cleavage dip angle show that the phenomenon driving strain orientation's variations has some influence out of scale of the study area (i.e. range is not reached, Figs.

3.2 right & 4.2),
-Adjustment and analysis of the trend of the variogram for shorter distances is of primary importance and in particular its initial slope characterizing the type of influence linked to the phenomenon. Exponential function (Fig. 3.2 left), for instance, displays a notably quick increase from its original. This corresponds to continuous spatial variations in the strain direction but with a rather limited size (or rapid decrease) of influence. In the synthetic example, such feature is attested to by the narrowness of triple junction's zones and rather high curvature of the trajectories there (Fig. 3.4). In contrast, the Gaussian theoretical variogram function displays a nearly null slope from its origin (Fig. 4.2) showing a very slow decrease of influence of the phenomenon in space. This is expressed in the strain field by low curvature and very smooth variations in the trajectories pattern (Fig. 4.4).

-Finally, a non-null ordinate value at the origin of the variogram -the nugget-effectreflects some variations in the strain field that cannot be considered regarding the given

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT
13 density of the data points. This includes a proportion of either background noise in the dataset (i.e. different kinds of uncertainties or errors made during the acquisition process), very local rheological effects (such as cleavage refraction in layered rocks for instance) that can be regarded as insignificant at the scale of the strain field considered, or both. In the synthetic case study presented here, the variogram computed for stretching lineations direction measurements displays a null nugget effect (original variogram Γ is zero, Fig.

4.2)

. This corresponds to an idealistic situation where neither measurement error nor very local deviation effects exist in the dataset and 100% of the cleavage direction variations can be modeled in the trajectories map computation. For natural case studies, a significant proportion of the variance may be included in the nugget and, thus, not modeled.

In order to detect potential anisotropy in strainfield variations, directional individual variograms can be calculated by selecting data pairs along predefined space direction classes. This technique is particularly efficient (i) to quantify the range of significant variations of values in the investigated directions and (ii) to highlight main structural anisotropies such as folds or sigmoidal trends induced by regional shearing [START_REF] Gumiaux | Geostatistics applied to best-fit interpolation of orientation data[END_REF], 2004[START_REF] Branquet | Synkinematic mafic/ultramafic sheeted intrusions: Emplacement mechanism and strain restoration of the Permian Huangshan Ni-Cu ore belt (Eastern Tianshan, NW China)[END_REF].

Data set

All the structural data available to date in the Axial Zone of the Pyrenees [START_REF] Cochelin | Vertical strain partitioning in hot Variscan crust: Syn-convergence escape of the Pyrenees in the Iberian-Armorican syntax[END_REF] have been treated. This extensive structural compilation (see supplementary materials) comprises 14 275 cleavage and foliation plane measurements, including 12 665 with dip information and only 1 623 mineral-stretching lineations. At first, analysis and interpolation of structures' orientation is performed for the entire dataset over the Axial Zone. Then, in order to deconvolute the finite strain field, three distinct structural sub-sets have been selected based on the type of geological units in which structural measurements have been made: (i) a first set for the upper

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT
14 crustal metasedimentary rocks called hereafter 'upper crust' comprising 5 366 cleavage planes and 272 stretching lineation measurements, (ii) a second set for the plutons comprising 3610 foliation planes (with a zoom on the Mt Louis pluton) and (iii) a third set for mid-lower crustal rocks called hereafter 'mid-lower crust' including 5 299 foliation planes and 1 351 mineralstretching lineation measurements (with a zoom on the Canigou gneiss dome). Spatial selection have been made using digitalized outlines of the corresponding units in the maps.

Cleavage/foliation data points are spatially duly distributed over the Axial Zone (see supplementary materials) and the strikes show a normal distribution centered on N96.5°E and a skewness factor of about 0 (nonparametric skew used here; Fig. 5a). Cleavage strikes from the mid-lower crust follow a very similar normal distribution (Fig. 5b) with a mean value at N94.6°E.

In contrast, the directional data from the upper crust are skewed to the left and show a mean value of N102°E (Fig. 5c). For the whole Axial Zone, the dip values of cleavage/foliation data show a mean value of 51° and a skewness around 0 (Fig. 5d). In the middle crust, dip values also show a skewness about 0 with a mean value of 47° (Fig. 5e) while data in the upper crust are skewed to the steep values and show a mean value of 51° (Fig. 5f). In the whole data set, flat-lying foliations (i.e. < 20°) represent less than 7% of the measurements. Therefore, their relative impact on interpolations of strike values can be considered as negligible in the present case study.

Stretching lineations have principally been measured in gneiss domes. Considering their data points' spatial distribution, geostatistics were applied to mid-lower crust and solely statistics for upper crust data.

Results

Axial Zone strain field

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT 15
The omni-directional variogram for the Axial Zone dataset displays a sill reached for a 10km general range value (grey dots/black curve, Fig. 6a). An abrupt change in slope at 1.8km indicates a second, intermediate range, which indicates the change to a different structuration, from local to more regional scale phenomenon effects. As a whole, this variogram can be modeled with (i) a nested (or imbricated) structure of two spherical functions of distinct range values and (ii) a global sill of 1095 square degrees and a nugget-effect value of 503 square degrees (corresponding to 46% of the total variance for this dataset).

The omni-directional variogram obtained for dip measurements over the Axial Zone shows an exponential trend, reaching a sill of about 480 square degrees at a 21 km range and showing 60% of proportional nugget effect (grey dots/black curve, Fig. 6b). The later suggests, on average, a significant variability of the dip values measured over short distances (i.e. distances shorter than 1 km).

Figure 6 here Based on the variograms (Fig. 6a), interpolation of cleavage trajectories at the scale of the whole Axial Zone was performed using a 10km range. In this case, the interpolated trajectories display an average N100°E parallel trend over most of the Axial Zone and rotates to an average of N130°E in eastern part of the Axial Zone (Fig. 7a). Based on the R parameter map (Fig. 7b), domains belonging to the upper crust show limited cleavage strike dispersions (with R values above 0.5). This is also shown by the regular trend of the trajectory map in the upper-crust domains and especially in a ~30km wide N-S corridor displaying sub-parallel trajectories in Central Pyrenees (Fig. 7b). This corridor is devoid of plutons (see Figs. corridor of sub-parallel trajectories referred to above (Fig. 7b). This variation corresponds to the fan of cleavage that had been early noticed by [START_REF] Zwart | The geology of the Central Pyrenees[END_REF] in this region. Calculation of the gradients between interpolated dip values [START_REF] Carrara | GIS techniques and statistical models in evaluating landslide hazard[END_REF] allows identifying zones or bands with rough variations of dip, and, in another hand, undisturbed domains such as in the upper crust level (Fig. 8b; note in particular the delineation of the ~NS corridor in the central part).

At second order, the strain field appears rather heterogeneous. Gneiss domes and plutons are very well individualized and evidenced by circular to concentric trajectories (Fig. 7a). Heterogeneity is also highlighted by converging and tightening of the trajectories underlining at least 5 longitudinal 100 to 150 km long shear zones. They mostly trend ~N100°E and locally N130°E, where they wrap around gneiss domes and plutons. Trajectories that converge in N100°E trending bands show rather symmetrical pattern with both clockwise and counterclockwise rotations on both sides of these bands. In contrast, N130°E trending bands induce convergence and bending of trajectories as smooth sigmoids. Such sigmoidal patterns (i.e. asymmetrical) for trajectories are compatible with dextral shearing (Fig. 7a) [START_REF] Ramsay | Strain variation in shear belts[END_REF]. These sigmoidal patterns form locally asymmetrical singular points (neutral or triple points) at western and eastern tips of elliptical bodies (e.g. around the Néouvielle pluton), which is also compatible with dextral shearing. Otherwise, cleavage singular points at eastern and western terminations of the domes and plutons appear mainly symmetrical in respect of the elongation axe of these domes and plutons. showing internal imbricated magmatic structures (e.g. Cauterets plutons, Fig. 7b). In the eastern part of the Axial Zone, high dispersion zones (i.e. with low R values) surround the gneiss domes. Some folded gneiss domes display internal dispersion too (e.g. the Aston gneiss dome, see [START_REF] Denèle | Decoupling between the middle and upper crust during transpression-related lateral flow: Variscan evolution of the Aston gneiss dome (Pyrenees, France)[END_REF]). Discrete dispersion bands, sometimes oblique to the cleavage trajectories, are also mapped within gneiss domes, like in the Canigou (Fig. 7b). This kind of dispersion could be related to post-Variscan faulting (see section 4.4).

The impact of gneiss domes and plutons on the regional strainfield is also well shown by cleavage dip variations (Fig. 8b) outlining concentric shapes and reflecting short-length gradients between the cores and envelops of these bodies (e.g. the Cauterets plutons). Other shortwavelength variations that are not directly related to domes or plutons display highly to slightly oblique trends to trajectories. For instance, N60° to N80°E trending bands are detected south of the Canigou gneiss dome. This reflects late tilting of the cleavage along Neogene normal faults as well as Alpine thrusts (Fig. 8b). The relative scarcity of such perturbations shows that the Variscan structures are only weakly affected by post-Paleozoic deformation phases within the Axial Zone (see [START_REF] Cochelin | Vertical strain partitioning in hot Variscan crust: Syn-convergence escape of the Pyrenees in the Iberian-Armorican syntax[END_REF][START_REF] Cochelin | Structural inheritance in the Central Pyrenees: the Variscan to Alpine tectonometamorphic evolution of the Axial Zone[END_REF].

As a whole, the overall strain field looks homogenous on a regional-scale and locally disturbed by the presence of gneiss domes and plutons, forming areas with curved cleavage trajectories and lower dip values (e.g. in eastern Axial Zone).

Upper crust

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT 18
The omni-directional variogram computed for the upper crust dataset (blue dots/curve, Fig. 6a)

shows an exponential trend with a significantly lower variance (849 square degrees) than that of the entire data set, but a similar nugget effect proportion (47%). The best-fitted theoretical variogram suggests an intermediate range at ca. 5 km. The four directional-variograms display a same trend than the omni-directional one (compare Figs. 6a & 9a,b) with similar range (5km), sill (810 square degrees) and proportional nugget effect (between 47 and 53%) values. The NE-SW (i.e. from N22.5° to N67.5°E) trending directional-variogram differs from the others, with a sinusoidal trend of periodicity around 18 km. This illustrates systematic regional scale undulations of directional data, along the NE-SW direction (i.e. forming perpendicular ~NW-SE striking bands like).

Figure 9 here Considering cleavage dip values in the upper crust, the omni-directional variogram is similar to that of the entire Axial Zone dataset but tends to show a greater instability (blue dots/curve, Fig. 6b). Directional-variograms are calculated from dip measurements parallel and perpendicular to the average cleavage strike (i.e. ~N100°E, Fig. 9c). The "parallel" N100°E directional-variogram has an exponential trend with a 5km range (much shorter than for the omni-directional variogram) and a proportional nugget effect around 52%. In contrast, the "perpendicular" N10°E oriented directional variogram (Fig. 9d) displays an exponential trend but is non-stationary (it does not reach a sill up to the maximum possible 40 km computational distance). The latter depends on the N-S width of the Axial Zone and on the critical number of data pairs available for calculation, which dramatically decreases beyond 40km and makes the variogram unstable (Fig. 9d). These two directional-variogram show that cleavage dip value does not have regional-scale variations along strike in the upper crust, whereas such variations are evidenced in the N-S Interpolations of cleavage directions have been performed up to 20 km interdistance, which corresponds to the largest range revealed by the variogram analysis of the upper-crust (i.e. in directional variogram, Fig. 9a). Regional trajectories computed from the sole upper crust are roughly parallel. They display an overall arc shape (Fig. 10a) from N110°E trends in the western part of the Axial Zone, longitudinal trends in the central Axial Zone to N130-140°E trends in the easternmost Axial Zone, which is quite similar to trajectories obtained in Fig. 7a. At second order, E-W elongated lens patterns can be deciphered in the Eastern Axial Zone. In the western Axial Zone, trajectories are weakly disturbed and only display a few variations over short (<10 km) wavelength undulations. The convergence of N110°E trending trajectories in the median part of the Axial Zone (e.g. between Ct and Nv) and disappearing to the east is the only apparent regional strain gradient revealed by the interpolation (Fig. 10a). Indeed, the rather consistent spacing of the trajectories indicates a more homogeneous strain pattern than obtained using the entire data set (Fig. 7a). Similarly as on the first trajectory map, singular points are symmetrically displayed at western and eastern tips of the 10 to 20 km elliptical bodies and aligned along the cleavage mean trend (Fig. 10a). The only asymmetrical converging pattern (i.e. sigmoidal) of trajectories are located, from west to east, around the Chiroulet-Lesponne domes (CL), the Néouvielle plutons (Nv) and the eastern part of the Canigou dome (Ca) (Fig. 10a). These sigmoidal trajectories again suggest dextral shearing within shear bands. The interpolation of cleavage dip values for the upper crust data set upholds a NS segmentation over the Axial Zone (Fig. 10b) with a northern part showing steep to vertical cleavage and a southern one where the cleavage is progressively tilted southward to reach moderate to low northward dips. In the Eastern Axial Zone, this segmentation of cleavage dip is still visible but is disturbed by the influence of gneiss domes and plutons (compare Figs. 1a & 10b). In central and western Axial Zone, the influence of Bossost and Chiroulet-Lesponne domes is also documented by two "anomalies" of shallowly-dipping cleavage inside the northern steeply-dipping domain (Fig. 10b). Even though only the upper crust structural dataset is considered in the present interpolations, both cleavage trajectories (showing lens shapes map structures) and interpolated dip map display anomalies surrounding gneiss domes and plutons (compare Figs. 7 & 10).

Whatever the cleavage dip, stretching lineations measured in the upper crust are, in most cases, characterized by high pitch values (>60° for 61 % of the set). As cleavage is mostly steeply dipping (Fig. 1b & 10b), lineations often have high plunges. In regional shear zones (Figs. 1 & 7a), lineations are associated with either dip-slip or reverse sense of shear with minor dextral component [START_REF] Cochelin | Vertical strain partitioning in hot Variscan crust: Syn-convergence escape of the Pyrenees in the Iberian-Armorican syntax[END_REF][START_REF] Cochelin | Structural inheritance in the Central Pyrenees: the Variscan to Alpine tectonometamorphic evolution of the Axial Zone[END_REF]. If plotted in a Dip-Pitch-Plunge ternary diagram (Fig. 11-left; see [START_REF] Chardon | Flow of ultra-hot orogens: A view from the Precambrian, clues for the Phanerozoic[END_REF], two distinct sub-populations are identified. Most of the dataset shows a trend from the vertical flow to the horizontal flow pole, which reflects a decrease of foliation dip with nearly constant high pitch values. The remaining part of the data is closer to the Orogen-parallel "strike-slip" flow pole, reflecting lower pitches on steep planar fabrics (Fig. 11left). At very first order, the upper-crust is thus characterized by (orogen-normal) horizontal shortening and (sub-) vertical stretching. In more details, variations from the vertical to the horizontal flow poles can be related to (i) the regional-scale N-S segmentation of the upper-crust strain in two main cleavage dip domains (Fig. 10b) or (ii) the shallowing of upper crust cleavage

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

21 around the domes (e.g. the Bossost dome and the Canigou dome, Fig. 10b). Furthermore, at a third order, the dispersion of the data between the vertical flow pole and the orogen parallel "strike-slip" flow pole (left axis of the triangle plot) reflects strain partitioning within strike-slip dominated transcurrent shear zones. In the Axial Zone, these shear zones -mainly N120-130°E trending -show dextral-reverse kinematics and are flanking gneiss domes and plutons [START_REF] Cochelin | Vertical strain partitioning in hot Variscan crust: Syn-convergence escape of the Pyrenees in the Iberian-Armorican syntax[END_REF]; e.g. the Mérens and Têt shear zones; Fig. 7a).

Figure 11 here

Upper crustal plutons: focus on the Mt Louis massif

The omni-directional variogram computed from the fabric direction measurements in the Mt Louis pluton (extracted from Anisotropy of Magnetic Susceptibility -AMS -foliation orientation extracted for each station of [START_REF] Bouchez | Two-stage deformation of the Mont-Louis-Andorra granite pluton (Variscan Pyrenees) inferred from magnetic susceptibility anisotropy[END_REF] is stationary (Fig. 12a): it reaches a 1142 square degrees sill at 3.4 km range and includes 54% of proportional nugget effect. The global cleavage direction variance for this pluton (corresponding to the sill) and its variability at very short distances (the nugget) are similar to those obtained for the global omni-directional variogram computed for the entire Axial Zone (Fig. 6a). This suggests that this massif -and certainly plutons in general -are of prime importance in the structuration of the strainfield. Yet, the Mt Louis experimental variogram displays a distinct range value and follows a Gaussian trend (with flat behavior close to the original) corresponding to particularly smooth tendency variations of the cleavage directions at the scale of the pluton. The regional interpolation (Figs. 7a & 12b) over the Mt Louis pluton and its wall rocks is characterized by longitudinal trajectories on the central-western part of the pluton with large and smooth undulations, whereas the concentric trajectories inside the pluton is only and roughly evidenced in its eastern part. As expected, the resulting interpolated trajectories at local-scale are more sensitive to local variations than those

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

22 obtained with larger range parameters (compare Fig. 12b &c). The sub-massifs are well individualized by concentric trajectories patterns, as exemplified in its western and eastern parts.

These concentric patterns are in agreement with the trajectories obtained manually by Bouchez and Gleizes (1995) (Fig. 12d). Furthermore, our interpolations looks delineating three NW-SE bands showing asymmetrical converging trajectories (s-shape sigmoidal structures) coinciding with dextral shear zones inferred by [START_REF] Bouchez | Two-stage deformation of the Mont-Louis-Andorra granite pluton (Variscan Pyrenees) inferred from magnetic susceptibility anisotropy[END_REF] (Fig. 12d). Even though interpolated sigmoids are compatible with apparent dextral shearing (Fig. 12c), trajectories display very smooth variations and non-localized/pervasive shearing.

Figure 12 here

The lower crust and focus on the Canigou dome

The omni-directional variogram of the mid-lower crust exposed in the domes displays, as a whole, a slightly higher variance (i.e. average sill) and higher proportional nugget effect values (50%) than for the entire Axial Zone dataset (orange dots/red curve Fig. 6a). It shows a nested structure with (i) an exponential evolution up to 6-7 km and, beyond, (ii) a sinus trend of 15-20 km periodicity around the average sill (Fig. 6a). Such a rhythmic variogram can be related to the peculiar spatial distribution of measurements, which form evenly spaced clusters centered on the gneiss domes. The periodicity may thus simply reflect the average size of these domes, which is known as a hole-effect.

Focusing on the Canigou dome, the local-scale omni-directional variogram computed for cleavage striking data is stationary. It displays a lower maximum range (5 km) and a lower proportional nugget value (~38%) than the regional-scale one (compare Figs. 6a & 13a). This indicates less variable cleavage orientations and thus a more consistent strain field when focusing

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

23 on the lower-crust. The variogram follows an exponential trend corresponding to a rapid decreasing effect of the phenomena with distance. The interpolated foliation trajectories reveal two half domes (Fig. 13b), with E-W hinges. Considering the R parameter (Fig. 13c), the contact zone between the two sub-domes coincides with a N50°E to N80°E trending band of highly dispersed cleavage directions. The high-R band actually coincides with the Neogene Py-Fillols normal faults, which cross-cuts the gneiss dome. Besides, the low-R domain at the south-eastern boundary of the dome also appears to be related to post-Variscan faults and reflects the local impact of the South Canigou Alpine thrust (Fig. 2) as suspected from the regional-scale map (Fig. 7b).

Figure 13 here The stretching lineations measured within (Fig. 14a) the Canigou dome display only limited variations in direction. The omni-directional variogram (Fig. 14b) is characterized by a high ~57% proportional nugget effect and is defined by a sinusoidal trend. The interpolated trajectories of stretching lineations exhibit an overall NE-SW trend (in red, Fig. 14c), with a very smooth pattern. In contrast to the cleavage measurements, lineations do not show disturbance that could be related to post-Variscan faulting. Indeed, the Py-Fillols faults are sub-parallel to the lineation trajectories and the latter thus look poorly affected by late brittle deformation. First order stretching trajectories are oblique on the hinge of the dome (compare Figs. 13b &14c). In greater details, trajectories display a systematic and progressive curvature, stretching becoming progressively more transverse along the dome's limbs than in its core. This defines a map radial trend from NE-SW in the lowermost structural levels of the dome to a NNE-SSW direction in the southern limb and up to NW-SE trend in its north-western part. This latter trend is not well expressed because the northern flank of the dome is actually truncated by the Neogene Têt-

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT 24
Cerdagne fault (Fig. 2). As a whole, interpolated cleavage and lineation trajectories within the Canigou massif thus look particularly smooth and continuous; cleavage is only disturbed by significant but quite local effect of recent faulting.

Figure 14 here Considering the overall Axial Zone, the average stretching direction is c.a. N98°E in the midlower crust (Fig. 15a) and, in a Dip-Pitch-Plunge diagram, data are clustered close to the orogenparallel horizontal flow pole (Fig. 11-right). Dispersion of lineations towards the orogen-parallel "strike-slip" flow pole (Fig. 11-right) must be due to the late amplification of gneiss dome by folding, producing steep to overturned southern limbs in some cases (see [START_REF] Denèle | The Hospitalet gneiss dome (Pyrenees) revisited: lateral flow during Variscan transpression in the middle crust[END_REF][START_REF] Mezger | Transpressional tectonic setting during the main Variscan deformation: evidence from four structural levels in the Bossòst and Aston-Hospitalet mantled gneiss domes, central Axial Zone, M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT 37 Pyrenees[END_REF][START_REF] Vilà | LP-HT metamorphism in a late orogenic transpressional setting, Albera Massif, NE Iberia: implications for the geodynamic evolution of the Variscan Pyrenees[END_REF][START_REF] Cochelin | Vertical strain partitioning in hot Variscan crust: Syn-convergence escape of the Pyrenees in the Iberian-Armorican syntax[END_REF]. In greater details, stretching lineations display two distinct direction families with regard to the gneiss domes' cartographic long axes. Longitudinal to slightly oblique stretching/flow is dominant, representing ca. 90% of the measurements within the domes (see the above described focus on the Canigou). Transverse stretching directions (i.e.

NS trends) are restricted to the domes' envelope and represent only 10% of the exposed middle crust (Figs. 15a & 15b). In details, oblique WNW-ESE to NW-SE trends are documented in the core of all gneiss domes from the northern and western parts of the Axial Zone (i.e. Aston, Lys-Caillaouas, Bossost, Chiroulet and Lesponne domes) and from the Albères massif to the East (Fig. 15b). Purely longitudinal stretching is restricted to the Hospitalet dome while the Canigou shows a specific ENE-WSW trend (Figs 15b). Transverse stretching directions are observed in the envelope of the Chiroulet, Aston and Canigou domes and can be noticed in the Bossost dome (Fig. 15b). In the Chiroulet and Aston, the transition between longitudinal and transverse direction of stretching is abrupt, whereas it appears to be more gradual from the core to the limbs of the Canigou, Hospitalet and the Bossost domes (Fig. 15b). We consider the direction of

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT 25
stretching in the core of gneiss domes as representative of the deep flow within the crust while direction of stretching in domes' envelope is related to their exhumation [START_REF] Cochelin | Vertical strain partitioning in hot Variscan crust: Syn-convergence escape of the Pyrenees in the Iberian-Armorican syntax[END_REF].

Discussion

Strain partitioning and crustal flow in the Variscan crust of the Pyrenees

The variogram calculated for the upper crust is stationary and shows the minimum of all the global variances (i.e. sill values) calculated in this study (compare variograms in Fig. 6). Such features reflect an overall steady strain direction throughout this domain, which bears the regional structural grain (Fig. 7a). Regional-scale interpolated cleavage trajectories for the upper crust are longitudinal and roughly collinear (Fig. 10a), which is typical of homogeneous deformation [START_REF] Choukroune | Soft lithosphere during periods of Archaean crustal growth or crustal reworking[END_REF]. Such cleavage trajectories pattern coupled with a dominant steep stretching (Fig. 11-left) is compatible with a bulk pure shear dominated deformation regime [START_REF] Tikoff | Stretching lineations in transpressional shear zones: an example from the Sierra Nevada Batholith, California[END_REF][START_REF] Teyssier | Fabric stability in oblique convergence and divergence[END_REF] with ~N-S horizontal shortening and subvertical principal stretch [START_REF] Cochelin | Vertical strain partitioning in hot Variscan crust: Syn-convergence escape of the Pyrenees in the Iberian-Armorican syntax[END_REF][START_REF] Cochelin | Structural inheritance in the Central Pyrenees: the Variscan to Alpine tectonometamorphic evolution of the Axial Zone[END_REF].

The partially molten mid-lower crust is observed in domes (Fig. 1a). The interpolated foliation trajectories, both at regional and local scales (see Figs. 7a & 13b), display circular to concentric trajectories within and around the domes. At regional scale, the ca. 15-20 km periodicity in the cleavage striking variations observed in the experimental variogram (Fig. 6a) shows that the domes display almost constant widths and could represent the periodical average map interdistance between them. We infer that this is related to a mechanical instability effect in the Variscan crust during original doming and exhumation of mid-lower crust. Smooth cleavage and lineation trajectories within gneiss domes (e.g. the Canigou massif) point to a rather homogeneous strainfield in the mid-lower crust with mainly longitudinal to slightly oblique stretching (Fig. 15) within domes' core.

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

Figure 15 here

In order to study the transitions from the upper crust to the mid-lower crust strains, one must consider the overall cleavage measurements dataset. At first glance, the nested structure of the regional-scale experimental variogram (Fig. 6a, grey dots/black curve) suggests that the Axial Zone includes distinct strainfields. Indeed, when considered separately, the cleavage measurements from the upper-and the mid-lower crust give two different variograms with distinct trends and parameters (Fig. 6a). Yet, if the two strain fields evidenced would display directions cross-cutting relationships (in particular within the contact/transition zones), the variogram computed for the overall data (in black, Fig. 6a) would display a significantly higher nugget value than the two others taken individually; which is not the case here. Besides, cleavage of the upper-crust has very steady orientations (i.e. both in strike and dip) and a dominant longitudinal collinear trend across the overall Axial Zone but some "disturbances" are linked to the gneiss domes. Thus, the finite strain orientation of upper-crust is dependent on the mid-lower crust fabric orientation when approaching domes' envelops, suggesting coeval deformation.

While shallow parts of the Variscan crust of the Pyrenees was affected by homogeneous thickening, the mid-lower crust was contemporarily affected by homogeneous longitudinal stretching (Figs 11 & 15) reflecting lateral horizontal flow, which is also compatible with NS horizontal shortening. This suggests vertical coupling during deformation [START_REF] Tikoff | Clutch tectonics and the partial attachment of lithospheric layers[END_REF].

As proposed by [START_REF] Cochelin | Vertical strain partitioning in hot Variscan crust: Syn-convergence escape of the Pyrenees in the Iberian-Armorican syntax[END_REF], such a vertical strain partitioning is typical of hot and ultra-hot orogens, where lithosphere has been significantly weakened (Cagnard et al., 2006a,b, Cruden et al., 2006[START_REF] Chardon | Flow of ultra-hot orogens: A view from the Precambrian, clues for the Phanerozoic[END_REF].

When considering the overall cleavage direction dataset -including upper-crust, plutons and gneiss domes -the interpolated trajectory map underline a quite heterogeneous pattern, with

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT
27 strain localization into regional scale, anastomosed longitudinal bands (Fig. 7a). These bands correspond to dextral-reverse shear zones observed in the field and mapped (e.g. [START_REF] A N U S C R I P T A C C E P T E D Accepted Manuscript 33 Carreras | Zooming on Northern Cap de Creus shear zones[END_REF][START_REF] Cochelin | Vertical strain partitioning in hot Variscan crust: Syn-convergence escape of the Pyrenees in the Iberian-Armorican syntax[END_REF]. The localization of these anastomosed shear zones seems to be intimately linked to the emplacement of gneiss domes and syn-kinematics plutons, as shear zones border these massifs (Fig. 7a). Domes and plutons must have acted as regional heterogeneities while the crust was affected by NS horizontal shortening. Therefore, we consider that these regional shear zones nucleated because of the emplacement of domes and plutons -themselves driven by mechanical instability during regional-scale shortening. Such an interpretation is further supported by geochronological and field structural studies arguing that steep shear zones were activated during and after the emplacement of gneiss domes and plutons [START_REF] Vilà | LP-HT metamorphism in a late orogenic transpressional setting, Albera Massif, NE Iberia: implications for the geodynamic evolution of the Variscan Pyrenees[END_REF][START_REF] Denèle | The Hospitalet gneiss dome (Pyrenees) revisited: lateral flow during Variscan transpression in the middle crust[END_REF][START_REF] Mezger | Transpressional tectonic setting during the main Variscan deformation: evidence from four structural levels in the Bossòst and Aston-Hospitalet mantled gneiss domes, central Axial Zone, M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT 37 Pyrenees[END_REF]Druguet et al., 2014;[START_REF] Cochelin | Vertical strain partitioning in hot Variscan crust: Syn-convergence escape of the Pyrenees in the Iberian-Armorican syntax[END_REF][START_REF] Van Lichtervelde | U-Pb geochronology on zircon and columbite-group minerals of the Cap de Creus pegmatites, NE Spain[END_REF].

At second order, directional-variograms computed for the upper-crust (Fig. 9a) document ~15-20km periodical undulations of the cleavage directions following a NE-SW spatial direction.

Likewise, the corresponding interpolated trajectories display asymmetrical strain gradients (sigmoidal undulations) compatible with a dextral component of simple shear along few N120-140°E trending bands (Fig. 7a), which are characterized by domains of more shallowly plunging lineations (Fig. 11-left). Nevertheless, these bands appear poorly significant at regional scale as the trajectories are mostly collinear and symmetrical. Thus, these undulations reflect horizontal strain partitioning in response to the emplacement of local heterogeneities in a more homogeneous general strainfield. Therefore, our results show that the dextral shearing component of deformation has been over-stressed in previous studies because most of them were focused on syn-kinematic plutons, and had had a minor impact on the strainfield.

The variogram analysis, a tool to investigate the rheological behavior of the crust

As a preliminary step, the variogram calculation and analysis is critical to detect and quantify significant variations (including the maximum correlation distance, variance and behavior at short distances) of the studied regionalized variable as well as potential anisotropy. We propose hereafter that, whatever the scale considered, the best fitted trend of variograms can be used to characterize the rheological behavior or the tectonic history of the geological object/domain investigated.

Amplitude of the variations in the measured values is evaluated by the sill (general variations; variance) and the nugget values in an experimental variogram. On the other hand, whatever those variations, the behavior to the origin of the variograms allows quantifying the way the variable varies with distance. In this study, gaussian or sinusoidal variograms have been computed for structures such as syn-kinematic plutons and gneiss domes. The very low original slope of these variogram functions points to limited and very progressive variations between structural orientation measurements, as exemplified by the Mt-Louis pluton and the Canigou gneiss dome (Figs. 12a & 14b). Indeed, the Mt Louis pluton, displays highly continuous and progressive strain variations and non-localized simple shear corridors. As shearing of the pluton happened during its cooling [START_REF] Bouchez | Two-stage deformation of the Mont-Louis-Andorra granite pluton (Variscan Pyrenees) inferred from magnetic susceptibility anisotropy[END_REF], we argue that the trend of the corresponding variogram can be diagnostic of such homogeneous/continuous deformation through time in case of weak mechanical behavior such as cooling magmas (Fig. 16). Similarly, in the Canigou gneiss dome, stretching lineations formed during pervasive metamorphism and crustal flow show continuous and smooth variations in trends, as illustrated by a sinusoidal variogram (Fig. 14b). To us, this phenomenon may also reflect a continuous ductile deformation at high temperature, which remains undisturbed during late Variscan folding and post-Paleozoic faulting (Fig. 16). 16. At regional scale, the cleavage orientations from the whole upper crust (metasediments) reflect interferences between a regional and consistent trend and the domes and plutons. For upper crust, the corresponding variogram clearly displays an exponential trend with a quite short range distance (Fig. 16). Domes and plutons can be regarded as geological discontinuities in the upper crust, explaining the exponential variogram trend for the upper crust dataset (Fig. 16) and the resulting interpolated trajectories showing local bending and reorientations around them (Fig. 10). Finally, the spherical trend of the omnidirectional-variogram that characterizes the whole Axial Zone reflects multi-scale interferences in a homogeneous strainfield recorded through time in the basement. This variogram corresponds to the imbrication, the convolution of all the variograms including local to regional variations induced by all the geological and structural objects/domains that compose the exposed crust (Fig. 16).

Concluding points

The present study shows that geostatistics is a powerful tool to model continuous strain trajectories from local ductile structure measurements. This inversion method can be considered -Separating (i) significant local to regional variations in strainfields from (ii) variations that cannot be considered and modeled from structural data. Variogram analysis allows determining the spatial scale (maximum distance) of influence of the tectonic and/or rheological controls on strain spatial variations. To the opposite, background noise or very local variations are integrated to the nugget effect and not modeled in computed trajectories for instance. At the scale of the Axial Zone of the Pyrenees, such nonsignificant variations make ~50% of the total variance of the dataset.

-Investigating strainfields at different scales, from a same structural dataset. By separating trends which are significant at a given scale from variations regarded as some noise or too "local", different sub-scales can thus be explored within a study area. When zooming on very local sub-domains, variations initially considered as "noise" in the regional-scale nugget effect can become significant and be interpreted. In the case of the Pyrenean Axial Zone, local interpolated map of the Canigou gneiss dome or Mt-Louis pluton display more detailed trajectories than at regional-scale, where internal variations can be regarded as irrelevant. b) Zoom on the regional trajectories obtained around the Mt Louis pluton (Fig. 7a); c) Foliation trajectories obtained using parameters deduced from the omni-directional variogram of the pluton; d) Foliation trajectories obtained manually by [START_REF] Bouchez | Two-stage deformation of the Mont-Louis-Andorra granite pluton (Variscan Pyrenees) inferred from magnetic susceptibility anisotropy[END_REF] on the basis of ASM measurement and field observations. dome (in red); Trajectories around the gneiss dome (in grey) are from the regional interpolation (Fig. 7a). c) Location of high dispersion of directional data in the Canigou gneiss dome combining of zoom of the R map (Fig. 7b), local trajectories (Fig. 13b) and location of the main fault zones (Fig. 2).

  . The Canigou and Mt-Louis massifs can be regarded as demonstrative examples of gneiss dome and pluton of the Axial Zone, analyses of local scale strainfields are made on these massifs in the present study.

  3 & 4), details of the different steps of our approach developed are given below(3.1)

  for circular data. Kriging weights are computed based on a given theoretical variogram has been adjusted -for each data set -to the parameters deduced from the detailed variogram analysis (parameters include the function type, the range value, the

  may be oversighted compared to those from steep cleavage domains as this analysis of cleavage direction was performed in 2D (i.e. map view). Indeed, strike value of subhorizontal fabrics (i.e. < 20°) may significantly vary in response to local orientation perturbations.

  Figure 3 here

Figure 5

 5 Figure 5 here

  Figure 7 here

  Figure 8 here

  to the fabric. Following the exponential best fit function on the 40 first kilometers, the range of the N-S directed variogram and the characteristic distance of influence of the corresponding phenomenon could be reached around 80-100km. This variation can be related to the regional N-S fan of schistosity (section 4.1).

  Figure10 here

  Figure16 here

  complementary to forward analog or thermo-mechanical modelling to investigate strainfields and the interplay of different mechanical controls leading to the development of finite strain in the crust. In particular, this study shows that the approach developed allows:

-Figure 1 :

 1 Figure captions

Figure 3 :

 3 Figure 3: Summary of the following steps to compute foliation trajectories. See in the text for

Figure 4 :

 4 Figure 4: Summary of the following steps to compute lineation trajectories, using lineation data

Figure 5 :

 5 Figure 5: a) Histogram of cleavage directional data, b) histogram of cleavage directional data in

Figure 6

 6 Figure 6: a) Omni-directional variograms computed for schistosity directional data, with the

Figure 7 :

 7 Figure 7: a) Cleavage trend trajectories computed from the interpolated cleavage direction in the

Figure 8 :

 8 Figure 8: a) Interpolated dip of cleavage in the Axial Zone. Gneiss domes are located in red and

Figure 9 :

 9 Figure 9: Directional variograms computed in the upper crust: a) for azimuth of cleavage along

Figure 10

 10 Figure 10: a) Cleavage trend trajectories computed from the interpolated cleavage direction in the

Figure 11 :

 11 Figure 11: Ternary diagrams (Pitch-Plunge-Dip) illustrating the crustal flow mode in orogens,

Figure 12

 12 Figure 12: a) Omni-directional variogram computed for directional data in the Mt Louis pluton;

Figure 13 :

 13 Figure 13: Zoom on the Canigou gneiss dome. a) Omni-directional variogram calculated for

Figure 14

 14 Figure 14: a) Stereogram of mineral-stretching lineations in the Canigou gneiss dome; b) Omni-

Figure 15

 15 Figure 15: a) Synthesis of direction of stretching within gneiss domes in the Axial Zone,

  Fig. 1b

  Fig. 2

Acknowledgments

This work was supported by the BRGM through the Référentiel Géologique de la France program (RGF). We warmly thank B. Bourgine (BRGM -France) for fruitful discussion about variogram analysis. We also thank two anonymous reviewers for helping to improve the manuscript and William Dunne for the editorial handling.