Dániel Kiss
email: daniel.kiss@unil.ch

Yuri Podladchikov

Thibault Duretz

Stefan M Schmalholz

Spontaneous generation of ductile shear zones by thermal softening: localization criterion, 1D to 3D modelling and application to the lithosphere

Keywords: Ductile shear zone, thermal softening, strain localization, shear heating, localization criterion, subduction initiation

The generation of ductile shear zones is essential for the formation of tectonic plate boundaries, such as subduction or strike-slip zones. However, the primary mechanism of ductile strain localization is still contentious. We study here the spontaneous generation of ductile shear zones by thermal softening using thermo-mechanical numerical simulations for linear and power-law viscous flow in one-dimension (1D), 2D and 3D. All models are velocity-driven. The 1D model exhibits bulk simple shear whereas the 2D and 3D models exhibit bulk pure shear.

The initial conditions include a small temperature perturbation in otherwise homogeneous material. We use a series of 1D simulations to determine a new analytical formula which predicts the temperature evolution inside the shear zone. This temperature prediction requires knowledge of only the boundary velocity, flow law and thermal parameters, but no a priori information about the shear zone itself, such as thickness, stress and strain rate. The prediction is valid for 1D, 2D and 3D shear zones in bulk pure and simple shear. The results show that shear heating dominates over conductive cooling if the relative temperature increase is > 50 • C.

Introduction

The spontaneous generation of shear zones in ductile rocks is fundamental for the formation of tectonic plate boundaries, such as subduction and strike slip zones, or the generation of tectonic nappes during orogenic wedge formation.

We refer here to spontaneous generation of a shear zone when the fundamental shear zone parameters, such as thickness, shear stress and strain rate, are not a priori prescribed by the natural or model configuration. We refer to ductile deformation when the deformation behaviour is described mathematically by a relation between stress and strain rate, such as by flow laws for diffusion, dislocation or Peierls creep (i.e. low temperature plasticity). The conversion of dissipative work into heat, the related local temperature increase and the associated decrease of temperature dependent rock viscosities has frequently been suggested as a cause of spontaneous strain localization and shear zone formation in the litho-sphere [START_REF] Yuen | Shear deformation zones along major transform faults and subducting slabs[END_REF][START_REF] Regenauer-Lieb | Rapid conversion of elastic energy into plastic shear heating during incipient necking of the lithosphere[END_REF][START_REF] Leloup | Shear heating in continental strike-slip shear zones: model and field examples[END_REF][START_REF] Kaus | Initiation of localized shear zones in viscoelastoplastic rocks[END_REF][START_REF] Takeuchi | Dynamic models of interseismic deformation and stress transfer from plate motion to continental transform faults[END_REF][START_REF] Thielmann | Shear heating induced lithospheric-scale localization: Does it result in subduction?[END_REF][START_REF] Duretz | Shear heating-induced strain localization across the scales[END_REF][START_REF] Jaquet | Dramatic effect of elasticity on thermal softening and strain localization during lithospheric shortening[END_REF][START_REF] Moore | Scaling of viscous shear zones with depthdependent viscosity and power-law stress-strain-rate dependence[END_REF]. We refer here to this thermally controlled strain localization mechanism as thermal softening. Despite its fundamental thermo-mechanical feasibility [START_REF] Hersey | Note on heat effects in capillary flow[END_REF][START_REF] Brinkman | Heat effects in capillary flow i[END_REF][START_REF] Gruntfest | Thermal feedback in liquid flow; plane shear at constant stress[END_REF], shear heating and thermal softening is still contentious as important softening mechanism causing strain localization in ductile rock [START_REF] Regenauer-Lieb | The initiation of subduction: criticality by addition of water?[END_REF][START_REF] Platt | Grainsize evolution in ductile shear zones: Implications for strain localization and the strength of the lithosphere[END_REF][START_REF] Bercovici | Mechanisms for the generation of plate tectonics by two-phase grain-damage and pinning[END_REF][START_REF] Ghazian | A numerical investigation of continental collision styles[END_REF][START_REF] Gueydan | Strain weakening enables continental plate tectonics[END_REF][START_REF] Platt | Influence of shear heating on microstructurally defined plate boundary shear zones[END_REF]. This is, for example, different from physics-based models of friction in rock where essentially all potential processes causing significant friction weakening are considered to be related to shear heating, such as "flash heating", thermal pressurization or temperature controlled chemical/phase changes, including melting and formation of pseudotachylites [START_REF] Sibson | Generation of pseudotachylyte by ancient seismic faulting[END_REF][START_REF] Fialko | Fusion by earthquake fault friction: Stick or slip?[END_REF][START_REF] Brown | melt welt'mechanism of extreme weakening of gabbro at seismic slip rates[END_REF][START_REF] Aharonov | A physics-based rock friction constitutive law: Steady state friction[END_REF]. For ductile strain localization, proposed alternative mechanisms not related to shear heating are, for example, grain size reduction [START_REF] Bercovici | Mechanisms for the generation of plate tectonics by two-phase grain-damage and pinning[END_REF][START_REF] Platt | Influence of shear heating on microstructurally defined plate boundary shear zones[END_REF], reaction-weakening caused by infiltration of fluids along precursor brittle faults [START_REF] White | Transformation-and reaction-enhanced ductility in rocks[END_REF][START_REF] Mancktelow | The control of precursor brittle fracture and fluid-rock interaction on the development of single and paired ductile shear zones[END_REF] or fabric development in rock with significant mechanical heterogeneities [START_REF] Montési | Fabric development as the key for forming ductile shear zones and enabling plate tectonics[END_REF]. Out of the different mechanisms proposed for ductile strain localization, shear heating and thermal softening (1) must occur in nature since dissipative deformation generates heat and rock viscosity is temperature dependent and (2) requires the least assumptions since no knowledge concerning grain size reduction and growth, fluid flow, reaction kinetics or mechanical heterogeneities required for fabric evolution is needed. Hence, ther-mal softening as mechanism itself is actually not contentious, but whether thermal softening alone can be significant enough to generate shear zones in ductile rock under natural conditions is debated, as well as its relative importance compared to other localization mechanisms.

A long-lived argument against the significance of thermal softening during ductile deformation in the lithosphere is that many natural shear zones with thickness ranging from hundreds of meters to several kilometers do not indicate a sharp change in temperature between the little-deformed wall rock and the highlydeformed shear zone. This argument persists, despite the fact that several thermomechanical studies have shown that even if a shear zone is caused by thermal softening, there are only small temperature gradients between the shear zone and the wall rock (e.g., [START_REF] Yuen | Shear deformation zones along major transform faults and subducting slabs[END_REF][START_REF] Takeuchi | Dynamic models of interseismic deformation and stress transfer from plate motion to continental transform faults[END_REF][START_REF] Schmalholz | Shear zone and nappe formation by thermal softening, related stress and temperature evolution, and application to the alps[END_REF][START_REF] Mako | Quantifying magnitudes of shear heating in metamorphic systems[END_REF]. Another argument against the importance of thermal softening is that the required shear stresses or the required strain rates are too large for typical lithospheric deformation conditions (e.g. [START_REF] Platt | Influence of shear heating on microstructurally defined plate boundary shear zones[END_REF].

To test the validity of the above arguments against thermal softening and to quantify thermal softening, we use a thermo-mechanical numerical model of ductile rock deformation based on the conservation equations of continuum mechanics and apply constitutive equations for ductile creep. We perform a scaling analysis with results of a one-dimensional (1D) model for which simple shearing is controlled by a boundary velocity and strain localization can be triggered by a temperature, hence viscosity, perturbation in the model center. We apply a temperature perturbation because such perturbation diffuses away if shear heating is not efficient. In contrast, a viscosity perturbation (i.e. perturbation of material properties) would remain even if shear heating is insufficient and would, hence, always generate a shear zone in the 1D model, with a thickness of the initial viscosity perturbation. The model configuration is based on the model of [START_REF] Yuen | Shear deformation zones along major transform faults and subducting slabs[END_REF] because for this configuration "No a priori assumption about slip-zone width or shear-stress magnitude is necessary; the thermal-mechanical structure of the slip zone evolves in time and all its characteristics are self consistently determined" [START_REF] Yuen | Shear deformation zones along major transform faults and subducting slabs[END_REF]. These model features are essential to study spontaneous generation of shear zones in homogeneous material. [START_REF] Yuen | Shear deformation zones along major transform faults and subducting slabs[END_REF] considered linear viscous flow laws only whereas we also consider power-law viscous flow laws to apply our results to lithospheric dislocation creep flow laws which exhibit power-law stress exponents typically between 2 and 4 (Table 1).

The aims of our study are to (1) quantify the temperature increase required for spontaneous shear zone generation, (2) quantify the relation between the width of the temperature variation across the shear zone and the width of the corresponding finite strain variation, (3) quantify stresses, velocities and strain rates required for shear zone generation, (4) derive an analytical formula which predicts the temperature inside a shear zone without a priori knowledge of the thickness, stress and strain rate of the shear zone, (5) compare the 1D model for bulk simple shear with 2D and 3D models for bulk pure shear and (6) evaluate the importance of thermal softening for ductile strain localization in the lithosphere.

Mathematical and numerical model

Governing system of equations

We assume incompressible viscous deformation in the absence of gravity and inertial forces. The governing system of equations is

- ∂v i ∂x i = 0 (1) -δ ij ∂P ∂x i + ∂τ ij ∂x j = 0 (2) ρc p ∂T ∂t - ∂ ∂x i λ ∂T ∂x i -τ ij ˙ ij = 0 (3) τ ij = 2µ eff ˙ ij (4) µ eff ( ˙ II , T ) = A -1 n ˙ 1 n -1 II exp Q nRT (5)
where equation ( 1), ( 2) and ( 3) are the equations for conservation of mass, linear momentum and energy, respectively, equation ( 4) is the creep flow law (constitutive equation) and equation ( 5) states the effective viscosity. The indices i and j correspond to coordinate axes 1, 2 and 3 and repeated indices imply summation.

In equation (3) we assume that all dissipative work is converted to heat (so-called

Taylor-Quinney coefficient is 1.0) since we do not consider grain size reduction which consumes typically only a minor fraction of the dissipative work [START_REF] Herwegh | From transient to steady state deformation and grain size: A thermodynamic approach using elasto-visco-plastic numerical modeling[END_REF][START_REF] Thielmann | Intermediate-depth earthquake generation and shear zone formation caused by grain size reduction and shear heating[END_REF]. x i are the components of the spatial coordi-

nates [m], t is the time [s], v i are components of the velocity vector [m.s -1 ], δ ij is the Kronecker delta, τ ij are components of the deviatoric stress tensor [Pa], ρ is density [kg.m -3 ], c p is heat capacity at constant pressure [J.K -1 ], T is tempera- ture [K], λ is thermal conductivity [W.K -1 .m -1 ], µ eff is effective viscosity [Pa.s],
˙ ij are components of the deviatoric strain rate tensor [s -1 ], ˙ II is the square root of the second invariant of the strain rate tensor [s

-1 ], n is the power law exponent [ ], A is the pre-exponential factor [Pa n .s -1 ], Q is activation energy [J.mol -1 ] and R is the universal gas constant [J.mol -1 .K -1 ].
Initial and boundary conditions are

T ( √ x i x i > r, t = 0) = T 0 and T ( √ x i x i ≤ r, t = 0) = T 0 + ∆T 0 (6) q i (x i = [0 or L i ], t) = 0 (7) v i (x j = 0, t) = 0 and v i (x j = L j , t) = ∆v i (8)
where T 0 is initial temperature [K], ∆T 0 is the value of the intial temperature perturbation [K], q i are components of the heat flux vector [W.m -2 ], L i is the total size of the model domain [m] in the different spatial directions and ∆v i is the farfield velocity difference [m.s -1 ] in the different spatial directions. The material parameters are homogenous and the initial temperature is constant except a small temperature perturbation, ∆T 0 , in a region around the model centre, x i = 0, whose size is specified with radius r (eq. 6, Fig. 1). This thermal perturbation mimics any kind of small variation of strength or thermal properties which are always present in natural rocks. The model is thermally insulated (eq. 7). For simple shear type deformation the model is kinematically driven by constant far-field boundary velocities (eq. 8). For pure shear type deformation in 2D and 3D only velocities normal to the boundaries are defined (i.e. i = j for eq. 8), otherwise free slip boundary conditions are used, so shear stresses are zero at the boundaries.

Numerical method

The system of non-linear equations (Eq. 1-5) is discretized on a regular Cartesian staggered grid. The problem is solved by a pseudo-transient iteration or relaxation scheme [START_REF] Versteeg | An Introduction to Computational Fluid Dynamics: The Finite Volume Method[END_REF][START_REF] Duretz | Resolving thermomechanical coupling in two and three dimensions: spontaneous strain localization owing to shear heating[END_REF]. The thermo-mechanical equations are recasted in the following form:

dP dω = ∂v i ∂x i , dv i dω = ∂τ ij ∂x j - ∂P ∂x i , dT dω = ρc p ∂T ∂t - ∂ ∂x i λ ∂T ∂x i -τ ij ˙ ij . (9

1D, 2D and 3D model configurations

The 1D model domain extends orthogonally across the shear zone and velocities are orthogonal to the model domain. The 1D model is driven by a velocity difference at the two model boundaries which imposes a bulk simple shear deformation (Fig. 1a). To test whether results of the 1D model are applicable to 2D and 3D shear zones, we perform also 2D and 3D numerical simulations for bulk pure shear. For the 3D model shortening occurs in one horizontal direction and extension in the vertical direction while the bulk extension in the second horizontal dimension is zero (Fig. 1c). The initial temperature perturbation has the shape of a quarter circle in the 2D model (Fig. 1b) and one eight of a sphere in the 3D model (Fig. 1c). To compare the 2D and 3D results with the 1D results we record the temperature and shear velocities along a profile line, with coordinates

x , which is orthogonal to the 2D and 3D shear zones. These results are directly comparable with the results of the corresponding 1D model (Fig. 1).

Fundamental features of 1D shear zone evolution

There are two end-member solutions of the numerical model: (1) The velocity field converges to homogeneous simple shear in the entire model domain, there is no strain localization and the temperature increases homogeneously in the model domain due to bulk shear heating.

(2) The temperature increases locally in the model center, which causes strain localization and the generation of a shear zone that is much thinner than the model domain. We show in the following fundamental features of solution (2) for representative simulations.

For simulations with linear viscosity the temperature increases in the shear zone during a transient stage and then reaches a constant temperature (Fig. 2a), in agreement with [START_REF] Fleitout | Thermal and mechanical evolution of shear zones[END_REF]. This temperature is independent on the initial temperature perturbation and model width (Fig. 2a). For simulations with power-law viscous flow laws the temperature also increases in the shear zone during a transient stage and then reaches a quasi-constant temperature (Fig. 2b). In contrast to the linear viscous model the temperature in the shear zone does not reach a strictly-constant value, but the temperature is slightly increasing with ongoing deformation, referred to here as quasi-constant (Fig. 2b). This quasi-constant temperature is also independent on the initial temperature perturbation and the model width (Fig. 2b). For such quasi-constant temperature in the shear zone center, shear heating must be locally balanced by thermal conduction (Eq. 3). [START_REF] Fleitout | Thermal and mechanical evolution of shear zones[END_REF] showed that constant boundary velocities guarantee that this balance is always reached. Hence, velocity-driven shearing of a dominantly viscous medium does not lead to a thermal runaway for which temperatures would increase exponentially, in an unbounded way.

The presented model results correspond to simulations with significant shear localization in the model center and show that a modest temperature rise of 100

o C can result in shear localization due to thermal softening (Fig. 2a andb).

The shear stress is spatially constant in the 1D model at each instant of time,

it is largest at the onset of deformation and decreases with time due to progressive temperature increase due to shear heating (Fig. 2c). When the maximal temperature reaches a quasi-constant value then also the stress reaches a quasi-constant value.

The characteristic width of the temperature variation across the shear zone, referred to here as thermal thickness, is not prescribed a priori but controlled by the thermo-mechanical process [START_REF] Duretz | Physics-controlled thickness of shear zones caused by viscous heating: Implications for crustal shear localization[END_REF][START_REF] Duretz | Shear heating-induced strain localization across the scales[END_REF]. After the maximal temperature has reached its (quasi-)constant value, the thermal thickness is increasing proportional with the square root of time due to thermal conduction (Fig. 3a). In the following, we distinguish between the thermal thickness and the finite strain thickness that is determined by the width of the finite strain profile (Fig. 4a). The thermal thickness is defined as the width of the temperature profile at this temperature, which is half between the maximal temperature in the shear zone center and the minimal ambient temperature far away from the shear zone (Fig. 3a). During the transient stage of temperature increase the thermal thickness is typically decreasing (Fig. 3b). The thermal thickness evolution is essentially unaffected by the initial temperature perturbation and model size, similar to the evolution of the maximum temperature (Fig. 3b).

The temperature profile is significantly wider than the corresponding finite strain (γ) profile across the shear zone (Fig. 4a), which agrees with results of [START_REF] Takeuchi | Dynamic models of interseismic deformation and stress transfer from plate motion to continental transform faults[END_REF] for strike slip zones and of [START_REF] Schmalholz | Shear zone and nappe formation by thermal softening, related stress and temperature evolution, and application to the alps[END_REF] for thrust-type shear zones. γ is calculated by time integration of the shear strain rates. The finite strain thickness is measured in the same way as the thermal thickness, that is, the width of the γ-profile at the value of γ half between the maximum value and the far-field value at the model boundary (Fig. 4a). The ratio of thermal to finite strain thickness increases during the transient phase of temperature increase. Once the temperature has reached its (quasi-)constant value this ratio converges towards a constant value. This shows that the finite strain and thermal thickness are linked, both are controlled by thermal conduction. After the transient phase, the thermal thickness is nearly one order of magnitude (factor 6 to 8) larger than the finite strain thickness (Fig. 4b). The example presented in figure 4a shows that a significant decrease of γ from ca. 17, in the shear zone center, to 1 is associated with only a minor temperature decrease from ca. 800 o C to 760 o C (i.e. only a 20% decrease).

Predictive scaling relationships and localization criterion

We performed a series of 1D numerical simulations with a geologically applicable range for all independent model parameters: ∆v, A, ρ, c p , λ, T 0 , n and Q.

For each parameter we used several representative values which were evenly dis-tributed within the chosen range (e.g.: n = {1, 2, 3, 5, 6}). To test the usefulness of several different sets of independent scales, we performed these simulations using the dimensional form of equations ( 4), ( 5) and ( 9). We run 1D simulations with all parameter combinations and recorded characteristic parameter values (e.g.: T max , µ min ) in regular intervals during shear zone evolution. We recorded data from more than 45'000 simulation stages (i.e. at specific times) from ca. 2'000 simulations.

Thermal thickness of shear zones

All simulations show that shear zones are widening proportional to the square root of time (Fig. 3) and that widening is controlled by heat conduction. Two fundamental types of conductive heat transfer between shear zone and surrounding region can be distinguished: (1) If there is no significant shear heating, then the initially higher temperature in the model center is decreasing with respect to the far-field temperature during shearing and the temperature evolution in the model can be approximated with an analytical solution for Gaussian cooling of an initial Dirac delta temperature profile (Fig. 5a). The spatial and temporal evolution of temperature can then be described by the equation:

∆T (x, t) = 1 √ 4πκt exp - x 2 4κt (10)
If the maximum temperature is in the model center, at x = 0, then the half width of the temperature profile is given by the value of x for which ∆T (x, t) = 0.5 ∆T max (t):

∆T max (t) 2 = 1 2 √ 4πκt = 1 √ 4πκt exp - x 2 4κt → (11) 0.5 = exp - x 2 4κt → x = -4 ln(0.5)κt ≈ 1.67 √ κt (12)
The corresponding full width of the Gaussian temperature profile, W G is then:

W G ≈ 3.34 √ κt ⇐⇒ W G / √ κt ≈ 3.34 (13)
(2) If there is significant shear heating, then the temperature in the model center reaches a (quasi-)constant value after a transient period (Fig. 2a andb). The temperature evolution in model can then be described with a half space heating model in which the temperature is kept constant at one side, representing the shear zone center, and the far-field temperature is the initial temperature at the model boundary. The analytical solution for the temperature evolution for such scenario quantifies the heating of a half-space, representing the region adjacent to the shear zone and is given by an error function solution (Fig. 5a):

∆T (x, t) = ∆T max erfc x 2 √ κt (14)
If the maximum temperature is in the model center, at x = 0, then the half width of the temperature profile is given by the value of x for which ∆T (x, t) = 0.5 ∆T max (t):

∆T max 2 = ∆T max erfc x 2 √ κt → 0.5 = erfc x 2 √ κt (15)
Using the approximation erfc(0.48) ≈ 0.5 yields:

0.48 ≈ x 2 √ κt → x ≈ 0.96 √ κt (16)
The corresponding full width of such temperature profile, W E , is then:

W E ≈ 1.92 √ κt ⇐⇒ W E / √ κt ≈ 1.92 (17)
If thermal evolution during shearing is dominated by conductive cooling, then the width of the temperature across the shear zone will grow according to W G , and according to W E if thermal evolution is dominated by shear heating. We plotted the dimensionless widths, scaled by √ κt, of the numerically calculated temperature profiles versus the maximal temperature difference (i.e. maximum temperature in the shear zone center minus initial temperature, ∆T) recorded in the numerical simulations (Fig. 5b). For insignificant shear heating, ∆T < ca. 20

o C, the temperature profile is widening according to W G (Fig. 5b). For significant shear heating, ∆T > ca. 100 o C, the temperature profile is widening according to W E (Fig. 5b). Between 20 and 100 o C for ∆T there is a transition zone where the thicknesses are in between W E and W G . The boundary between the two heat transfer domains occurs at ∆T ≈ 40 o C. The results, hence, indicate that a temperature increase of at least 40 o C in the shear zone is required so that shear heating dominates the heat transfer across the shear zone.

Maximum temperature of shear zones

After a transient phase the temperature in the shear zone is (quasi-)constant.

For such quasi steady state, heat production and conduction are essentially balanced in the shear zone, that is:

0 ≈ ρc p ∂T ∂t x=0 = λ ∂ 2 T ∂x 2 x=0 + µ ∂v ∂x 2 x=0 . ( 18
)
For this quasi steady state, we want to determine a scaling relationship between the term representing diffusion (with λ in Eqn. 18) and the term representing shear heating (with µ). Such relationship can be of the form:

λT c ≈ a µ c v 2 c ( 19
)
where the subscripts c indicate a characteristic value of the corresponding parameter, and a is a proportionality constant. The characteristic length scale has been dropped because it has the same power in both terms for diffusion and heat pro-

T max ≈ ∆v e µ min Q λR ( 20
)
where e ≈ 2.72 is the Euler number, and e -1 is the proportionality constant a (Fig. 6a). All parameters in this formula correspond to a specific time during shear zone evolution. Equation 20 is useful in applications where a shear zone viscosity (µ min as a function of T max and ˙ IImax or τ II ) can be constrained, for example, for rock deformation experiments. For most natural shear zones viscosities cannot be easily constrained. This is because, for example, for power-law viscous flow knowledge of the strain rate is required to determine the effective viscosity. We approximate the strain rate by the ratio of ∆v/ √ κt assuming that √ κt provides a representative value for the shear zone thickness. Using then the dataset from all the 1D simulations we determine a formula to fit the shear zone viscosity:

µ min ≈ 1.28 e 2 λQ ∆v 2 n 2 R ln ∆v 2 nR λQ A -1 n ∆v √ κt 1 n -1 + 1.1 -2 . ( 21
)
For n = 1 the term with the approximate strain rate disappears. If we substitute the approximation of the shear zone viscosity in equation ( 20) we get:

T max ≈ -1.13 Q nR ln ∆v 2 nR λQ A -1 n ∆v √ κt 1 n -1 + 1.1 -1 (22)
Crosschecking with the numerical results provides the minus sign for taking the square root of µ min in Eqn. 20. Equation 22predicts the maximum temperature in all numerically simulated shear zones with a maximal error of < 50 o C and with a root mean square error of only 20 o C (Fig. 6b). The prediction of T max using equation ( 22) does not require any a priori knowledge of the shear zone thickness, the stress, the strain rate and the effective viscosity in the shear zone.

The great advantage of equation ( 22) is, hence, that T max inside a shear zone can be estimated exclusively with flow law parameters (n, A and Q), thermal parameters (λ and c p ), the density (ρ), the applied boundary velocity difference (∆v) and the duration (t) of shearing.

Localization criteria

A possible criterion for shear localization is that shear heating must dominate the heat transfer between shear zone and the surroundings. Based on the results discussed in the previous section we suggest T max -T 0 > 50 o C as localization

criterion. An alternative criterion can be derived by separating the variables and the constant e in equation ( 20) and squaring both sides, which yields:

∆v 2 µ min Q λRT 2 max ≈ e 2 (23)
If the maximum temperature is replaced by the smaller initial temperature, T 0 , and the minimum viscosity by the larger initial viscosity, µ 0 , then a modified criterion for shear localization is:

Br 1 = ∆v 2 µ 0 Q λRT 2 0 > e 2 (24)
The dimensionless number on the left hand side is a particular version of the Brinkman number (Br 1 ). Several authors suggested different versions of the Brinkman number as criterion of shear localization, based on scaling analyses (e.g. [START_REF] Brinkman | Heat effects in capillary flow i[END_REF][START_REF] Gruntfest | Thermal feedback in liquid flow; plane shear at constant stress[END_REF][START_REF] Yuen | Shear deformation zones along major transform faults and subducting slabs[END_REF][START_REF] Brun | Strain heating and thermal softening in continental shear zones: a review[END_REF].

Another typical version of the Brinkman number (Br 2 ) is Br 1 divided by the Arrhenius exponent (Q/RT 0 ). The corresponding localization criterion is then:

Br 2 = ∆v 2 µ 0 λT 0 > 1 (25)
We plotted the values of both Br 1 and Br 2 versus the viscosity decrease in the shear zone center (µ min /µ 0 ) for all simulations. Both numbers are proportional to the viscosity decrease (Fig. 6c), hence they are useful criteria for strain localization. We prefer using the criterion T max -T 0 > 50 o C due to its simplicity, or the criterion based on Br 1 because it has been directly derived from the analytical formula (eq. 20).

Comparison of 1D, 2D and 3D shear zones

Equations 20 and 22 for predicting the temperature in the shear zone are based on a 1D model, which is driven by far-field simple shear. We apply this prediction to shear zones that develop in 2D and 3D models, which are driven by far-field pure shear, in order to test the general applicability of the temperature prediction (Fig. 7). The rate of temperature increase in the shear zone is the largest for the shows that the thickness of the shear zones are essentially identical. Therefore, 1D, 2D and 3D shear zones caused by thermal softening under both far-field pure and simple shear exhibit the same thermo-mechanical characteristics.

Application to dislocation creep flow laws

We apply equation 22 to predict the maximum temperature in shear zones using typical flow laws for rock-forming minerals relevant to the lithosphere. Equation ( 22) depends on the duration of deformation. The typical, observed time scale of deformation varies as a function of shear velocity. To make the results for different velocities comparable we assume a characteristic shear strain of 20, which is the ratio of displacement and shear zone thickness (W γ ). As a first order estimate we use one tenth (Fig 4b) of the characteristic thermal width W T = 2 √ κt (Fig. 5b) as a shear zone width. With these relationships we can determine a representative characteristic time of the deformation to reach a shear strain of 20:

20 = γ c = t c ∆v 0.2 √ κt c → t c = 0.2 γ c √ κ ∆v 2 (26)
which yields t c ≈ 5 Ma for a velocity of ≈ 1 cm.yr -1 .

For typical plate tectonic velocities in the order of few centimeters per year, the flow laws for wet and dry olivine, and dry plagioclase are associated with maximum temperatures between 500 and 700 o C (Fig. 8). Such temperatures correlate well with typical temperatures of highly sheared basement nappes outcropping in orogens such as the Alps. Typical metamorphic peak temperature ranges between 500 and 650 o C [START_REF] Keller | Structural and metamorphic evolution of the camughera-moncucco, antrona and monte rosa units southwest of the simplon line, western alps[END_REF][START_REF] Manzotti | Exhumation rates in the gran paradiso massif (western alps) constrained by in situ u-th-pb dating of accessory phases (monazite, allanite and xenotime)[END_REF] and reach up to 800 o C in the Lepontine dome [START_REF] Nagel | Tertiary subduction, collision and exhumation recorded in the adula nappe, central alps[END_REF]. [START_REF] Takeuchi | Dynamic models of interseismic deformation and stress transfer from plate motion to continental transform faults[END_REF] provided a thorough study of the temperature anomalies around the San Andreas strike-slip fault. They used heat flow and surface deformation measurements to constrain their models. They conclude that a temperature increase of 160 to 375 o C, dependent on rheology, is expected at 20 km depth for a 4 cm.yr -1 long term average velocity difference. These values are in a good agreement with our prediction for dry anorthite and olivine and wet olivine flow laws. Flow laws for wet quartzite and Westerly granite provide maximal temperatures < 300 o C in the same velocity range (Fig. 8).

A recent study of [START_REF] Chu | Ultrafast eclogite formation via melting-induced overpressure[END_REF] provides well constrained information about the duration of deformation and maximal temperature of eclogite shear zones of the Taconic orogenic belt (New England). The eclogite bodies are hosted in feldspar rich felsic paragneiss. The authors conclude that the P -T history of the shear zones can be best explained by shear heating. Using the known displacement, related to the known deformation time, the shear velocities can be constrained to be between 25 -70 cm.yr -1 . For the inferred velocity range the applicable flow law for dry anorthite yields a good fit (Fig. 8).

Rocks that are commonly considered to result from significant shear heating are pseudotachylites. They are often associated to earthquakes, having typically slip velocities on the order of a m.s -1 [START_REF] Bizzarri | Analytical representation of the fault slip velocity from spontaneous dynamic earthquake models[END_REF]. Such fast deformation processes are commonly considered to be dominated by frictional deformation and, hence, frictional heating. However, recent progress in understanding of the physics of friction suggests viscous creep on grain contacts and asperities as the mechanism for velocity weakening of the friction coefficient, reported at high shear-velocity (≈ 1 m.s -1 ) rock deformation experiments for various rock types [START_REF] Aharonov | A physics-based rock friction constitutive law: Steady state friction[END_REF]. Moreover, a recent experimental study of high shear-velocity (≈ 1 m.s -1 ) deformation of calcite reports that such fast shear deformation is characterized by an initial frictional deformation followed by (quasi) steady-state viscous creep [START_REF] Pozzi | A new interpretation for the nature and significance of mirror-like surfaces in experimental carbonate-hosted seismic faults[END_REF]. In these experiments most of the strain is generated by (quasi) steady-state viscous creep. Since the (quasi) steady-state temperature developing in our models is path independent (Fig. 2), we can apply our viscous model result to estimate the temperatures in such small-scale and high-velocity shear zones regardless of which deformation mechanism dominates initially. To test the model-based temperature estimations we consider natural pseudotachylites from Corsica [START_REF] Andersen | Stress release in exhumed intermediate and deep earthquakes determined from ultramafic pseudotachylyte[END_REF]. The reported peak metamorphic temperature is at least 1750 o C. There are no reported constraints on shear velocities and, therefore, we apply peak slip velocities of 4 m.s -1 (that is typical for seismic events with a displacement around 1 m), as an upper limit [START_REF] Bizzarri | Analytical representation of the fault slip velocity from spontaneous dynamic earthquake models[END_REF]. To estimate a lower bound, we use the width of the main pseudotachylite vein (W = 1.23 cm), the displacement (d = 1 m, yielding γ c = d/W ≈ 80) along it and the scaling relationship between characteristic width of shear zones and the duration of deformation (eq. 16). Reordering equation (16) yields t = W 2 /(1.92 2 κ) ≈ 41 s (assuming κ ≈ 10 -6 m 2 .s -1 ). Using this estimate of duration we can estimate the slip velocity v = d/t ≈ 0.02 [m.s -1 ].

As such estimates have typically an order of magnitude uncertainty we take a five times lower value as a lower velocity bound (i.e. 4 mm.s -1 ). Within this wide velocity range, the predicted temperatures for most flow laws agree with the reported peak temperature (Fig. 8). Clearly, there are many uncertainties and simplifications related to this temperature estimate, but this estimate nevertheless indicates that peak temperatures reported for the considered pseudotachylites potentially could have been generated in viscous shear zones for typical slip velocities, in the order of 1 m.s -1 .

Whether shear heating causes shear localization depends on the initial temperature of the rock because localization will not occur if the ambient rock temperature at the onset of shearing is larger than the predicted maximum temperature.

We consider typical lithospheric geotherms and temperatures for the upper crust between 200 and 400 o C, for the lower crust between 400 and 600 o C and for the mantle lithosphere > 600 o C. For typical plate tectonic velocities, we calculate the predicted maximum temperature for different initial temperatures representing the ambient temperature at the onset of deformation (Fig. 9a). The temperature difference, ∆T , between predicted maximal temperature and initial, ambient, temperature indicates the intensity of shear heating and, hence, shear localization by thermal softening. Shear heating is significant for ∆T > 50 o C since for such values of ∆T the heat transfer between shear zone and wall rock is dominated by shear heating (Fig. 5b). For plate tectonic velocities of a few cm.yr -1 , shear heating is always important in all three lithospheric units (Fig. 9a). As expected, for the same velocity shear heating is always more intense in the upper and colder regions of the lithospheric units. For example, for a velocity difference of 3 cm.yr -1 the expected temperature increase in a lower crust made of ca. 400 o C hot anorthite is between 100 and 150 o C (Fig. 9a).

Shear heating is even more important if we consider a thinned thermally relaxed continental lithosphere, for example, at a passive continental margin. This is because the temperatures at the top of the lower crust and mantle lithosphere are colder than for a normal continental lithosphere (Fig. 9b).

We also analyze the initial stresses for configurations for which shear heating and strain localization is significant (Fig. 9c andd). We consider scenarios for which ∆T > 50 o C and for which initial shear stresses, τ xy , are < 1 GPa (Fig. 9c andd). The initial shear stresses are the largest stresses during shear zone formation since stress magnitudes decrease during shear zone formation due to thermal softening (Fig. 2c). We assume a velocity difference of 3 cm.yr -1 and vary initial bulk strain rates by varying the 1D model size, L. For flow laws of wet anorthite shear heating is significant for ambient temperatures between 380 and 470 o C and for strain rates, ˙ , between 10 -16 and 10 -13 s -1 . For example, for typical tectonic strain rates ˙ = 10 -15 s -1 shear localization by thermal softening is significant for shear stresses between 200 and 400 MPa for ambient temperature between 420 and 470 o C (Fig. 9c). For dry olivine and for ˙ = 10 -15 s -1 thermal softening is significant for shear stresses between 200 and 400 MPa for ambient temperature between 540 and 570 o C (Fig. 9d).

Our results indicate that the shortening of a tectonic plate, for example around a thinned passive continental margin, can likely generate significant shear heating and associated spontaneous shear zone generation by thermal softening. Such shear zone generation can take place in the ductile regime without reaching a brittle-plastic yield stress. Ductile strain localization by thermal softening could cause the generation of subduction zones, which is supported by numerical simulations (e.g. Thielmann and Kaus 2012).

Discussion

Localization criterion and Brinkman number

We already suggested the use of T max -T 0 > 50 o C or Br 1 > e 2 as localization criterion. Different versions of the Brinkman number (e.g. [START_REF] Brinkman | Heat effects in capillary flow i[END_REF][START_REF] Gruntfest | Thermal feedback in liquid flow; plane shear at constant stress[END_REF][START_REF] Yuen | Shear deformation zones along major transform faults and subducting slabs[END_REF][START_REF] Brun | Strain heating and thermal softening in continental shear zones: a review[END_REF] have been proposed and are also known under different names, for example, Gruntfest number (Gr).

Using the relations ∆v = ˙ 0 L and τ 0 = µ 0 ˙ 0 , several Brinkman numbers can be formulated:

Br 2 = ∆v 2 µ 0 λT 0 = ∆v 2 µ 0 λT 0 L 2 L 2 = ˙ 2 0 µ 0 L 2 λT 0 = ˙ 0 τ 0 L 2 λT 0 = τ 2 0 L 2 µ 0 λT 0 = Gr (27)
The right-most version with the square of the stress is often termed Gruntfest number, Gr. All parameters with the subscript 0 are initial, bulk values before the occurrence of strain localization or shear zone formation. The parameter L is the model size and not the thickness of the shear zone. The same exercises can be repeated by dividing equation ( 27) with the dimensionless Arrhenius term Q/RT 0 , and it would result in several forms of Br 1 (equation 24). A particular localization criterion would be

Q RT 0 L 2 ˙ 0 τ 0 λT 0 > e 2 , (28)
which is identical to the criterion of [START_REF] Karato | Deformation of Earth Materials: An Introduction to the Rheology of Solid Earth by Shun-ichiro Karato[END_REF], if e 2 on the right hand side is replaced by 1. The different versions of the Brinkman number are useful for different deformation scenarios. For example, if the deformation is driven by an applied shear stress, then the version with the square of the stress, i.e. Gr, is useful. If the thickness of the shear zone is a priori defined by the model or experimental configuration, then a version including L is useful whereby L then represents the pre-defined shear zone thickness. In general, for kinematically driven models, we prefer versions without any length scale L, because the model size does not affect the shear zone evolution (Fig. 2a and2b).

Our results show the applicability of three different localization criteria. As example, we use Br 2 > 1 (Fig. 6c). Applying a typical plate tectonic velocity of 3 cm.yr -1 , an effective viscosity of 2 × 10 23 Pa.s, a thermal conductivity of 3 W.m -1 K -1 and an ambient temperature of 500 o C (773 K) yields Br 2 ≈ 78. Our results show that for this value of Br 2 shear zone generation by thermal softening can occur (Fig. 6c). Using a typical tectonic strain rate of 10 -15 s -1 , and the applied effective viscosity of 2×10 23 Pa.s generates a shear stress of 400 MPa, which is a feasible flow stress for the upper and colder regions of the mantle lithosphere or the lower crust. Indeed, the spontaneous generation of km-scale shear zones by thermal softening was demonstrated in 2D thermo-mechanical simulations of lithospheric shortening for viscoelastoplastic rheology (Jaquet et al., 2017;Jaquet and Schmalholz, 2017).

Here, we do not consider viscoelastic effects. However, it was shown that elasticity can significantly impact thermally-induced strain localisation (Regenauer-Lieb and [START_REF] Regenauer-Lieb | Rapid conversion of elastic energy into plastic shear heating during incipient necking of the lithosphere[END_REF][START_REF] Duretz | Shear heating-induced strain localization across the scales[END_REF][START_REF] Jaquet | Dramatic effect of elasticity on thermal softening and strain localization during lithospheric shortening[END_REF] so that strain localization can be even more significant than predicted by our localization criteria.

Shear zone thickness

A ductile shear zone is commonly observable in the field, or experiment, by the significant variation of finite strain across the shear zone. The width of the varia-tion of finite strain across the shear zone is nearly one order of magnitude smaller than the corresponding width of the temperature variation (Fig. 4b). For significant shear heating the finite strain thickness, W γ , is ca. 2 √ κt/7 (Figs. 4b and5b).

We assume that the observable width of a shear zone is determined by values of finite strain > 1. Based on Fig. 4a this observable thickness is approximately 2 to 3 times thicker than W γ . Assuming a typical thermal diffusivity of 10 -6 m 2 .s -1 yields W γ = ca. 1.5 km and, hence, an observable thickness of 3 to 4.5 km for a shear zone which is active for 1 Ma. For a shear zone which is active for 4 Ma W γ = ca. 3 km and the observable thickness is 6 to 9 km. Since W γ only depends on time and is independent on the applied shear velocity it can be applied to any shear velocity and displacement. For example, a shear displacement of 100 km for a shear velocity of 2.5 cm.yr -1 requires 4 Ma, for which the observable thickness is 6 to 9 km. A shear zone with such thickness, velocity and displacement is likely typical for major lithospheric shear zones related to subduction zones.

Such thickness relation only applies to depth levels in the lithosphere for which thermal softening controls the strain localization. These durations of shear zone activity and corresponding predicted shear zone thicknesses agree also with those formed by thermal softening in 2D thermo-mechanical numerical simulations of lithospheric shortening (Jaquet and Schmalholz 2017;Jaquet et al. 2017). Based on the same relationship we expect sub-mm thickness for all shear zones that have been active for less than a few seconds.

If natural shear zones would have been formed by thermal softening with moderate temperature increase of 75 to 150 o C, then there would be no significant temperature variation between the shear zone and its wall rock because natural shear zones are observable due to the significant finite strain variation. This difference between finite strain and temperature variation explains why many ductile shear zones do not exhibit a significant observable temperature variation. The lack of a sharp and observable temperature variation, for example expressed by variation in metamorphic grade, is not a sufficient argument against the importance of shear heating and thermal softening.

Thermal softening and grain size reduction

There is still ongoing dispute concerning the primary mechanism of ductile strain localization in the lithosphere. Alternative to thermal softening, grain size reduction in combination with mechanisms, such as pinning, that prohibit grain growth (generally referred to as damage) is often proposed as primary strain localization mechanism. Clearly, in nature both mechanisms act simultaneously. We argue that thermal softening is a suitable mechanism for spontaneous strain localization in essentially homogeneous material whereby only minor heterogeneities can trigger strain localization. Grain size reduction can assist thermal softening and grain size reduction is likely a mechanism that is important during progressive shear zone evolution and can decrease the widening rate of the finite strain profile due to heat conduction. For example, [START_REF] Thielmann | Intermediate-depth earthquake generation and shear zone formation caused by grain size reduction and shear heating[END_REF] studied numerically the formation of shear zones by thermal runaway using a combined approach of thermal softening and grain size reduction. They showed that grain size reduction reduces the stress required for thermal runaway and hence assists ductile shear zone formation by thermal softening. Currently, different grain size evolution models are applied, for example, [START_REF] Thielmann | Intermediate-depth earthquake generation and shear zone formation caused by grain size reduction and shear heating[END_REF] apply the so-called paleowattmeter model in which grain size is a function of flow stress, strain rate and temperature, whereas [START_REF] Platt | Influence of shear heating on microstructurally defined plate boundary shear zones[END_REF] applies a piezometer in which grain size depends on flow stress only. To reliably quantify the impact of grain size evolution better constrained grain size evolution models for various rock types are needed.

Conclusions

A ductile shear zone which is generated spontaneously by thermal softening during a velocity-driven bulk deformation exhibits the following fundamental features: (1) After a transient period of temperature increase the temperature in the shear zone remains constant for linear viscous flow and quasi-constant for powerlaw viscous flow.

(2) The shear stress in the shear zone is largest at the onset of shear zone formation and subsequently decreases towards a (quasi-)constant value associated with the establishment of a (quasi-)constant temperature.

(3) The width of temperature variation across the shear zone is 6 to 8 times wider than the variation of the corresponding finite strain. Therefore, the shear zone does not exhibit a sharp, and hence easily observable, temperature variation between highly-strained shear zone and little-strained wall rock. (4) The shear zone is continuously widening during shearing due to thermal conduction between shear zone and wall rock.

(5) Shear heating starts to dominate the heat transfer between shear zone and wall rock once the temperature increase in the shear zone is > ca. 50 o C.

Different versions of the Brinkman number can predict the onset of shear zone generation by thermal softening. However, the Brinkman number cannot quantify the temperature increase inside the shear zone and, hence, the intensity of thermal softening. We derived a new analytical formula that predicts the maximal temperature inside the shear zone. This temperature prediction requires only information on the bulk deformation, such as far-field velocity, flow law and thermal parameters, and, therefore, no a priori knowledge of the shear zone itself, such as thickness, flow stress and strain rate. Temperature predictions across the scales of geological velocities show first order agreement with several natural shear zones including Alpine basement nappes, eclogite shear zones and pseudotachylites. We show with 1D, 2D and 3D numerical simulations that this temperature prediction is valid for shear zone generation under both bulk simple and pure shear.

Our results indicate that shear zone generation by thermal softening likely occurs during lithosphere deformation in the continental lower crust and the mantle lithosphere for typical lithospheric velocities of few cm.yr -1 or bulk strain rates between 10 -16 and 10 -14 s -1 . For these deformation conditions, shear stresses of few hundred MPa can already cause shear zone generation by thermal softening.

Based on our results and their application to lithospheric flow laws and deformation conditions, we argue that spontaneous shear zone generation by thermal softening is a feasible and likely the primary mechanism for spontaneous lithospheric scale shear zone generation. Thermal softening is probably a key constituent of subduction initiation, for example, at a thinned passive continental margin. Additional processes, such as grain size reduction, fabric development or fluid-related reactions can cause additional softening during progressive shear zone evolution and likely intensify the strain localization.

Acknowledgements

Lithology

A [Pa -(n+r) s -1 ] n f H2O [Pa] r A eff [Pa -n s -1 ] Q [J.mol -1 ] λ [W.K -1 .m -1 ] ρ [kg.m -3 ]
Wet quartzite 1 6. Table 1: Rheological and thermal parameters for the used lithologies. A eff = F Af r H2O d -p is an effective pre-exponential factor ('A' in the main text) that incorporates grain size (d) and water fugacity (f H2O ) dependence. In all cases we use dislocation creep therefore the grain size exponent is p = 0. All of these flow laws describe stress and strain rate relationship in uniaxial compression experiments. In order to convert them into strain rate dependent invariant forms we need to introduce a geometry factor, which is F = 2 n-1 3 (n+1)/2 for all presented cases (see e.g. [START_REF] Gerya | Introduction to numerical geodynamic modelling[END_REF]).

The rest of the parameters are: power law exponent (n), water fugacity exponent (r), activation energy (Q), thermal conductivity (λ), density (ρ) and finally heat capacity is constant for all (c p = 1050 J.kg -1 K -1 ). The sources of the rheological parameters are: In all models a thermal perturbation (red) is in the model center. Due to the symmetry of the problem we solve the 2D and 3D models only for the positive coordinate region. To compare 1D results with 2D and 3D results, the results of the 1D model are rotated so that they correspond to the direction x in the 2D and 3D models which is orthogonal to the shear zone. 2) is based on a dry peridotite flow law (see Table 1). In both models we applied 3 cm.yr -1 velocity difference. Different lines correspond to models with different model size and initial perturbation (see legend, which applies to all three panels). After a transient stage 2b and2c. The temperature profile is significantly wider than the finite strain profile. The thickness of both profiles is measured at the vertical value which is half between the maximum and minimum value of the profile. b) Evolution of the ratio of thermal thickness to corresponding finite strain thickness with progressive time for different model configurations (colors are corresponding to Fig. 2b and2c). After a transient stage the ratios approach values between 6 and 8 showing that the temperature variation is nearly one order of magnitude wider than the corresponding finite strain variation. For comparison, also the ratio of thermal thickness to corresponding instantaneous strain rate thickness is displayed with gray dashed lines, because this thickness ratio is constant and ca. 4 for all shown simulations.

  1D model and the smallest for the 3D model. This is because in the 1D model the initial thermal perturbation is at the position of the future shear zone whereas in the 2D (Fig.7dto f) and 3D (Fig.7a to c) models the initial thermal perturbation is present only in a fraction of the future shear zones. Also, the background temperature increase due to bulk shear heating of the model domain is the largest in the 3D and the smallest in the 1D model (Fig.7g to i). Nevertheless, equation22(associated with the 1D results) accurately predicts the temperature inside the 2D and 3D shear zones after the transient stage of temperature increase. The results also confirm that the initial temperature perturbation applied in the 1D model has no impact on the maximum temperature in the shear zone because the temperature of the 2D and 3D shear zones are unaffected by the initial thermal perturbation. A comparison of profiles of the velocities parallel to the 1D, 2D and 3D shear zones

1

  Hirth et al. 2001, 2 Carter and Tsenn 1987,3 Rybacki and Dresen 2004, 4 Hirth and Kohlstedt 2003.

Figure 1 :

 1 Figure 1: Model configurations for 1D simple shear, and 2D and 3D pure shear bulk deformation.

Figure 2 :

 2 Figure 2: Representative results of 1D shear zone models. Time evolution of maximum temperature at the shear zone center for linear viscous (a) and power-law viscous (b) flow. Model (1) is based on a dry anorthite diffusion creep flow law with A eff = 0.14 Pa.s -1 (d = 0.3 mm) and Q = 467 kJ.mol -1 , model (2) is based on a dry peridotite flow law (see Table1). In both models we

  the maximum temperature converges to a constant (a) or (quasi-)constant (b) temperature. c) Shear stress evolution for power law simulations, shown in panel b). The shear stress always decreases with progressive shear zone evolution and converges to a quasi-constant value. evolution of different models

Figure 3 :

 3 Figure 3: a) Representative time evolution of a temperature profile during shear zone formation. The model setup is equivalent with the one of the ∆T 0 = 25 o C, and L = 1000 km from Fig. 2b and 2c. The dashed horizontal line indicates the thermal thickness which is measured at the temperature which is half between the maximum and minimum temperature of the corresponding profile. b) Evolution of thermal thickness for different representative simulations (same colors are used in Fig.2b and 2c). Each line shows the result of a simulation with different initial temperature perturbation and model size. After a transient period, the thickness evolution for all simulations is linearly increasing with the square root of time.

  ratio of temperature and strain profiles

Figure 4 :

 4 Figure 4: a) Temperature and corresponding finite strain profile for a simulation with a dry olivine flow law after 3.5 Myr. The model setup is equivalent with the one of the ∆T 0 = 25 o C, and L = 1000 km from Fig. 2b and 2c. The temperature profile is significantly wider than the finite strain

Figure 6

 6 Figure 5: a) Representative temperature evolution for cooling of an initial Gaussian temperature profile, left curves, and for heating for a constant temperature in the model center, right curves. Lines for t 1 to t 4 display temperature profiles at progressive times. The dashed horizontal lines indicate the half-width of a specific temperature profile. Horizontal x-coordinates and temperatures are dimensionless and temperatures have been scaled so that the initial temperature is identical. b) Plot of the dimensionless widths of temperature profiles determined from 1D numerical simulations versus the corresponding temperature increase in the shear zone centre, ∆T max . The blue and red horizontal lines indicate the theoretical dimensionless width for cooling of an initial Gaussian temperature profile (equation 13) and for heating for a constant temperature in the shear zone (equation 16), respectively. A value of ∆T max ≈ 40 o C indicates the transition between the two types of heat transfer and for ∆T max > 40 o C shear heating (SH) is significant.

Figure 7 :

 7 Figure 7: Comparison of temperature and velocity profiles across 1D, 2D and 3D shear zones. a) to c) shows the temperature field for three stages of 3D shear zone formation. d) to f) shows the temperature field for three stages of 2D shear zone formation. The profile lines indicated with x in both 2D and 3D models are used for the comparison with the 1D model. g) to i) show 1D temperature profile and the 2D and 3D profiles along the x profile for three different times. j) to l) show the 1D velocity profile and the 2D and 3D profiles along the x -profile for three different times. The displayed velocity magnitudes are normal to the profile orientation and, hence, parallel to the shear zone orientation.

  duction (right side of Eqn. 18). Because the quasi steady state occurs only in the shear zone center, it is reasonable to chose characteristic values that are representative for this location. There are several formally correct and reasonable choices for the characteristic values, but after testing the scaling relationship (19) with the numerical results, we found that T c = RT 2

max /Q, v c = ∆v and µ c = µ min provides the best fit. T max and µ min are always the maximal temperature and minimum viscosity, respectively, in the shear zone center. Based on equation (

19

) and these characteristic values, T max can be predicted with:

We thank an anonymous reviewer for constructive and helpful comments. This work was supported by SNF grant No. 200020-149380 and the University of Lausanne. We thank Ludovic Räss for assistance in performing the 3D numerical simulations.

 Table 1). The same constant finite shear strain, γ c = 20, is assumed for all shear zones and the corresponding duration of the deformation is calculated with equation 26. For pseudotachylites the typical values of finite strain are higher. The expected temperature range for γ c = 80 is indicated by the light gray area. The quartzite flow laws are not displayed for high velocities because for those the argument of the logarithm is approaching exp(-1.1) ≈ 0.3, where the prediction starts to significantly deviate from the solution.


domain: shs shs.archi spi

text:
David Lo Buglio
email: david.lobuglio@gamsau.map.archi.fr

Vanessa Lardinois
email: vanessa.lardinois@gamsau.map.archi.fr

Livio De Luca
email: livio.deluca@gamsau.map.archi.fr

Revealing shape semantics from morphological similarities of a collection of architectural elements The case study of the columns of Saint-Michel de Cuxa

Keywords: Architectural representation, morphological analysis, semantic characterization, image Processing

Over the last three decades, the introduction of digital technologies in the field of architectural documentation has profoundly changed the tools and the acquisition techniques. Most of the developments concerns metrical and colorimetric characteristics of the objects studied. These developments, surrounding the practice of architectural survey, tend to respond primarily to the requirements of completeness. In this context, it seems necessary to assess the impact of these instruments on the cognitive value of architectural representation. With a strong technological presence, the study of the built heritage is facing a problem of "overload" of information. Indeed, it fails to strengthen the representation in its role as vehicle of knowledge. Confronted with the intelligibility deficit, this article propose an original approach for reading morphological features of an artifact by using a bottom-up approach: the meaning of the elements comes from the analysis of low-level geometric properties of a collection of instances related to the same theorical model. The idea is to rely on data accumulation in order to make emerge high-level semantic features from the comparative analysis of common low-level geometric features. The introduced principles are illustrated by the comparative analysis of 31 columns of the cloister of the abbey of Saint-Michel Cuxa, aimed to identify a common semantic layout.

I. INTRODUCTION

Today, it seems that the three-dimensional representation is an effective medium for the documentation, the study and the heritage enhancement. In this context many research projects focus on tools and methods of architectural digitization. But, beyond the technological impact, the understanding of their contribution for the production and the cognitive value of representations constitute a major challenge for research.

Over the past three decades, the introduction of digital technologies in the field of architectural documentation has profoundly changed the instruments and techniques of acquisition. These developments have increased the accuracy and completeness of architectural survey. The discrimination of information provided by the operator doesn't seem to be decisive anymore in a survey process based on quantitative requirements [START_REF] Collectif | Le Relevé En Architecture Ou L'éternelle Quête Du Vrai. Lieux Dits[END_REF]. Whether it is the lasergrammetry or the photogrammetry, all the techniques being developed for heritage documentation tend to respond primarily to the completeness requirements.

If there is a desire to use the digital model for structuring, organizing and archiving heterogeneous information, the threedimensional representations of buildings are mostly exploited as a means to document the existing state of an object. By continually pushing the level of completeness of the digital models, nowadays one can observe a strong metric and visual coherence between architectural artifacts surveyed and their digital instances.

However, the growing body of point clouds, 3D models or un-interpreted data does not seem to increase the level of intelligibility or the cognitive dimension of the representations produced. For example, human analysis is still crucial to the understanding of the different conceptual articulations of an architectural element. Faced with this overload of "information" and the intelligibility deficit resulting, it is necessary to provide a methodological answer to articulate, within the same graphical space, the concepts of "multiplicity" and "intelligibility".

Today, the majority of responses to this issue are based on a "top-down" approach, i.e. an approach that uses the prestructured knowledge of the domain [START_REF] Treisman | A Feature-integration Theory of Attention[END_REF]. Although there are several support tools for the organization, the segmentation and the semantic characterization of 3D data, most of these solutions depend strongly on the level of the operator's understanding of the object [START_REF] Spina | Point Cloud Segmentation for Cultural Heritage Sites[END_REF]. One of the challenges of this research is to provide an alternative to the traditional approaches (subjective) with the aim to offer (objective) analysis solutions to increase the cognitive value of representations. So, this paper explores how the statistical study and the data accumulation allow achieving an objective estimation of the morphological properties of an architectural object. In line with this ambition, the second challenge is to evaluate how the creation of a morphological signature (common to a set of elements) can be used to the shape classification. Through the study of the columns of the cloister of the abbey of Saint-Michel Cuxa, we will assess what the examination of morphological similarities can bring to the semantic characterization of architectural elements.

By focusing on the epistemological issues surrounding the survey and the representation of architectural artifacts, this article tries to:

• locate issues related to the intelligibility deficit.

• mention some approaches implemented to answer at this first issue. It should be noted that these approaches focus mainly on a description "from above" (topdown), ie a description that takes advantage of the prestructured knowledge of the domain.

• propose an alternative to these approaches by focusing on the possibilities of a "bottom-up" study which relies on the "objective" estimation of the morphological characteristics of an object.

• finally, evaluate the relevance of this new approach from the results obtained on a specific case study: the columns of the cloister of the Abbey of Saint-Michel Cuxa.

II. BETWEEN COMPLEXITY AND INTELLIGIBILITY OF ARCHITECTURAL REPRESENTATIONS

With their ability to combine within the same graphical space the architectural codification and the figuration of the object, the canonical representation systems of architecture (plan, section, elevation and perspective) are, since the Italian Renaissance, the privileged vehicle of the shape knowledge. Nevertheless, this paradigm of communication, common to all disciplines concerned with the study of built heritage, is implicitly challenged by the arrival of digital technologies [START_REF] Mediati | L'occhio sul mondo per una semiotica del punto di vista[END_REF]. Indeed, many tools and methods in the field of architectural documentation have been modified by the introduction of digital solutions. The architectural photogrammetry has reached a level where access and cost are able to fill many needs at the scale of the built heritage. The development of multi-stereo correlation systems allows registration of dimensional and colorimetrical properties. Therefore, It is no longer uncommon to get scans describing the architectural object with millions of colored points (coordinates) [START_REF] Pierrot-Deseilligny | Automated Image-Based Procedures for Accurate Artifacts 3D Modeling and Orthoimage[END_REF].

If this technology craze has helped to strengthen the architectural survey, it did not consolidate the intrinsic cognitive dimension to its practice. Its analytical value is mainly due to the work of interpretation and at the choice of information collected on the object. In a context where the level of visual and metric coherence between the object and its digital instance continues to grow, the mass of un-interpreted data raises questions. Is the multiplicity of acquired data is able to increase the level of intelligibility of architectural representation?

A study published in 2011, on the cognitive dimension of architectural surveying, offers an observation of the intelligibility level of a set of representations depending on the degree of automation of the tools and methods mobilized [START_REF] Buglio | Critical Review of 3D Digitization Methods and Techniques Applied to the Field of Architectural Heritage: Methodological and Cognitive Issues[END_REF]. We mean with "intelligibility", the ability of the 3D digitization to refer to a world of knowledge. One understands that the informative richness of digitization is partly linked to the cognitive engagement of the operator and not the amount of data they contain. In order to achieve representations with a high degree of intelligibility, the action of the operator seems crucial for the analysis, classification and interpretation of data.

Current approaches (or automated process) effectively guarantee multiple data processing (color and geometry) without improving the intelligibility of representations produced. In fact, they respond only partially to the needs of communities concerned with the morphological analysis of the building heritage [START_REF] Blaise | Informative Modelling[END_REF]. In this context, the only research on the development of tools for the acquisition and data processing cannot respond structurally to the intelligibility deficit observed. It should be noted that this criticism is not intended to call into question the relevance of technological advances in the field of architectural representation but to understand their limits and deficits.

The main issue of this work is whether the introduction of a new methodology for shape analysis and statistical analysis would permit an approach where the data flow would support the intelligibility of representations produced. Thus, this paper presents an approach to complete the acquisition protocols. This research, based on the use of multiple data, aims to produce a semantic description from a low level morphological analysis.

III. FROM THE HIGH LEVEL SEMANTIC DESCRIPTION TO THE LOW LEVEL GEOMETRIC DESCRIPTION

For some years now, the communities involved in the heritage documentation seem gradually take over the problems related to intelligibility deficit. There are now numerous studies concerned with the exploitation of non-interpreted data. Indeed, the geometric description and the analysis of architectural shapes require a reduction of multiplicity to intelligibility. Nowadays, there are several applications for the understanding of geometric shapes. Research around the automatic and semi-automatic segmentation, the semantic annotation, the parametric modeling or the image-based modeling from geometric primitives is symptomatic of the need to organize the 3D digitization around structured knowledge [START_REF] Oxholm | Reassembling Thin Artifacts of Unknown Geometry[END_REF] [9] [START_REF] Luca | La photomodélisation architecturale. Relevé, modélisation et représentation d'édifices à partir de photographies[END_REF].

Thus, some recent research provides satisfactory answers to the problems of intelligibility. Most of them can be considered as a "top-down" approach because the operator leads the analysis from the top by taking advantages of the prestructured knowledge of the domain. This approach requires a semantic enrichment process where human analysis is decisive for the steps of segmentation, 3D reconstruction and / or semantic annotation.

It was also at the heart of two relevant research approaches. The first one focused on the formalization of moldings based on the rules governing the combination of geometric atoms [START_REF] De Luca | A Generic Formalism for the Semantic Modeling and Representation of Architectural Elements[END_REF]. This is a formalization of high level since the rules emanate from an architectural semantics. The second one focused on the development of a parametric modeling tool (Generative Modeling Language) for the formalization of architectural concepts from a low-level geometric vocabulary. In this case, it is a lower level since the reconstruction articulates a geometric semantics [START_REF] Havemann | Generative Parametric Design of Gothic Window Tracery[END_REF].

If these two examples differ in their analysis scale, architectural semantics on one side and geometric semantic on the other side, they both rely on a grammar of shapes. Therefore, we can consider that the two works offer a "topdown" approach. But beyond this first orientation, it seems necessary to complete this first step with an approach that would take advantage of the multiplicity of acquired data. This implies that the meaning of the elements would come from a "free" morphological estimation and not from a shape thesaurus. The "free" estimation corresponds to the opportunity to study without bias on the results, the point aggregate, the curvatures of surfaces and volumes of the architectural elements.

Beyond the study of singular elements, the use of statistical analysis would permit to find correspondences between different geometric entities. Morphological comparisons that can be done on a set of shapes offer opportunities to study similarities in a collection of architectural elements. Indeed, by establishing common geometric observation criteria and multiplying comparisons between elements, it is possible to identify morphological signatures that can contribute, in a second step, at the semantic characterization of a collection. This approach is therefore based on the self-expansion of a large amount of data.

The medium-term challenges of this approach are:

1. to establish new shape classification criteria. The comparison of architectural objects on the basis of geometric similarities aims to create libraries of architectural shapes organized around pre-structured knowledge.

2. to question architectural treaties with the classification work (through the geometrical analysis of the typical shapes, the composition rules, the principles of scale, the positioning constraints, the orientation, etc...).

3. to study the variability of architectural elements. In the field of built heritage, the morphological signature (or theoretical model) of a collection of elements can serve as the basis for the study of their variability in space and time. In addition, it can also have an impact on the understanding of the evolution of a style. In the continuation, the use of self-extension mechanisms can provide useful answers to observe unclassified and uninterpreted shapes.

IV. LOW LEVEL DESCRIPTION : THE CASE OF THE COLUMNS OF SAINT-MICHEL DE CUXA In order to experiment this approach, a first study was conducted on the columns of the cloister of the abbey of Saint-Michel Cuxa. They have the dual advantage of being in a good state of conservation and belong to a group composed of 31 columns present on site. Indeed, the "bottom-up" approach is only relevant if it is part of a cumulative study. Therefore the challenge is to assess the common denominators / geometric similarities in the columns collection. It is also interesting to note that the many theories that have existed since the Renaissance make the column a relevant case study [START_REF] Palladio | The Four Books of Architecture[END_REF].

The identification of a morphological signature could lead us towards a common theoretical model to all of the columns of the cloister. It could emerge a 32nd theoretical column from the 31 columns that still exist. It may be recalled that this approach aims to assist human analysis on the shapes classification and comparison.

The survey (used in this study) combines the latest representation systems with the architectural image-based modeling developments. Indeed, at the level of architectural elements documentation, image-based modeling (survey, 3D modeling and representation from images) has reached a point where costs, accessibility and quality of results can fill many needs at the built scale. The columns have been digitized during a photogrammetric survey campaign in December 2010. The digital instances of each column (obtained after the acquisition on site and the processing of the different point clouds) express a certain level of completeness "Fig. 1". But beyond a high geometric fidelity, our concern is the degree of intelligibility of the digitization, and their ability to build a theoretical model.

Before describing the technical principles of the approach (and the establishment of a first morphological signature), we will observe the tools currently available at the level of shape analysis. Thus, this chapter will focus on:

• the conversion of 3D digitization into a twodimensional representation (most suitable for image analysis).

• the existing algorithms in the field of image analysis and thus on the relevant shape descriptors for the creation of morphological signature.

• the study of the major geometric discontinuities of each column, the evaluation of their average and their re-projection in the form of relevant profiles (in order to build a first morphological signature).

• the link between each digital instance and its theoretical model.

A. Acquisition and extraction of depth maps

In order to use existing tools in the field of image analysis, one of the issues is to convert in two dimensions the complex geometric articulations of the columns. The polygonal models of each column are projected onto a two-dimensional depth map. This map is obtained from the projection of the polygonal mesh on a UV cylinder surrounding it. The use of a cylinder is not intended to characterize the geometric complexity of the object. The ambition is to refer to the concept of the column by exploiting its geometric archetype "Fig. 2".

The interest of the approach is that each pixel of the cylinder is encoded as a grayscale value determined from the distance separating it from the polygonal mesh. The unfolded projection cylinder (or the two-dimensional depth map) contains all the geometric complexity of the digital instance. From this map or from other computed maps (like the normal map and the ambient occlusion map), it is therefore possible to lead a precise morphological study based on image analysis tools [START_REF] Papadakis | Panorama: A 3D Shape Descriptor Based on Panoramic Views for Unsupervised 3D Object Retrieval[END_REF].

B. Morphological analysis and description of low level

At the level of shapes study, image analysis algorithms offer many opportunities for contour extraction, identification of lines and analysis of similar regions within an image [START_REF] Zhang | Review of Shape Representation Anddescription Techniques[END_REF]. To understand the benefits and the limitations of these tools for the low level morphological study, all the known solutions have been systematically experimented and rigorously inventoried.

If, some algorithms are an appropriate solution for image analysis in the industrial or the medical field, their diversion for the study of heritage artifacts presents some difficulties. They require a hard work of setting and adaptation to meet specificities of built (automatic identification of rhythms, geometric discontinuities, similarities,...).

A first experiment has permitted to assess some shape analysis operators directly on the depth maps coming from the columns "Fig. 3". It is important to note that the operators are estimated on their contribution to the low level geometric description, e.g. on their ability to identify some geometric breaks in a column. The three families of description methods and image process methods mobilized were:

• Threshold-based segmentation • Contour-based segmentation • Region-based segmentation
Nevertheless, the analysis operators only make sense if they can isolate, bring together or discriminate some of information within the image. We mean by relevant descriptors, those able to locate specific morphological characteristic (scale variation, symmetry, degree of curvature, surface discontinuity,...).

Before going further, it is important to remember that the description of the vast majority of shapes can be done by the knowledge of: 1. the geometry, which refers to the understanding of volumes and surfaces of the object 2. the topology, which focuses on neighborhood relations within a shape 3. and finally, the organization or the spatial composition of the basic physical entities

These three points are an appropriate observation triangle at different levels of granularity (scale). Since our interest is on the entity as a whole (column), this work focuses primarily on the geometric level. This level involves the characterization of the morphology of an element according to two aspects: its curves and its discontinuities.

Indeed, the geometric study of a shape implies the knowledge of its curves. This is explained by the fact that any surface can be defined by "n" curves. For each one corresponds a unique profile whose mathematical representation is similar to a set of functions. The calculation of curvature maps is therefore an essential aspect for the understanding the surfaces [START_REF] Floriani | Shape Analysis and Structuring[END_REF].

But beyond this first aspect, a shape can also be described by its major discontinuities. Is meant by major discontinuities, each surface change marked by the rupture of a curvature. In the particular case of the study of architectural elements, the geometrical discontinuities provide the understanding of the conceptual articulations present in the object. Each passage of a concept to another (e.g. from the shaft of the column up to the astragalus) results in a strong discontinuity of surface. Considering the interest that can represent this last aspect for the communities involved in the study of the built heritage, the rest of this chapter is primarily dedicated to the results obtained through the study of geometric discontinuities.

If this paper focuses, up to now, on the study of singular elements, it is necessary to remember that the final goal of this research is to achieve a common morphological signature to the 31 columns.

C. Study of morphological discontinuities and 3D reprojection of the average profiles

To summarize, the general issue of this work is to provide a response to the lack of intelligibility mentioned before, a deficit that the technological developments do not solve. Rather than consider the enrichment of representations by a "top-down" approach, the method proposed here is based on the accumulation of data. Through the accumulation, the idea is to bring out common geometric characteristics of to a set of architectural elements. In this case, our attention focuses on the study and on the exploitation of discontinuities with the aim to establish a first common morphological signature of the columns. Its implementation is organized around a development in three stages:

1) First time: line detection

The observed image analysis tools provide useful technical tracks for the implementation of shape descriptors. Indeed, the study of discontinuities requires to locate and isolate significant lines present in the depth maps. They correspond to abrupt changes between shades of gray on a portion of the image.

Among the operators for the contour-based segmentation, the standard Hough transform (SHT) is a pattern recognition technique to automatically detect the horizontal and vertical lines present in the image.

In order to optimize the detection operator, it is preferable, in a first step, to apply a contour detection like the type Sobel [START_REF] Jähne | Principles of Filter Design[END_REF]. It provides a much lighter image with less detail. The contours are sharper and the shapes are properly closed. The threshold value is adapted to obtain an optimal display of contours. This pre-treatment is necessary because the Hough operator is preferably applied on binary images.

At the level of the lines detection operator, the general principle of the Standard Hough Transform is to admit the existence of an infinite number of lines passing through all the points of the image. Each segment can be characterized with two parameters: an angle θ (Theta) and a length ρ (rho) perpendicular to the line, and passing through the origin of the system. The lengths and the angles obtained in the Cartesian system (x, y) are then represented in the form of coordinate vectors RHO-THETA in Hough space. They are (in this space) sinusoidal curves that correspond to each point of the image. The intersection of two curves marks the linear link between two points in the Cartesian system [START_REF] Duda | Use of the Hough Transformation to Detect Lines and Curves in Pictures[END_REF].

Due to its use on architectural elements, sometimes strongly altered, it may happen that the Hough transform fails to correctly identify some lines of the image. Although the majority of the line segments detected are horizontal, the transform also recognizes oblique lines (not adapted for future use). Solutions of discrimination have been created in order to correct these errors of interpretation of the Hough transform. Some oblique (whose inclination does not exceed a discriminator value) are adjusted to obtain usable horizontal, i.e. a unique coordinate "y" for the discontinuity. This adaptation is necessary to meet the requirements of the morphological study of the column.

As presented, the Standard Hough Transform has been adapted and configured in order to automatically bring out all the linear networks present in the depth maps. These networks correspond to the visible representation of the major geometric discontinuities (both horizontal and vertical) of the object. The analysis of low-level (or bottom-up), as seen here, seems to provide a satisfactory response to the automatic detection of conceptual articulations present in the column.

The observation of the results shows that all the discontinuities detected actually correspond to the transition between two architectural concepts. One recognizes in the illustrations below "Fig. 4", the transitions between:

1. the plinth and the lower torus Before studying the average value whose aim is the creation of a morphological signature common to all the columns, this point focuses on the exploitation of the discontinuities at the level of the object. Although it is possible at this stage to conduct an analysis (for each depth map) on the heights and the distances between the discontinuities, it is also possible to achieve the simplified geometric reconstruction of each 31 columns. Indeed, each discontinuity on the map corresponds to the position of a relevant profile for the simplified formalization of the object.

These profiles can be obtained by the projection of the gray values of each pixel located along the discontinuity. By considering that these shades are the graphical representation of the distances separating the digital instance from the cylinder of projection (see above), this implies that all the pixels contain the information of their position. These values (coded from 0 to 255) are first converted into a metric scale and then projected perpendicularly to the plane of the depth map (here along the Z axis) "Fig. 5".

To find similar profiles of a column section, it is necessary to transfer the result of the "planar" projection on a "cylindrical" projection. To solve this problem, the profiles obtained in a Cartesian coordinate system are transposed in a polar coordinate system where the value of each point in Z corresponds to the radius (ρ) of the new profile and the value X to the radians (θ).

With this proposal, we can already consider the analysis of profiles and discontinuities as a relevant way to describe certain geometric characteristics of the column (e.g. rhythm, minimums, maximums and average).

3) Third step: study of average values

If the previous step provides interesting avenues for the reconstruction of simplified geometries, the study of average offers opportunities for the creation of morphological signatures. It is conceivable that the "hybridization" of the values of discontinuities (from a collection of columns) can bring out the generic characteristics of the shape. In this perspective, this point focuses on the creation of a morphological signature based on the geometric discontinuities of the 31 columns processed.

One might think that the average calculation of all the depth maps of the 31 columns may be sufficient to formalize a valid general average. However, such an approach would lead to the construction of a 32nd theoretical column shapeless, unable to guarantee the intelligibility of its various articulations.

To bypass the geometric aberration that would result, it is possible to separate the study of average in three distinct parts:

1. the study of the average number of discontinuities observed for each depth map, 2. the study of the average vertical position of discontinuities for each of the 31 columns, 3. the study of the average gray values between discontinuities of the same level (i.e. between discontinuities in the same place on different maps). These different steps in the analysis have been implemented in a tool that is based on the developments of the previous point "Fig. 6". The mean values (in the form of profiles) are projected in 3D in order to represent graphically a first morphological signature. It is interesting to note that the theoretical model developed expresses all conceptual articulations present in the majority of the columns. The observation of all the columns shows that some of them have altered parts. For example, the column n°24 of the west wing has a heavily damaged base (especially at the level of the plinth). But with the accumulation of data, the degradation (even important) does not appear in the theoretical model. Although it is only the representation of the average discontinuities, the morphological signature described a first common semantic skeleton for all the columns, the 32nd theoretical column.

V. CONCLUSIONS AND PERSPECTIVES

In a context where the architectural surveying is dependent on tools and technologies that it mobilizes, the heritage documentation is facing a problem of "overload" of information. Once collected, they fail to strengthen the representation as a vehicle of knowledge. In front of the intelligibility deficit, this research has provided methodological orientation based on a "bottom-up" approach. The semantic characterization based on morphological similarities of elements offers new opportunities for the communities concerned with the study of built heritage. This work also shows that the enrichment of representations with prestructured knowledge of the domain could be complemented by the presence of a geometric description based on the study of similarities.

As we could see, the depth maps are an ideal medium for the analysis of discontinuities in a column. By their ability to highlight the conceptual articulations, these breaks help to ensure the reconstruction of simplified geometries, geometries that maintain their level of intelligibility. But beyond these aspects, one also observed that the calculation of averages constitutes an interest for the definition of morphological signatures. It opens avenues for the formalization of a 32nd theoretical column which interest is to serve as a reference for the study of altered elements.

However, one of the current limitations of this research is the lack of shape descriptors taken into account for the creation of this first morphological signature. In the current state of this work, the descriptors focus exclusively on the study of discontinuities. An answer would be to include other shape attributes. For example, one can think of the study of lines of curvature on surfaces.

In conclusion, the possibility of linking the concepts of multiplicity and intelligibility opens many tracks. If the simplified reconstruction of complex geometries shows interesting avenues for the field of 3D modeling, the formalization of morphological signatures constitute also an important issue. Their use could meet the following challenges:

• the visualization of the gap between the digital instance and the theoretical model (in order to better understand the alterations of buildings).

• the study of similarities between objects. The idea is to compare their discontinuities to find geometric correspondences in a collection of blocks, for example. This would permit to support the research in archeology, especially for the anastylosis work. • the classification of shapes according to low level morphological criteria. The use of theoretical models could provide an opportunity to better understand the propagation of an architectural style in time and space. This last point makes legitimate the reformulation of shape thesaurus in the light of morphological signatures.

Finally, to conclude, this approach avoids the failings of the representation of the territory 1:1 scale that Borges evoked [START_REF] Borges | De la rigueur de la science[END_REF], where the practice of survey is reduced to the faithful reproduction of the visual attributes of space.

Fig. 1 .

 1 Fig. 1. Masks and point clouds from the photogrammetric acquisition of the columns of Saint-Michel de Cuxa

Fig. 2 .

 2 Fig. 2. Generating a depth map from the polygonal mesh, analysis of geometric discontinuities and profiles from the map.

Fig. 3 .

 3 Fig. 3. Threshold-based segmentation and contour-based segmentation on the depth map from the west column No. 25.

2 .

 2 the lower and upper torus 3. the upper torus and the shaft 4. the shaft and the astragalus 5. the astragalus and the historiated capitals 6. the historiated capitals and the echinus 7. the echinus and the abacus 2) Second step: extraction and 3D reprojection of the relevant profiles

Fig. 4 .

 4 Fig. 4. Depth map, sinusoidal of the Hough transform, detection of image lines, extrapolation of lines and discrimination of certain oblique, view of the mean value (on "y").

Fig. 5 .

 5 Fig. 5. 3D reprojection of the discontinuities.