Sensitivity analysis for crustal rheological profiles: examples from the Aegean Region

We tested the validity of rheological profiles and their sensitivity to variations in the input parameters, with particular emphasis on the resulting BDT (brittle-ductile transition) depth and corresponding strength and temperature. For this purpose, we selected two test-sites from the Aegean Region, one in an extensional tectonic setting and the other in a strike-slip regime, and carefully realized the corresponding “reference” rheological profiles based on literature data and specific geological constraints. The reference envelopes have been then compared with a set of different profiles realized by varying the input parameters of the constitutive equations of the brittle and ductile behaviours within reasonable ranges. Firstly, tests were performed by changing the value of only one input parameter per time, with the aim of quantifying and comparing its influence on the BDT properties. The parameters exerting the greatest control are the activation energy, the power-law exponent and the surface heat flow (through its influence on the geothermal gradient), for the creep behaviour. As regards the brittle behaviour parameters, the friction coefficient and the pore fluid pressure could play a significant role especially in determining the maximum strength. In a second phase, we simultaneously varied all the input parameters in order to consider the possible synergistic effects on the resulting rheological profiles and to verify the likelihood and consistency of the reference models. For the statistical approach, one hundred thousand random combinations of the analysed parameters have been generated. The particular care spent on selecting the range of values of each parameter is reflected in the results of the statistical analyses, which show a good agreement with the reference profiles and allow estimating the overall uncertainties. Finally, the obtained strength envelopes have been compared with the accurately relocated depth distribution of recent seismic sequences that affected the two test areas. In both cases, the depth corresponding to the 95% of the total released energy nicely fits the BDT depth obtained from the rheological modelling, therefore confirming that this parameter could represent a reasonable and reliable approximation of the seismogenic layer thickness.

As concerns the ductile behaviour, it is classically modelled as a power-law creep [Sibson, 1977;Ranalli and Murphy, 1987;Kohlstedt et al., 1995; among others]: (2) where έ represents the strain rate (s -1 ), n is the exponent of the power-law, A is the so called "power-law parameter" (Pa -n ·s -1 ), E is the activation energy (J·mol -1 ), R the Boltzmann constant and T is the temperature (K), which is a function of depth (z). We will come back in a subsequent section on the law governing the thermal vertical gradient (geotherm), as far as its estimate represents a crucial issue of the methodological approach.
The transition from the brittle to the ductile behaviour at the BDT depth can hardly be considered as a sharp change from one prevailing rheology to another. Indeed, this rheological change is more likely to occur along a transitional zone where mixed behaviour is expected to characterize in space and time the rock volume [Sibson, 1977;Kohlstedt et al., 1995 among others]. Nevertheless, the simplified approach described by equations (1) and (2) could be regarded as an acceptable approximation at least for regional and continental scale studies. This is well documented by its extensive application in the literature [see, among others, Brace and Kohlstedt, 1980;Ranalli and Murphy, 1987;Cloetingh and Burov, 1996].

Rock density
Assuming as a reasonable first approximation that gravity acceleration is constant within the depth range of the brittle behaviour, only rock density has been tested. It primarily depends on the lithological composition of the overlying rocks. Major changes occur at the boundary between principal crustal layers, though minor lithologicalcompositional vertical variations could occur due to temperature-and compaction-induced changes. By including also possible sources of error and/or uncertainty associated to laboratory measurements, a possible variability up to ±10% with respect to a reference density has been considered in the numerical modelling.

Tectonic parameter
The parameter (equation (1)) depends on both tectonic regime ( ) and friction coefficient ( ); we thus analyse the two factors separately. The former has the effect of reducing and increasing the strength in extensional and contractional regimes, respectively [e.g. Sibson, 1974;1977]. In a purely Andersonian regime =0 indicates a uniaxial compression, while the value of =1 is associated with uniaxial extension.
For most common crustal rocks, the friction coefficient μ has been investigated since the 1950s [see the classic reviews by Byerlee, 1968;. He observed that the maximum friction coefficient linearly correlates the normal stress to the shear strength, independent of rock type. The statistical fit is quite good and suggests a progressive decrease of μ from ~0.85 to 0.6 around 200 MPa of normal stress, which generally occurs between 5 and 15 km. On the other hand, the experimental data collected by Byerlee [1978] also indicate that in case of thick fault gouges the friction coefficient could be quite lower. This feature has been documented in different geological settings where phyllosilicate minerals occur, leading to a lubricating effect [e.g. Abers, 2009;Lockner et al., 2011;Ikari et al., 2011;Middleton and Copley, 2014]. For example, evidences for low friction coefficient have been found for both lithospheric [Lachenbruch and Sass, 1980;Zoback et al., 1987;Carpenter et al., 2011;Lockner et al., 2011 among others] and smaller-scale faults [Collettini and Holdsworth, 2004;Abers, 2009;Haines et al., 2014], as well as in subduction settings [Remitti et al., 2015].
In crustal rheological modelling a value of 0.70-0.75 is commonly considered [e.g. Ranalli and Murphy, 1987; Ranalli.1995; among many others]. However, since Byerlee observations were based on intact samples, these values are likely suited for cratonic and (almost) tectonically stable regions, while in strongly deformed settings, pervasively characterized by mechanical discontinuities representing weak rock volumes, a lower value of the friction coefficient Δ = · · · ( ) ¹ is expected. In general, the larger are i) the fault surface, ii) the cumulative slip and iii) the number of seismic events that reactivated a fault, the smoother will be the rupture surface and the sliding path, and the richer in fine-grained elements is the cataclastic zone. Accordingly, an older and/or more active fault could be considered mechanically weaker than the surrounding intact rocks. As a consequence, crustal volumes containing major active faults could be reasonably characterized by a lower friction coefficient. In summary, the tectonic parameter is proportional to μ and a decrease of the friction coefficient decreases its value and vice versa.

Pore fluid pressure
In dry conditions the pore fluid pressure is nil, while in hydrostatic conditions, considering common crustal rock densities, it usually varies in the range 0.35-0.45. However, in peculiar geodynamic and tectono-stratigraphic conditions within the uppermost crust, for example within accretionary wedges, suprahydrostatic conditions ( e > 0.5) affecting large sedimentary volumes are much more likely to occur. In such a geological setting poorly consolidated ocean-bottom sediments enter the subduction zone while still being highly permeated by circulating fluids. As the subduction process advances, these sediments undergo progressive compaction due to the rapid increase of the gravity stress component [Caputo, 2005b]. This leads to considerable pore reduction and fluid pressure increase due to a poor water seepage. The Skempton coefficient could thus locally increase up to unity. In deeper crustal rocks, such extreme hydraulic conditions are rare, and the coefficient could be reasonably assumed to be roughly hydrostatic. Nevertheless, slightly suprahydrostatic conditions could be generated in rock volumes rapidly heated up due to the induced thermal expansion of the fluids and the consequent pressure increase, or in case of magma-related fluid injection. From a mechanical point of view, an increase of the pore fluid pressure leads to a decrease of the effective stress, which in general promotes failure. This diminishes the peak shear strength and thus it slightly deepens the BDT.

Creep behaviour parameters
The power-law creep behaviour (equation (5)) implies a strong dependence of the shear strength on temperature and its gradient. Accordingly, we have also analysed, in the dedicated section 4.3, the parameters that are commonly used to estimate the distribution at depth of the temperature. When dealing with the ductile behaviour, it should be also emphasized the paramount importance of the selected lithology. Unlike the brittle one, where just density is lithology-dependent, most parameters used in equation (5) can markedly vary from one lithotype to another.

Power-law parameters
The three parameters A, n, and E are commonly referred to as "power-law parameter", "power-law exponent" and "activation energy", respectively. The values proposed for several lithologies have been generally obtained by fitting the experimental data [among others, Shelton and Tullis, 1981;Hansen and Carter, 1982;Kirby, 1985] to the power-law equation, based on a least-squares regression approach. The experiments, mostly conducted in the 1980s, were carried out at temperatures between ~400 and ~1000 °C, with confining pressures up to ~1 GPa and strain rates between 10 -4 and 10 -8 s -1 . Such setting conditions require paying specific attention when extrapolating the observed data to real geological conditions. In any case, considering a dry quartzite rock, variations of the power law parameter A show a negligible effect on the BDT depth and shear strength. The power-law exponent n exerts a higher weight (as expected for simple mathematical reasons): a 10% increase (decrease) of n causes a shallowing (deepening) of the BDT in the order of 2-3 km and a decrease (increase) in peak strength of some tens of MPa.
From a physical point of view, the activation energy E represents an estimate of the amount of energy necessary to initiate the dislocation glide creep phenomena at the crystalline scale. An increase in E leads to an increase in strength, as intuitively much more energy (viz. stress) must be stored and applied on that volume to allow deformation begin.
Its variability plays a major role on the BDT depth determination and on the resulting strength. For example, a +10% variation could cause a BDT deepening of some kilometres and a strength increase in the order of 10 MPa.

Strain rate
The rate at which deformation takes place is one of the most relevant factors for the rheological modelling of the ductile behaviour (see equation (2)). The strain rate represents the long-term ductile deformation velocity and values in the order of 10 -15 s -1 are commonly considered in the literature (see Ranalli, 1995;Sibson, 1982 among others). A larger strain rate causes general strengthening associated with a BDT deepening and a peak strength increase. This is a direct effect of inhibiting the ductile deformation (rocks are more viscous) thus fostering brittle failure propagation at greater depths. As concerns the related uncertainties, we considered values in the order of ±10%, which do not induce substantial changes in both BDT depth (say, few hundred meters) and shear strength (few MPa).
Lateral variations of ± one order of magnitude, depending on the resolution of the interpolation, are normally documented by geodetic investigations carried out at the regional scale and represented by smoothly variable maps [Hollenstein et al., 2008;Floyd et al., 2010;Devoti et al., 2011;Kreemer et al., 2014;Mastrolembo and Caporali, 2017]. Estimates of strain rates from geodetic measurements, if deprived of any non-tectonic disturbance or coseismic slip effect, are considered to be representative of the long-term, interseismic deformation rates and rheological behaviours at depth in the crust and in the upper mantle [see Burgmann and Dresen, 2008, for a review].
More recently Doglioni et al. (2015) have also used space geodesy to estimate strain rate values below the BDT.
Other variations of the strain rate also certainly occur during the seismic cycle. For example, even if during the interseismic period the bulk of the deformation could be characterized by a low value, say in the order of 10 -15 s -1 , during the seismic and immediate post-seismic stage a strain rate increase of even 10 orders of magnitude can be temporarily experienced, at least by the uppermost portion of the plastosphere in correspondence of the deepest tip of the reactivated fault zone [Scholz, 1990]. Due to the transient nature of this phenomenon, attention should be paid when investigating short-term seismogenic processes. However, since we are primarily interested in the long-term rheology and its relationship with the thickness of the seismogenic layer during the longer interseismic period, the possible time-dependent short-term variations of the strain rate have been disregarded in the present paper. The reason for this choice lies in the fact that, when the seismic cycle is close to its end and a seismic sequence is about to occur, the depth of the BDT effectively corresponds to the one stabilized during the interseismic period.
In this view, it is sensible to consider the interseismic BDT depth as a constraint for the thickness of the seismogenic layer and to compare the rheological transition with the depth distribution of the seismic events, especially at the start of the sequence and in the relation with the depth of the mainshock.

Thermal gradient
Variation of temperature with depth is a crucial issue in modelling rheological profiles, as the thermal gradient represents the primary factor controlling the ductile strength and consequently the BDT depth [e.g. Sibson, 1984;Doser and Kanamori, 1986]. Unlike oceanic lithosphere, whose thermal structure is quite simple and well described by several numerical models [e.g. McKenzie et al., 2005], for the continental lithosphere we need to consider different thermo-rheological layers. For the purpose of the present paper, we distinguish sedimentary cover, upper crust, lower crust and upper mantle. Additional factors, like the occurrence of radioactive elements and their spatial distribution (especially in the crust), as well as the presence of hot igneous intrusions or peculiar structural/tectonic settings [Lucazeau and Le Douaran, 1985;Mareschal and Jaupart, 2013] could also play a significant role. All these possible phenomena and features should be properly considered for a correct and rigorous thermal modelling. At this regard, we use the 1D equation from Çermak [1982] as the one describing the temperature distribution over depth, as first proposed by Lachenbruch [1968]: (3) where T(z) is the temperature (°K) at depth z, T 0 is the temperature at the top of the considered layer, q 0 is the surface heat flow density (W/m 2 ), A 0 is the radioactive surface heat production (W/m 3 ), D is an exponential decay constant, which is also called characteristic depth and has the dimensions of a length (m), and finally k is the thermal

Massimiliano Maggini and Riccardo Caputo
conductivity (W/(m·K)). Given the importance of the temperature gradient on the strength envelopes, each parameter of the thermal equation is analysed individually in the following section due to the potential impact on the modelling.

Surface heat flow
Surface heat flow q 0 is probably the most important factor due to its effective weight on temperature (equation (3)) and hence on the ductile strength (equation (2)). It is the most readily measurable parameter on the surface that acts as a proxy for deep thermal conditions [Ranalli, 1995]. As a common rule, the surface heat flow can be viewed as an indicator of the tectonothermal and geological history, together with the geodynamic setting of a given area.
Caution should be paid when treating heat flow data from literature, as they should be thoroughly analysed and interpreted within their specific geological frame and scale of observation. Indeed, local measurements do not always represent a proxy of the actual thermal gradient at depth, as surface recordings may be affected by very localized and shallow heat sources, such as hydrothermal vents, fracture-related conduits for deeper hightemperature fluids or even the occurrence of shallow magma chambers, sills and laccoliths (see Beardsmore and Cull, 2001 for a thorough review on the possible errors and bias related to heat flow measurements and estimates). One should also consider that whenever interpolated maps are used to derive the surface heat flow, the scale and the informatic tool used for the map interpolation, together with the degree of peaks filtering and/or smoothing, play relevant roles on the resulting heat flow values. Not to mention the possible bias related to the fact that heat flow measurements are more commonly recorded by technicians and scientists where manifest surface evidence of thermal anomalies are observed, which tend to coincide with very localized heat sources, especially in case of hydrothermal and volcanic settings. Notwithstanding all these possible complications, heat flow maps still remain one of the most readily available data providing a crucial source of information for thermal (and hence rheological) modelling. In general, high heat flow values increase the thermal gradient, thus decreasing the ductile strength and leading to a remarkable shallowing of the BDT depth.
Following the above arguments, a ±10% of uncertainty for the heat flow has been assumed in the modelling, which could induce temperature variations of ±20-35 and ±80-85 °C at 15 km or at a typical Moho depth, respectively. Such differences could therefore lead to a shallowing/deepening of the BDT of some kilometres and a peak strength decrease/increase of few tens of MPa.

Thermal Conductivity
The thermal conductivity k is the coefficient linking the thermal gradient (ΔT/dz) with the surface heat flow q 0 , where the relationship is expressed by the Fourier's law for the thermal conduction. Since we are mainly interested in the crustal lithospheric modelling, we only consider conduction as the unique truly effective process of heat transfer, as far as convection and diffusion start being efficient only in the asthenosphere [see for example Ranalli, 1995;Turcotte and Schubert, 2014]. From a physical point of view, thermal conductivity represents an estimate of how "easily" and efficiently heat is transferred to the surface, given a certain thermal gradient. Although a minor dependence on temperature and density of the surrounding rocks has been documented [e.g. Birch and Clark, 1940], for thermo-rheological modelling it could be disregarded [see Beardsmore and Cull, 2001 among others]. High k value mean that heat is efficiently transferred to the surface and therefore a relatively low thermal gradient is required to fit the surface heat flow. This would result in generally cooler and hence stronger rocks and a consequently deeper BDT. After having realized several tests with different thermal conductivities and geothermal gradients, we observed that a variation of ±10% of k could induce a temperature difference of ca. ±10-30 and ±25-40 °C at 15 km or in correspondence of the Moho, respectively.

Internal radiogenic heat production
The radiogenic heat production is one of the primary sources of heat in rocks, especially for the continental, felsic upper crust. The process is due to the decay of unstable isotopes of some highly incompatible elements (Th, U, K among others) particularly abundant in the differentiated, mostly felsic, continental crust. Such a genetic process for the distribution of the radioactive elements also explains why the mantle, and to a lesser extent the lower crust, being depleted in these unstable isotopes, are characterized by noticeably lower levels of radiogenic heat production. The radiogenic heat production depth distribution was first described by Lachenbruch [1968], whose equation states that where D is the so-called characteristic depth describing how rapidly the radiogenic production decreases with depth.
In most cases, the commonly assumed value is ~8.5 km ±1.5 km [e.g. Pollack and Chapman, 1977]. A relatively high value of A 0 means that a considerable portion of the total heat transferred to the Earth surface is produced locally (and at shallow depths) by radioactive decay processes and therefore even a relatively high value of q 0 does not necessarily imply a high geothermal gradient. In other words, at a fixed depth, for example the Conrad discontinuity, in case of a high A 0 value, a lower temperature is expected with respect to the case of a low A 0 value. Accordingly, the introduction of higher A 0 values in the rheological modelling causes a BDT deepening and a ductile strength increase, though induced variations are minimal, if uncertainties of ±10% are considered.

Case studies from the Aegean Region
As already mentioned, two case studies of rheological profiles are here proposed and used as starting points for the quantitative sensitivity analyses carried out in this research. They are both from the Aegean Region and represent two Andersonian regimes (tensile and transcurrent, respectively) that largely pervade this area of the Mediterranean realm ( Figure 1). Indeed, almost the totality of the seismogenic faults in the internal Aegean domain (with respect to the Hellenic subduction arc) is characterized either by pure normal and strike-slip kinematics or by a combination of the two [see for reference the GreDaSS database from Caputo and Pavlides, 2013]. Active thrusting is limited to the central and western sectors of the Hellenic subduction zone and to the offshore regions external with respect to it. Consequently, in these regions, information on some variables such as the heat flow or the lithological stratification is even more biased with possible uncertainties and errors due to the increasing difficulty in retrieving reliable data offshore. Moreover, a precise knowledge of all the offshore active thrusts is lacking, also because of the remarkable paucity of precise seismological data regarding seismic sequences occurred in historical times. In addition, there is an intrinsic larger error in the localization of seismic events occurring offshore, due to the only partial azimuthal coverage resulting from uneven distribution of the seismic stations.
All these caveats strongly limit the possibility to obtain a reliable rheological model for the thrust-dominated region that would also need a dedicated rheological characterization and model for the subduction process, which is beyond the scope of this paper. Accordingly, we did not select any test site for the thrusting regime. In the following, the parametric selection of the two sites is firstly discussed to provide our preferred models, which are assumed as the 'reference' ones in the subsequent sensitivity tests.

Kallidromo area
The first selected case study is from the Kallidromo area, Central Greece, along the southern border of the Sperchios Basin ( Figure 2a) located in the eastern sector of the External Hellenides fold-and-thrust belt. The tectonostratigraphic and geological evolution of the area has undergone some major deformation phases, which could be synthetically described as follows: a) the mostly Cenozoic nappe stacking caused by the Alpide orogeny [Brunn, 1956;Aubouin, 1959;Doutsos et al., 1993;2006;van Hinsbergen et al., 2005;Mountrakis, 2006]; b) the mainly Pliocene NE-SW trending post-orogenic collapse [Mercier, 1976;Mercier et al., 1987;Caputo and Pavlides, 1993;Doutsos et al., 1994] or gravitational spreading [Le Pichon and Angelier, 1979], causing the thinning of the previously thickened crust and c) since Middle Pleistocene, the general rearrangement of the stress field within the Aegean domain [Mercier et al., 1979;1989;Angelier et al., 1982; among others], still pervading the whole region.
The latter tectonic regime is characterized by a ca. N-S crustal stretching and is generating new, mainly E-W transitional/oceanic sequence has a thickness of ~1 km [Smith and Moores, 1974], a reasonable value for the whole sedimentary layer thickness is ~2.5 km. The total crustal thickness can be retrieved from the regional estimates of the Moho depth. The latter in the Kallidromo massif is ~31-32 and 33-34 km, respectively according to gravity data [Tiberi et al., 2001] and tomographic inversion [Zelt et al., 2005], while Makris et al. [2013], using gravity data constrained by seismic profiles, propose a value of 32-33 km. We thus selected 32.5 km as the most likely value to represent the Moho depth at this test site. No data are available to discriminate between upper and lower crust for the investigated area. However, on a global scale, a thickness ratio of lower-to-crystalline total crust of ca. 0.4-0.45 has been proposed [e.g., Rudnick and Fountain, 1995;Rudnick and Gao, 2003]. Applying this approach to the Kallidromo area and considering an average topography of 0.5 km (as obtained from SRTM30), the assumed thicknesses for the upper and lower crust are 17.5 and 13 km, respectively (Table 1).

Massimiliano Maggini and Riccardo Caputo
Based on the geodynamics of the broader area and the local geological conditions and considering also the availability of laboratory tested rocks, the most likely lithologies for representing the three upper layers (viz. sedimentary cover, upper and lower crust) in the rheological modelling are metasediments, dry quartzite and wet diorite, respectively (see Table 1). In order to verify whether these choices are reasonable, we used the compilations by Christensen and Mooney [1995] relating seismic velocities at depth with lithologies. Two deep seismic soundings performed slightly north of the Kallidromo test site, in the Evoikos Gulf [Makris et al., 2001], have been considered, showing seismic velocities in the range 5.9-6.3 and 6.5-6.8 km/s for the upper and lower crust, respectively.
According to Christensen and Mooney [1995], both intervals are compatible with the selected quartzite and diorite lithologies showing typical velocity of ~5.9-6.0 and ~6.5 km/s at ca. 10 and 20 km-depth, respectively. On the other hand, dry quartzite and wet diorite have been already proposed for this region [Tesauro, 2009]  strength envelopes have been taken as reference models during the sensitivity analyses ( Figure 2).
All values are reported in Table 1.
A value of 0.70-0.75 for the friction coefficient μ is commonly assumed in rheological modelling (see Section 3.2).
However, due to the cumulative slip (at least since Middle Pleistocene) in the order of some hundreds of meters, a slightly lower friction coefficient (0.60) has been assumed as the preferred value for the purpose of the present paper. As far as in the Kallidromo region the tectonic regime is purely tensile, is unity and the parameter could be readily calculated following Ranalli [1995], obtaining a value of 0.68.
As concerns the pore fluid pressure, and as far as no particular reasons exist for supra-or sub-hydrostatic conditions to occur within the crust of the Kallidromo area, hydrostatic conditions could be reasonably assumed, with the Skempton coefficient e (equation (1)) being conventionally posed equal to 0.4.
As concerns the surface heat flow, we analysed the works of Cloetingh et al. [2010], Fytikas and Kolios [1979], Hurter and Haenel [2002] and Taktikos [2001]. The latter two have been however discarded. Indeed, Hurter and Haenel [2002] show a systematic anomalously high heat flow with respect to the other studies, likely due to some bias in the interpolation. As for the map of Taktikos [2001], it is characterized by the presence of some single spikelike markedly high values probably suffering i) from the scattered geographical distribution of the source data, ii) a possible bias due to a concentration of measurements at (hot) spring sites and iii) a likely unsuitable interpolation method. We thus averaged the values proposed by Fytikas and Kolios [1979] and Cloetingh et al. [2010], assigning a greater weight to the former due to their focus on the Aegean Region. In summary, the inferred surface heat flow value for the Kallidromo test site is 75 mW/m 2 .
Based on the above constraints, estimates and assumptions at the base of the preferred parametric model (Table 1), the local rheological profile was reconstructed, using a purposely developed Matlab® code (see Section 6). The results for the Kallidromo test site show the occurrence of the BDT at a depth of 12.3 km, where the temperature is 303° C and the corresponding shear strength is 136 MPa (Figure 2b). The strength envelope shows the predominance of the frictional sliding behaviour for the sedimentary layer and the shallower-medium depth portion of the upper crust, where the brittle-ductile transition occurs; in contrast, ductile creep characterizes all the remaining layers.
The shape of the envelope implies that the maximum strength along the entire profile occurs at the BDT depth.

Cephalonia area
The second case study corresponds to the area of the Cephalonia Island (Figure 3a), which is the major of the Ionian archipelago. From a tectonic point of view, this region is dominated by the presence of the Kephallinia Transform Fault (KTF), which exhibits prevailing dextral strike-slip kinematics with a slight reverse component [e.g., Scordilis et al., 1985;Louvari et al., 1999]. The KTF separates the External Hellenides accretionary wedge, still evolving along the Hellenic Arc, from the Adriatic/Apulian block. This crustal-scale shear zone is characterized by ~20 mm/a of cumulative slip rate [Serpelloni et al., 2005;Ganas et al., 2016]. indicate the occurrence of a sedimentary cover (mainly consisting of marls, sandstones, limestones, cherts, dolomites and anhydrites) up to 3.5-4 km-thick, though Makris et al. [2013] suggest a larger overall thickness of the sedimentary layer, say ca. 6.5 km, considering the nearby presence of the Mediterranean ridge accretionary prism.
The latter value has been accordingly assumed as more representative for this site. The estimates of the Moho depth are provided by several authors [Sodoudi et al., 2006;Raykova and Nikolova, 2007;Makris et al., 2013;Grigoriadis et al., 2016] using different methods, like P-and S-waves receiver functions, surface waves tomography, gravity data or gravity constrained by seismic profiles. Likely due to the different methodological approaches and original datasets, they provide a quite wide range of values between ~27 to ~40 km.
For the purpose of our rheological modelling and considering the proximity to the down going, transitional/oceanic, Ionian/Nubian lithosphere, an average Moho depth of 32.5 km has been thus assumed.
We then assigned the thicknesses of the upper and lower crust following the same approach described for the and 7·10 -15 s -1 . Also, in this case and for the same reasons, the upper bound value of 7·10 -15 s -1 has been assumed.
As concerns the coefficient of friction μ, we consider that the KTF system contains major crustal faults characterized by high slip-rates and likely quite "smoothed", suggesting a reasonably lower value (μ = 0.5).
The Cephalonia test site is affected by a transpressional tectonic regime for which a more appropriate stress ratio of ~0.4 is estimated (thus leading to = 0.98).
As regards the pore fluid pressure, considering that the investigated crustal volume belongs to the accretionary wedge being obliquely affected by the KTF, a slightly suprahydrostatic value of e = 0.7 has been considered.
The Cephalonia area is located far away from the Hellenic volcanic arc, it is characterized by a recently thickened crust and, as expected, the surface heat flow is quite low. Indeed, the local heat flow value is ~42 mW/m 2 , as obtained from the same approach of the other test site.
Based on the above constraints, estimates and assumptions providing the preferred parametric model (Table 1), the local rheological profile for the Cephalonia test area yields a BDT depth of 16.4 km, where the temperature is 246° C and the corresponding shear strength is ~125 MPa (Figure 3b). Differently from the Kallidromo test site, here two additional brittle layers occur at depths greater than the first BDT (respectively in the uppermost portion of the lower crust and in the upper mantle). The brittle layer in the lower crust has a thickness of only ~2 km while that in the mantle is much thicker (~25 km). Consequently, the maximum value of the strength (~500 MPa) along the profile is attained at the lowermost brittle-ductile transition within the mantle.

Methodological approach
As above mentioned, the principal aim of this work is to assess and analyse the sensitivity of rheological profiles to the several input parameters (equations (1) to (3)). Accordingly, each parameter has been varied within the uncertainties provided by the authors that carried out the laboratory experiments, or considering the constraints provided by a careful geological, tectonic and evolutionary analysis of the selected lithospheric volumes (see also Section 5). If the uncertainty cannot be defined otherwise, the variability is assumed to be ±10% with respect to the preferred value. We also tested parametric variations of ±5% and ±15%, but we did not find any significant impact on the sensitivity results relative to the ±10% variation. For both test sites, the rheological profiles were firstly produced based on the preferred values and subsequently recalculated using the extreme uncertainty values separately for each parameter, while keeping fixed all the others. All the uncertainties and the results in terms of BDT depth, corresponding temperature and strength are reported in Tables 2 and 3, where the obtained values are compared to the reference model, in order to emphasize the influence of the different parameters on the rheological modelling.
As a reasonable compromise between the need for precise, high-resolution profiles and computations running time, the code for implementing the rheological modelling uses 100 m steps along the vertical. The strength envelopes based on the reference model for the two test sites are respectively represented in Figures 2b and 3b, where the strength is on the right side of the graphs while temperature is on the left side.
Along the vertical profile, the BDT is pinpointed where ductile strength becomes smaller than frictional one and the prevailing deformation mechanism switches from stick-slip, potentially seismogenic behaviour to slow, continuous and distributed creep.
In case of occurrence of a ductile layer embedded within two brittle layers, the BDT has been selected in correspondence of the shallower brittle-to-ductile transition, provided a sufficient thickness of the ductile layer is interposed. At this regard, seismological evidences [e.g. Rolandone et al., 2004] and mechanical-rheological modelling [Beeler et al., 2018] have proved that the embrittlement caused by the temporary increase of the strain rate after a major seismic event may result in a sudden and transient deepening of the BDT, with respect to the interseismic period depth, in the range of 1-3 km. Indeed, aftershocks are sometimes recorded at depths slightly greater than those of the mainshock and of the background seismicity during 'quiescence' periods. Considering that the seismotectonic setting discussed by Rolandone et al. [2004] and Beeler et al. [2018] strongly differ from our selected test sites, as far as i) our faults are inclined and not vertical, ii) the maximum expected magnitude is

Single parameter effects
As concerns the effects of the single-parameter variations on the BDT depth of the Kallidromo area, the most influencing ones belong to the constitutive equation of the power-law creep deformation mechanism. In particular, the activation energy causes the strongest BDT depth variation (±2.3 km), while the uncertainty on the exponent n also induces a +0.8/-0.9 km change. Such effects are likely due to the high reference value (10 5 J/mol) and the exponential role in equation (2) geothermal gradient, induce BDT changes of about 1 km and up to 1.4 km. As discussed in a previous section, this is in line with the weight exerted by the thermal properties on the depth distribution and layering of brittle versus ductile behaviours. As concerns the other tested parameters, and especially those used for describing the brittle behaviour (equation (1)), none of them seems to have a relevant influence on the BDT depth, inducing variations in the calculated values of only few hundred meters (maximum 0.4 km).
As regards the strength at the BDT depth, a major role is conversely played by the choice of the parameters describing the frictional sliding behaviour, and particularly the friction coefficient that may induce a strength variation of 25 MPa if a low value (μ = 0.4) is assumed. Similar changes (-20/+18 MPa) are caused by ±0.1 variation of the pore fluid factor, while the tested ranges of , έ and have a limited effect by inducing less than 10 MPa difference in critical differential stress. The relevance of friction and fluid pressure on the strength is actually an expected feature, since these factors directly affect the capability of rocks to build up and sustain a certain amount of stress. However, also the parameters governing the 'thermal-ductile' behaviour of rocks (E, q 0 , n and k) could induce comparable variations of several MPa from ±9 and up to +24/-25 MPa (e.g. activation energy).
A further property we focused on for assessing the relative importance of all constitutive parameters for the thermorheological modelling is the temperature at the BDT depth, which can also give insights on the mineralogical and mechanical characteristics of the typical rocks at the corresponding depths. Once more, the most influencing parameter is the activation energy, whose ±10% 'uncertainty' relative to the reference value leads to temperature changes of almost 50° C. However, since the activation energy and to a less extent the power-law exponent, though affecting the ductile creep strength, are not directly controlling the geothermal gradient, their influence on the temperature is simply due to the vertical shift of the BDT itself. Perhaps quite unexpectedly, the tested uncertainty of the other "thermal" parameters (namely, the heat flow, the thermal conductivity and the radiogenic heat production) only causes very slight changes of the BDT temperature, in the order of a few degrees. For these parameters, the apparent lack of influence is likely due to the fact that, even though the geothermal gradient is substantially modified, at the same time also the BDT depth is shifted, and the combined effect ultimately acts to produce an almost zero net variation of the BDT temperature. None of the other parameters produce temperature variations at the BDT within the tested ranges.
The Cephalonia test site is characterized by a deeper BDT with respect to Kallidromo (16.4 vs 12.3 km), though exhibiting a comparable strength at the BDT (125 vs 136 MPa). This is likely due to the different lithologies and, above all, fluid pressure conditions at depth.
Analogous to the Kallidromo case study, the most relevant BDT depth variations are obtained when testing the uncertainties on the activation energy, inducing +4.3/-4.0 km changes. Also, the heat flow induces a sensible BDT depth variation (-1.8/+2.4 km). On the other hand, n and k, among the other 'creep' parameters, and μ and e among the parameters describing the brittle behaviour have a limited influence (in the range -1.3 and +1.4 km, but with opposite roles). The effects of the other parameters are instead negligible.
As concerns the critical differential stress, the greatest weight is exerted by the Skempton coefficient (-36/+66 MPa) and by the friction coefficient (+38/-41 MPa). In particular, pore fluid pressure could locally attain markedly variable values within an accretionary wedge setting [e.g. Suppe, 2014] such as that of Cephalonia and consequently it plays a major role in defining the stress tensor distribution at depth and the total strength of rocks. On the other hand, also the friction coefficient could be characterized by possible vertical (and lateral) variations due to the presence of asperities and irregularities along the fault surface of a major and well developed structure such as the Kephallinia Transform Fault; this effect should be carefully taken into account in computing the BDT strength when modelling the rheology of crustal volumes containing major shear zones. Similar to the Kallidromo case study and likely for the same reasons, also in terms of critical stress, the tested range of the activation energy causes important differences in the calculated strength (+32/-31 MPa). Although to a less extent, also the uncertainty on the heat flow and secondarily on the power-law exponent and factor could induce less than ±20 and ca.

Combined effects
We also analysed the combined effects on the BDT depth and the associated values of strength and temperature, by simultaneously varying all tested parameters. In order to achieve this goal, a statistical approach has been followed.
Firstly, we randomly generated a set of 10,000 values for each tested parameter, in which the mean and the 95% confidence interval are based respectively on the preferred value of the reference model (Table 1) and the extremes of the uncertainty range reported in Tables 2 and 3. For the distribution of the values, a Gaussian normal distribution has been generally assumed with the only exceptions of the friction coefficient (for both test sites), the Skempton coefficient e for Cephalonia and the tectonic regime factor for Kallidromo. Indeed, being characterized by asymmetric ranges with respect to the reference model values, they are better represented by skewed normal distributions, or in the case of for Kallidromo, by a truncated normal distribution.
Secondly, assuming the variability of each parameter as basically independent from the others, we randomly combined the above sets (of 10,000 values each) for calculating as many rheological profiles for each case study. The results for Kallidromo and Cephalonia test sites are represented in Figures 4 and 5, respectively, as probability density functions graphs of the three main descriptive indicators, namely, the BDT depth, strength and temperature with bins of 0.5 km, 5 MPa and 5° C, respectively. In order to verify the robustness of the statistical results, for each analysed indicator, the procedure has been applied to ten 10k populations and as many curves are plotted in Figures   4 and 5 (thin lines). The light grey thick line represents the average of the ten datasets (viz. 100k random combinations), which confirm the stability of the outcomes.
It is noteworthy that the Kallidromo distribution for the BDT depth ( Figure 4a) follows a quasi-normal geometry, slightly positively skewed. For each of the ten plots the mode, expressed in terms of the most populated depth class having a dimension of 0.5 km, has been selected and the ten values averaged, showing that such a mean value basically coincides with the BDT depth obtained from the reference model (bin 12.0-12.5 km versus 12.3 km). The standard deviation for the BDT depth over the whole 100k random parametric combinations is quite limited, being ca. ±1.5 km, thus indicating that all the values are focused around the mean and occur in the upper crust layer.
The frequency distribution of the strength at the BDT depth ( Figure 4b) is slightly more positively skewed than the BDT depth (quasi-normal) distribution. For this parameter the most frequent class corresponds to 130-135 MPa, in good agreement with the reference model value (136 MPa). The slightly longer tail on the greater values may be due to the effect of the factor, whose assumed distribution is obviously truncated at 1 (the reference model value) and results in slightly greater strengths for all the other values of the distribution. The standard deviation is up to ~22 MPa, mostly because of the uncertainties in the pore fluid pressure, the friction coefficient and the abovementioned tectonic regime factor . Also, as regards the BDT temperature distribution (Figure 4c), this resembles a normal one, with the peak in the bin 295-300° C, being compatible with the expected values for the onset of plasticity in quartz, typically placed at 300 ±50° C (e.g. Scholz, 1989). The mean value of the 100k-combinations is ~303° C (standard deviation of almost 27° C), perfectly coinciding with that obtained from the reference model (303° C).
As concerns the Cephalonia test site, the BDT depth frequency distribution basically follows a quasi-normal shape (Figure 5a), though the deeper tail in the upper crust is slightly truncated due to the lithological stratification (upper/lower crust boundary), while a minor cluster occurs around bin 26.0-26.5 km-depth (i.e. in the lower crustal layer). However, as far as this cluster represents less than 4% of the total random combinations, it could be considered statistically meaningless and therefore the distribution should be not effectively considered as bimodal.
The calculated average value of the 10 10k-distributions most populated classes falls in the bin 16.0-16.5 km and the 100k-combinations mean is 16.5 km (while the standard deviation is ±3.1 km), in perfect agreement with the reference model value of 16.4 km.
The distribution of the peak strength is instead clearly unimodal and associated to a roughly normal distribution, even though with a sensible positive skewness (Figure 5b). This difference could be possibly, but not exclusively, related to the few BDTs occurring in the lower crust, since the strength is not solely dependent on depth, as it is instead the case for the BDT temperature. The average of the most populated classes for the ten 10k-datasets falls in the bin 115-120 MPa, with a 100k-combinations mean of 136 MPa (standard deviation of ca. ±43 MPa), in reasonable agreement with the value obtained on the basis of the reference model (125 MPa).
As previously mentioned, the temperature is directly dependent on depth; therefore, as expected, the overall distribution of the BDT temperatures obtained for the Cephalonia test site closely resembles that of the BDT depth, showing a statistically negligible peak around bin 340-345° C (Figure 5c), also in this case representing less than 4% of the results. The mean value of the 100k-combinations and the associated standard deviation are respectively 247 and ±29° C, while the most populated class is the range 245-250° C. All statistical indexes are in good agreement with the temperature value obtained on the basis of the reference model (246° C). The lower BDT temperature with respect to the Kallidromo test site (ca. 300° C), even though the BDT is deeper in the Cephalonia area, is likely a consequence of the very low geothermal gradient characterizing the Ionian region with respect to Central Greece.

Massimiliano Maggini and Riccardo Caputo
18 is undergoing crustal stretching, thinning and warming at least since the Pliocene.
We also tested the respective dependence of BDT depth, strength and temperature for both test sites considering all the ten 10,000 combinations, by analysing their distributions as a function of each other. In Figure   6 and Figure    Cephalonia one. We argue that this is due to the higher geothermal gradient and the fact that all BDTs occur in the upper crustal layer for the Kallidromo test site, thus limiting the spread for the values of the BDT strength and temperature. It can be also noted that the depth versus temperature plots in both cases show the greatest correlation. This is intuitively attributed to the direct proportionality between temperature and depth (temperature as a general rule always increases with depth), which instead does not characterize the depth versus strength and strength versus temperature relations.

Comparison with seismicity
In order to further verify the reliability of our rheological modelling and sensitivity analyses for both test sites, the BDT depths have been compared with the local relocated seismicity data available from the literature. For this purpose, we considered the distribution at depth of relocated events from two recent seismic sequences, occurred in August 2013 and during January 2014 in the Kallidromo area and Cephalonia Island, respectively.
As already mentioned, the brittle-ductile transition could be considered as a reasonable approximation of the maximum depth at which seismogenesis occurs. More precisely, no major earthquake is expected to nucleate below the BDT, even though coseismic rupture during strong events may propagate below the transition [e.g., Scholz, 1988]. Indeed, in exhumed crustal shear zones the presence of pseudotachylites have been documented embedded within a mylonitic fabric matrix; the former clearly documenting large amounts of coseismic frictional sliding (and heat production high enough to induce local melting), while mylonitic rocks are typical of distributed (viz. ductile) deformation [e.g. Sibson, 1980;Passchier, 1982;Fabbri et al., 2000]. Bearing in mind the possibility that the rupture propagation on major crustal faults during strong events could extend somehow deeper, the BDT depth as inferred from rheological modelling could nevertheless represent a crucial information to constrain the width of the seismogenic faults, especially in case their geometry at depth is still debated.
It is noteworthy that the aftershocks spatial distribution commonly highlights and focuses on the boundaries of the main patch ruptured during the mainshock and if the rupture reached the bottom of the seismogenic layer, some of the subsequent foci will be located below the long-term (viz. interseismic) BDT.
Based on local velocity structure models and relocation procedures of the hypocentres (such as HypoDD) (Waldhauser and Ellsworth, 2000), detailed seismological and tectonic analyses have been carried out for the Kallidromo and Cephalonia areas affected by two recent seismic sequences [e.g. Ganas et al., 2014;Karastathis et al., 2015]. The use of precisely relocated seismicity arises from the need to compare our results with an independent data source, namely the hypocentral distribution, with appropriate accuracy and resolution.
The 2013 Kallidromo sequence occurred on a secondary fault antithetic to the main fault of the Tithorea seismogenic source [GRCS434 in Caputo and Pavlides, 2013], which represents one of the major normal faults characterizing the Sperchios rift area. The magnitude of the mainshock was 5.4 and though the maximum depth reached by some of the aftershocks [from Ganas et al., 2014; 90 days period from the mainshock of 7 August 2013] is up to 15 km (Figure 8), it is common practice in this kind of analyses to consider the depth including the 95% of the events, which in our case corresponds to 12.6 km. If we instead apply the empirical relationships between magnitude and energy [e.g. Kanamori, 1983], the 95% of the energy release has occurred within the first 13.2 km.
Both values are in reasonable agreement with the results from the rheological modelling based either on the preferred values (12.3 km; Table 2) or on the most likely depth according to the statistical analysis (12.3 km; Figure 4).
As concerns the 2014 Cephalonia sequence, it was characterized by two mainshocks of similar magnitude (M W 6.0 and 5.9) occurred on the 26th of January and 3 rd of February, respectively [e.g. Papadopoulos et al., 2014;Karastathis et al., 2015]. The causative faults are represented by antithetic (westward steeply dipping) and subvertical splay structures [Briole et al., 2015] belonging to the crustal-lithospheric boundary represented by the Kephallinia Transform Fault which separates the External Hellenides from the Apulian Foreland. According to the distribution of the relocated hypocentres [from Karastathis et al., 2015; time period from 26 January to 15 May] the maximum locally recorded depth is about 18 km, while the 95% of the events occurs within the first 13.4 km (Figure 9). If we consider the energy distribution, the 95% has been released above 16.3 km, being in perfect agreement with the BDT depth calculated on the basis of the rheological modelling (16.4 km; Table 3 and Figure 5).   Figure 3. On the right column are shown the statistical distribution of hundred thousand random combinations discussed in the text (see also Figure 5).

Concluding remarks
In the present paper we revised the largely used approach for rheological modelling based on the frictional sliding and the power-law creep behaviour of rocks [e.g. Sibson, 1982;Ranalli, 1995]. The major goal of this investigation was to analyse the sensitivity of the modelled results relative to the most important parameters included in the constitutive equations (1) to (3). The analyses were performed by means of several numerical tests. At this regard, a dedicated Matlab code has been implemented for providing 1D rheological profiles.
We firstly analysed two case studies from the Aegean Region, belonging respectively to a tensile and a mainly transcurrent tectonic regime (Kallidromo and Cephalonia Island areas, Figures 2 and 3), both typical of this Mediterranean sector. The "preferred" rheological profiles of these two test sites provide a BDT depth, strength and temperature of 12.3 km, 136 MPa and 303° C, and 16.4 km, 125 MPa and 246° C for the two investigated areas, respectively (Figures 2b and 3b).
In order to analyse the effects of the uncertainties in the values of the relevant parameters, we carried out numerous tests by varying a single input parameter per time within a reasonable range, as discussed in the text and reported in Tables 2 and 3. The results of the tests on the single parameter effects show a general strong influence of the activation energy E and secondarily of the power-law exponent n on all three principal outputs (BDT depth, strength and temperature) for both case studies. This is likely due to the high reference value (10 5 J/mol) and the exponential position in equation (2) of the two parameters, where even small differences could play a sensible role. Also the uncertainty on the conductivity k could have some minor effect on the BDT depth and associated strength. Only future laboratory researches on all these parameters could potentially reduce these uncertainties, though some implicit variability will always remain due to unpredictable lithological variations within the averaged crustal layer properties. On the other hand, the tested ranges for the radioactive surface heat production A 0 and the so-called power-law parameter A have no significant influence on the final rheological models.
Moreover, the choice of the heat flow q 0 has some considerable influence on the BDT depth and strength, but not on the associated temperature. In this case, however, the reduction of uncertainty for this parameter could be obtained by dedicated fieldwork on the investigated areas or by means of an improved generation of regional heat flow maps.
Indeed, the available ones for the Aegean Region [Fitikas and Kolios, 1979;Taktikos, 2001;Cloetingh et al., 2010] are commonly based on relatively few data and show marked differences among them, likely due to the different datasets used by the authors and the different approaches applied for the areal interpolation.
Among the parameters describing the brittle behaviour of rocks, only the coefficient of friction μ and the Skempton coefficient e seem to have some influence on the BDT properties (in the case of the Cephalonia case study), though limited to the critical shear stress (Tables 2 and 3). Secondary and almost negligible effects are associated with the rock density and the stress ratio .
We also analysed the combined effects and the potential synergistic role of the different parameters on the rheological modelling, by applying a statistical method based on the generation of 10 random datasets, for each tested parameter, each one formed by 10,000 values. The probability density functions graphs of all the three main descriptive indicators (i.e. BDT depth, strength and temperature; Figures 4 and 5), for both case studies, resemble normal slightly skewed distributions, where the values of the most populated classes nicely agree with those obtained on the basis of the preferred model.
The final objective of all these sensitivity tests is, on one side, to gain a more precise control on the reconstruction of rheological profiles and, on the other hand, to highlight the statistical weight of each parameter, whose variability can potentially affect the rheological modelling. A better knowledge of their quantitative influence could, on its turn, result in more suitable and adequate rheological profiles for lithospheric-scale geodynamic studies and for comparison with seismic distribution at depth or seismotectonic analyses.
As a further validation of the followed approach and the obtained results, the rheological models have been compared with the relocated hypocentral distribution of two seismic sequences that recently affected the two investigated areas Karastathis et al., 2015]. It is thus observed the good fit between the depth of the rheologically modelled BDT (either based on the preferred model and the statistical distribution of the ten 10k random combinations) and the depth within which 95% of seismic energy has been released. This comparison and particularly the latter observation confirm the rheological modelling as a powerful tool for constraining one of the principal seismotectonic parameters, namely the maximum depth of the seismogenic sources. Based on the dip angle of major crustal faults, the maximum fault width could be constrained, representing a crucial information in seismic hazard assessment analyses.