Severine Furst
email: severine.furst@gm.univ-montp2.fr

Michel Peyret

Jean Chéry

Bijan Mohammadi

LITHOSPHERE RIGIDITY BY ADJOINT-BASED INVERSION OF INTERSEISMIC GPS DATA, APPLICATION TO THE WESTERN UNITED STATES

Keywords: GPS, interseimic velocity, effective rigidity, global optimization, San Andreas Fault system, uncertainties

While vertical motion induced by long-term geological loads is often used to estimate the flexural rigidity of the lithosphere, we intend to evaluate the shear rigidity of the lithosphere using horizontal motion. Our approach considers that the rigidity of the lithosphere may be defined as its resistance to horizontal tectonic lateral forces. In this case, a spatial distribution of an effective shear rigidity can be estimated from the analysis of the interseismic velocity fields. We consider the Western United States zone where weakly strained areas (e.g., the Sierra Nevada) are connected with areas of large strain rate (e.g. San Andreas Fault system).

By inverting interseismic strain distribution measured by geodetic methods, we infer the effective shear rigidity of the lithosphere. The forward problem is defined using the equations of linear elasticity. The inversion relies on the minimization of the sum of a quadratic measure of the differences between measured and modelled velocity fields. The functional also includes regularization terms for the parameters of the model. The gradient of the functional with respect to the minimization parameters is computed using an adjoint formulation. This permits the treatment of large dimensional minimization problems. Finally, a measure of the uncertainty of our inversion is illustrated through the covariance matrix of the parameters at the optimum. The optimization chart is validated on two synthetic velocity distributions.

Then, the effective shear rigidity variations of the Western United States are estimated using the CMM3 interseismic velocities. The inversion displays low effective rigidities along the San Andreas Fault system, the Mojave Desert and in the Eastern California Shear Zone, while rigid areas are found in the Sierra Nevada and in the South Basin and Range. Finally, we discuss the differences between our strain rate and rigidity maps with previously published results for the Western United States.

INTRODUCTION

Geological strain occurring over millions of years results from the continuous accumulation of anelastic processes in the crust and in the lithosphere in response to plate motion. Active deformation areas are identified by seismicity and geodetic deformation. In active deformation area, the comparison of plates motion from geology and geodesy, at these two different time scales, provides a fair agreement in term of horizontal velocities [START_REF] Sella | REVEL: A model for Recent plate velocities from space geodesy[END_REF]. Geologic and geodetic comparisons can also be made across active faults using standard models for interseismic strain [START_REF] Mccaffrey | Block kinematics of the Pacific-North America plate boundary in the southwestern United States from inversion of GPS, seismological, and geologic data[END_REF][START_REF] Meade | Block models of crustal motion in southern California constrained by GPS measurements[END_REF][START_REF] Savage | Geodetic determination of relative plate motion in central California[END_REF]). It appears that most of the documented faults display a close agreement between geodetic and geologic strain rates [START_REF] Vernant | What can we learn from 20years of interseismic GPS measurements across strike-slip faults?[END_REF].

From a mechanical viewpoint, the close agreement between short and long-term strain rates (i.e. time scales from 10 yrs to 1 Myrs) probably reflects the stability of the stress balance in the lithosphere under the action of slowly evolving remote forces associated to subduction, basal drag and, more generally, the plate system gravitational potential energy. Under the action of these forces, strain distribution is mostly controlled by the lithospheric strength. By strength we mean the maximum force sustainable by the lithosphere. Like lithospheric stress, lithospheric strength cannot be determined precisely with depth, unless with crude rheological yield strength envelope models [START_REF] Tesauro | Ductile crustal flow in 265[END_REF]. Indeed, a precise strength estimate with depth would require a detailed knowledge of the temperature profile with depth, lithology and water contents, as well as friction law in the brittle domain and temperature dependent viscous laws in the ductile crust and mantle. Therefore, the lithospheric strength can only be approached through its integral measure along depth, with numerical models of the lithosphere. These solve stress equilibrium using elasto-visco-plastic laws with prescribed boundary conditions [START_REF] Bird | Computer simulations of California tectonics confirm very low strength of major faults[END_REF][START_REF] Chéry | An integrated mechanical model of the San Andreas Fault in central and northern California[END_REF]. However, a simplified version of lithospheric strength is embedded in the concept of effective elastic thickness (EET) applied to plate flexure. Indeed, it has been shown that plates submitted to topographic and other internal loads display vertical motions controlled by plate rigidity [START_REF] Watts | Isostasy and flexure of the Lithosphere[END_REF].

Combined analysis of topographic and gravimetric signals allows for computing effective elastic thickness and its variation at continental scale [START_REF] Lowry | Flexural rigidity of the Basin and Range-Colorado Plateau-Rocky Mountain transition from coherence analysis of gravity and topography[END_REF][START_REF] Pérez-Gussinyé | Spatial variations of the effective elastic thickness, Te, using multitaper spectral estimation and wavelet methods: Examples from synthetic data and application to South America[END_REF]. A fair agreement is generally found between heat flow and EET where small values of EET correspond to high heat flow zones.

Both lithospheric strength and effective elastic thickness are commonly associated with the long-term behaviour of the lithosphere. However, these concepts can be adapted in order to interpret interseismic geodetic measurements [START_REF] Chéry | Geodetic strain across the San Andreas fault reflects elastic plate thickness variations (rather than fault slip rate)[END_REF]. For a typical time of 10 years of geodetic observation and in the absence of large earthquakes, a linear evolution of GPS motion is often observed. Therefore, a collection of GPS velocities may be used in order to compute strain rate maps at plate scale [START_REF] Kreemer | A geodetic platemotion and Global Strain Rate Model[END_REF]. Even if this latter analysis is purely kinematic, the resulting geodetic strain rate must satisfy stress equilibrium over the time of observation. In such a problem, the unknown is the incremental lithospheric strength.

One example is the spatial variation of the stress change integrated over the depth over the time of GPS observation. The problem can be simplified assuming that lateral strength variation is modulated by geodetic plate thickness [START_REF] Chéry | Geodetic strain across the San Andreas fault reflects elastic plate thickness variations (rather than fault slip rate)[END_REF]. The integrated value of the shear-stress at depth is what we call the effective shear rigidity. It is conceptually similar to the flexural rigidity: the effective shear rigidity expresses the resistance of the lithosphere to lateral forces (unit is N), while the effective flexural rigidity is related to the resistance of the lithosphere to vertical bending (unit is N.m).

In [START_REF] Chéry | Plate rigidity inversion in southern California using interseismic GPS velocity field[END_REF], we proposed a global optimization approach to estimate effective plate rigidity maps by the inversion of a GPS velocity field. The inversion provides a rigidity field realizing a RMS between the observed and modelled velocity fields close to 2 mm/yr for a dataset in southern California. However, we faced difficulties to properly fit high velocity gradients in the vicinity of the San Andreas Fault system. This is because the method did not allow the consideration of large inversion problem and therefore the local spatial density of our model parameters was too low. Moreover, a priori velocity boundary conditions were necessary and no uncertainties estimated.

In this paper, we present an enhanced version of the method to address the previous issues: the number of optimization variables can now be arbitrary thanks to the use of an adjoint formulation of the forward problem. This permits high spatial resolution for the rigidity. boundary conditions are not anymore prescribed but now treated as optimization variables as well.

uncertainties are calculated for optimal rigidity value.

The paper is organized as follows: (1) we describe the new features of the method and we state the differences with respect to [START_REF] Chéry | Plate rigidity inversion in southern California using interseismic GPS velocity field[END_REF], [START_REF] Chéry | Geodetic strain across the San Andreas fault reflects elastic plate thickness variations (rather than fault slip rate)[END_REF] we demonstrate the efficiency of our new approach on a synthetic dataset that mimics a strike slip fault locked at depth, (3) we propose a refined rigidity map of southern California and we study the sensitivity of the solution of the inversion problem with respect to the location of domain boundaries. Finally, (4) we compare and discuss our results with those already published, both in terms of strain rate maps and effective elastic thickness.

GOVERNING EQUATIONS AND FORWARD MODELLING

Geophysical laws provide the mathematical framework to compute the outcome of some physical processes: this is called the forward problem. In other words, the model and its inputs are known and specific data (e.g. seismic, geodetic, or magnetic) are sought thanks to the equations linking the physical parameters to the solution at the observation location. Most of the time, we only have access to the consequences of a physical process (e.g. the geodetic measurements). These consequences need to be inverted to determine the physical properties of the Earth interior (Tomography: e.g. [START_REF] Montelli | Global P and PP traveltime tomography: rays versus waves[END_REF][START_REF] Tanaka | P-wave tomography of the mantle beneath the South Pacific Superswell revealed by joint ocean floor and islands broadband seismic experiments[END_REF] Volcanoes and geothermal zones: e.g.; [START_REF] Anderson | Bayesian inversion of data from effusive volcanic eruptions using physics-based models: Application to Mount St. Helens 2004-2008[END_REF][START_REF] Dzurisin | A comprehensive approach to monitoring volcano deformation as a window on the eruption cycle[END_REF][START_REF] Mossop | Volume strain within The Geysers geothermal field[END_REF]. Application to reservoirs: e.g. [START_REF] Hesse | Joint inversion in coupled quasi-static poroelasticity[END_REF]. In some cases, there are analytical s the observations. For most geophysical problems, the limited amount of data used to reconstruct a model with infinite degrees of freedom leads to the non-uniqueness of the solution. Consequently, the inverse problem only provides one of the many models that explain the data and has uncertainty because the real data are subject to uncertainties and errors.

The effective elastic thickness of the lithosphere can vary laterally due to both elastic properties and the rheological failure properties that limit elastic strength. Flow strength depends on other factors than the temperature. Also, part of the variation imaged by the geodetic technique is probably due to the limits of frictional strength on faults [START_REF] Bird | Computer simulations of California tectonics confirm very low strength of major faults[END_REF]. The thermal plate regime probably exerts a large influence due to the sensitivity of the effective plate rigidity with respect to its temperature profile [START_REF] Watts | Isostasy and flexure of the Lithosphere[END_REF]. Here, we model these rigidity variations as lateral variations of the elastic properties of a plate with constant thickness. Thus, our forward model is made of a domain ( ) symbolizing a 2-D plate, which can deform according to linear and isotropic elasticity (Fig. 1). Along the boundary of the domain , we apply Dirichlet conditions (i.e., in-plane velocities ) and assume free normal traction at the surface of the plate (plane stress assumption). This hypothesis means that strain perpendicular to the plane can occur. The forward model is therefore composed of three equations, the stress equilibrium (Eq. 1), a constitutive equation linking the strain rate to the stress rate for a 2-D plate (Eq. 2) and boundary conditions (Eq.

3): on ( 1) on ( 2)

on (3)
where and are the strain-and stress-rate tensors, is the Kronecker delta function, and is assumed to be constant and equal to 0.25. distribution of (Fig. 1). Eventually, we try to fit with the observations .

INVERSION METHOD

Running the direct problem requires the prescription of the velocity on the boundary nodes and the rigidity for each mesh element. Contrary to the approach proposed in [START_REF] Chéry | Plate rigidity inversion in southern California using interseismic GPS velocity field[END_REF], the boundary conditions are not imposed anymore in the inverse problem and are treated as optimization parameters. We associate one rigidity parameter to each mesh element leading to a very large optimization problem. Our global optimization algorithm requires the gradient of the functional. We consider an adjoint formulation of the forward model to access this gradient with respect to all the model parameters simultaneously.

Cost function

We want to invert observed data and determine the model parameters minimizing the distance (here -norm) between the observed data and the predicted field inside the domain :

(4)

where is the cost function to minimize, and subscript means that the -norm is weighted by the inverse of the covariance matrix of the geodetic measurements.

Geophysical inverse problems are usually ill-posed and need to include a subjective degree of regularisation to achieve relevant geophysical solutions (e.g. [START_REF] Zaroli | An objective rationale for the choice of regularisation parameter with application to global multiple-frequency S-wave tomography[END_REF]. We therefore introduce to the cost function two Tikhonov regularization terms to control local fluctuations of the parameter vector [START_REF] Tarantola | Inverse Problem Theory and Methods for Model Parameter Estimation[END_REF][START_REF] Tikhonov | On the stability of inverse problems[END_REF]. We separate the regularization of the parameters along the boundaries and those associated to the rigidity inside the domain:

(5) where and are regularization operators. The former acts over the domain and controls the regularity of the rigidity distribution, while the latter monitors the regularity of the boundary conditions. Both are particular forms of non-linear Laplace-Beltrami operators with a local control of the level of regularization [START_REF] Mohammadi | Applied Shape Optimization for fluids, 2nd Editio[END_REF]. The weights have to be chosen by the user. Series of different optimizations have been run to highlight the effect of on the inversion. By doing so, our goal is to adjust to the data , while preserving some degree of regularity on both the rigidity inside the domain and the velocities along the boundaries. However, for each simulation, we only have the values of the velocities (and not rigidity) to compare with. Hence, adjusting the regularity of the rigidity is largely subjective and we found that using no regularization ( ) for rigidity leads to acceptable spatial rigidity gradients. Hence, in this study, we only consider the regularization term of the boundary conditions.

In order to choose , we explore the trade-off between the residual data misfit and the regularization term . This is featured in a trade-off or Pareto curve, which gathers all feasible solutions that cannot be improved in any of the objectives without degrading the other objectives (e.g. [START_REF] Vassilvitskii | Efficiently computing succinct trade-off curves[END_REF]. The selection of an optimally regularized solution depends upon the requirements of a particular study. We will illustrate the impact of the regularization over the boundary conditions for the rigidity inversion in Southern California.

Global optimization

We apply a global optimization algorithm [START_REF] Ivorra | Design of code division multiple access filters based on sampled fiber Bragg grating by using global optimization algorithms[END_REF] to iteratively invert interseismic geodetic data. Global optimization is necessary as we have no information on the convexity of the cost function and several local minima can be present. The global optimization strategy is meant to improve the initial condition for classical gradient-based methods looking for an initialization in the attraction basin of the global optimum [START_REF] Mohammadi | Applied Shape Optimization for fluids, 2nd Editio[END_REF].

In addition to the Tikhonov regularization mentioned above, the gradient of the functional is smoothed (Mohammadi andPironneau, 2004, 2009) in order to control the regularity of the parameters. The optimization algorithm ends when the functional or the variations of the gradient are smaller than some user-defined thresholds. A synthetic flow chart of the inverse problem is given in Fig. 2.

Functional derivatives computation is done using an adjoint formulation of the forward model (e.g. [START_REF] Plessix | A review of the adjoint-state method for computing the gradient of a functional with geophysical applications[END_REF]. In most of the inverse problems in geophysics, the cost function cannot be analytically linearized. If a finite difference approach is adopted, the number of forward computations for assessing the gradient of the functional is proportional to the number of parameters. Let us briefly recall the adjoint technique. The gradient of the functional with respect to the model parameters can be expressed as follows:

(6)

where is the functional, the parameters and the velocity calculated at each node of the mesh. From the equilibrium equation , we incorporate the rigidity matrix and the stress vector in (Eq. 6):

(7)

Defining the adjoint variable as the solution of the system (because is self-adjoint in our case), we obtain:

(8)
Consequently, the amount of computation needed to obtain the gradient of the functional mostly corresponds to the solution of one forward model, by opposition to a finite difference scheme which needs a number of forward model solutions equal to the number of parameters.

Model Parameters initialization

Real GPS datasets present large spatial variations of density measurements. GPS stations are usually set to observe the velocity gradient around fault zones. Therefore, we usually expect null strain in geographical areas where measurements are sparse. This information can be used to define the initial guess for the lithosphere rigidity in the optimization procedure. This is similar to what is done in topological optimization (e.g. [START_REF] Allaire | Structural optimization using sensitivity analysis and a level-set method[END_REF] where the initial structural rigidity is set to the maximum admissible value. Optimization then aims at making the structure softer and softer. A common problem in mechanical structure design is to optimize the topology of an elastic structure given certain boundary conditions. Optimality implies to minimize the weight, but at the same time, the structure needs to be as strong and rigid as possible. The rigidity of each element is hence reduced at each iteration of optimization when requested.

Synthetic and real cases presented in this paper involve rigidity reaching very large values in areas that exhibit little internal deformation. Thus, the rigidity amplitude ranges from a given minimum to infinity in no-deformation zones. This semi-open variation domain is not suitable for numerical search. Consequently, we choose a parametrization design using the compliance, , of the material instead of the rigidity. The compliance is defined over the interval , the lower bound corresponding to a quasi-rigid body. We have considered different values of . It appears that a value of which corresponds to two order of magnitude admissible variation for the rigidity is sufficient to fully capture the range of most strain-Appendix B). This use of compliance insures greater stability of the inversion process. For ease of understanding and interpretation, we express our results in terms of rigidity after the inversion is completed.

Model parameters uncertainty

GPS observations are plagued with uncertainty due to various factors: instrumental noise, field measurement procedure, the skill of the operator and local environmental motions. These uncertainties affect in a complex way the GPS time series and generate a coloured noise on positions [START_REF] Mao | Noise in GPS coordinate time series[END_REF]. But, these also induce uncertainties on the model parameters determined through our optimization process. For that reason, it is essential to quantify the impact of data uncertainty propagating through the inversion. Hence, the resulting rigidity distribution is complemented with a sensitivity map.

To determine the model parameters uncertainties, we link the covariance matrices of the parameters and data. Let us consider the observation (geodetic velocities) as a sum of a with zero variance (i.e. and an uncertain quantity :

. For this sum, the covariance matrix is given by: (

Because is deterministic, and are independent (i.e.

. Therefore, the covariance matrix reduces to:

(10)

We consider a linear relationship between and :

(11) which leads to:

(12) where is the Jacobian matrix made of the derivatives of the velocities at geodetic measurement locations with respect to the model parameters. Similarly to (Eq. 9), we define the covariance matrix of the predicted parameters as:

(13) where (the actual values of the parameters) is assumed deterministic. Again, and are assumed independent, and therefore:

(14)
Finally, equation (Eq. 12) becomes:

(15) (16)
So, this equation formulates the uncertainty propagation from geodetic measurements to the model parameters via the Jacobian matrix . The construction of this matrix can be performed in two different ways. The simplest approach consists in expressing it analytically as a function of the gradients that have been evaluated during the resolution of the adjoint problem. Indeed, can be explicitly deduced from the equation . This approach is straightforward and theoretically correct, but it is numerically unstable since it involves the inversion of singular matrices. Consequently, it is more robust to build from finite difference computations. This consists in perturbating one parameter around its optimum (typically by 10%), and then computing the perturbation of the predicted velocity at all geodetic measurement locations. This approach is numerically robust because it involves no matrix inversion. With the second member of Eq. 16 in hand, we can now provide an estimation of the variance (diagonal of the covariance matrix) of the optimization variables.

In this study the parameter is the compliance with its standard deviation . We define dissymmetric upper and lower bounds around the optimum for the rigidity parameter . To represent the uncertainty on the rigidity ( ) we use the fact that is the inverse of the compliance and therefore:

(17)

Determination of effective rigidity for a synthetic case

Before running our optimization scheme on real cases, we evaluate its efficiency to recover a given rigidity distribution (target rigidity) associated with a specific 2-D velocity field .

Surface strain across a locked fault zone can be interpreted either using the concept of a slipping fault zone beneath a locking depth [START_REF] Savage | Geodetic determination of relative plate motion in central California[END_REF] or by assuming a shear-rigidity variation perpendicular to the fault [START_REF] Chéry | Geodetic strain across the San Andreas fault reflects elastic plate thickness variations (rather than fault slip rate)[END_REF]. Differences and similarities between these models are discussed in this latter paper. According to the variable rigidity hypothesis, we define a target given by: (

Where is a non-dimensional rigidity, is the distance to the fault and is a characteristic dimension. Solving force balance within such a plate leads to the following fault-parallel velocity field:

(19)
Therefore, such a velocity distribution is the solution of the spatially variable function of Eq.

18 but can also be associated to a screw dislocation at depth (e.g., [START_REF] Savage | Geodetic determination of relative plate motion in central California[END_REF].

In the case of active fault systems, is generally associated to a physical locking depth which can be estimated using geodesy and seismology. In the case of the San Andreas Fault system, values of d range from 6 to 22 km depending on the location along the fault and the method of determination (e.g. [START_REF] Smith-Konter | Locking depths estimated from geodesy and seismology along the San Andreas Fault System: Implications for seismic moment release[END_REF].

We conduct two tests to verify the ability of the method to retrieve the rigidity distribution given by Eq. 18 for different GPS data sets. We also test different values of (from 2 km to 17 km) in order to generate velocity fields commonly observed on the San Andreas Fault. The specific case of fully-creeping faults ( ~0 km) is discussed in Appendix A, with application to the SAF segment located North of Parkfield. We focus here on the consequences of processing two different spatial distributions of GPS data: (1) evenly spatially distributed and

(2) concentrated near highly strained zones. This corresponds to on site situations (Fig. 3).

Both distributions are made of about 120 GPS velocities vectors.

Several experiments have been conducted to define an optimal mesh size. On the one hand, the computational time is related to the mesh size. The finer is the grid, the longer will be the optimization (about 60 times longer for a mesh 3.6 times finer). On the other hand, the grid needs to be fine enough to capture the variations of the velocity field, notably close to high velocity gradient areas such as the creep zones of the Parkfield segment. Eventually, a spatially adaptive mesh should be implemented. We choose to work with a mean constant spacing of 20 km. This configuration is generally a good compromise between the number of available geodetic measurements, the number of parameters that need to be adjusted and the size of the object we want to study.

For the first case, GPS measurements (black arrows on Fig. 3) are uniformly distributed over the domain, with a constant spacing of 20 km. For the second distribution, we mimic GPS network by producing a velocity field whose spatial density decreases with the distance to the fault. In both cases, the domain is a 200-km square with a 20-km mesh size (394 elements). The dextral strike-slip fault (green line on Fig. 3) has a slip rate of 30 mm/yr and is locked at 10 km depth. The admissible values for the non-dimensional relative rigidity range from 1 to 100 (see Appendices A and B for a discussion of such a choice).

We apply our optimization algorithm to invert the two velocity fields (Fig. 4 and5). In the case of the uniform dataset (Fig. 4 a-f), we first compare the synthetic velocities (red dots on Fig. 4e) to the modelled ones (grey dots on Fig. 4e) along a profile (white dashed line on Fig. 4a) perpendicular to the fault (black dashed line on Fig. 4a-c). The dataset from the inversion almost perfectly matches the characteristic shape of a 2-D arctangent velocity field given by Eq. 19 (Fig. 4e). The misfit between predicted and theoretical (Fig. 4d) permits to estimate the tendency in over or under estimating in our inversions. Besides, the difference between synthetic and modelled velocities, hereafter called residual velocities, is lower than 0.25 mm/yr over the whole domain (Fig. 4c).

Contrary to real cases where the true effective rigidity distribution is unknown, synthetic cases allow for testing the efficiency of our inversion method to retrieve the quadratic rigidity field given by Eq. 18. Fig 4a shows the rigidity distribution over the whole domain of analysis, while Fig. 4b shows the uncertainty distribution map and Fig. 4f focuses along one transect across the fault. We can notice that, as expected, the code predicts a low rigidity zone (90% of the elements ranging between 1 and 3) along a 40 km-wide area centred on the fault.

Also, increases rapidly with the distance to the fault to reach high values (>30) 60 km from the fault. Associated with these rigidity values, we find uncertainties that are very small where rigidity is small but quite high when the opposite occurs (Fig. 4b). This mainly comes from the predominance of the (squared) rigidity term in Eq. 17. This expresses the fact that, in areas that do not deform significantly, very large values of rigidity are admissible (up to infinity) without modifying significantly the local velocity field. Since our search interval for rigidity is bounded, our optimal solutions tend to underestimate the real rigidity in non-deforming areas. This can be seen far from the fault in all the synthetic cases presented in this study (Fig. 4, 5, A.1 and A.2). Finally, we find that, within the uncertainties estimated by our method (paragraph 3.4), our predicted rigidity distribution fits its theoretical value. This is clearly true along the transect crossing the fault on Fig. 4f.

For a data set whose density decreases with distance to the fault (Fig. 5a-f), we observe the same ability for the optimization algorithm to retrieve an arctangent-shaped velocity field (Fig. 5f) and this despite sparse data away from the fault. In this synthetic case, some elements of the grid contain more than one velocity, making the capture of very local velocity gradients difficult if not impossible using one single rigidity parameter over each mesh element. Therefore, residual velocities (Fig. 5c) are generally higher for case 2 than for case 1, with 10% of residual vectors greater than 1 mm/yr (the 1 uncertainty associated with the data being 2 mm/yr) mainly located in the vicinity of the fault where each element of the grid contains several GPS measurements. As a result, local gradients are more difficult to estimate than for case 1 and this could explain the distribution of residual we observe in Fig. 5c.

Despite these moderate residuals, the mean residual velocity over the whole dataset is as low as 0.85 mm/yr which is lower than the 1 uncertainty of the data. This situation is typical of real dataset with a high density of GPS installed in highly deformed areas.

Finally, as for the uniform case, we find that our inversion leads to a distribution of rigidity that fits well its theoretical model within the predicted uncertainties (Fig. 5f). Indeed, considering the 40-km band width around the fault, (Fig. 5a) shows that 66% of the elements show low rigidities (between 1 and 3) while 29% present moderate ones (between 3 and 10).

As described above, in very few deforming areas, the optimal solution underestimates the real rigidity but the uncertainty associated with these high values of rigidity tends to be quite high.

Moreover, the uncertainty values also depend on the local density of geodetic measurements.

Consequently, even when the optimization leads to fairly correct values of low-to-moderate rigidity close to the fault, their uncertainties may be large (Fig. 5b) as one can see along the transect between 20 and 100 km especially if the data distribution is random (Fig. 5f). Again, we present the misfit between and Overall, the satisfactory results of this experience lead us to keep this dimensioning of the grid (triangles with about 20 km edges) for the real case application below.

EFFECTIVE RIGIDITY OF WESTERN USA

Tectonic context and GPS data

The tectonic of the Western United States mostly occurs in response to the relative motion between the Pacific plate and the North American plate. Two main zones accommodating the deformation are the San Andreas Fault system zone and the Basin and Range. In northern California, the relative motion between the Pacific plate and the Sierra Nevada reaches a differential rate of 30 mm/yr and results in large earthquakes. East of the Sierra Nevada, a significant part of the deformation (~10 mm/yr) occurs within the Basin and Range over a broad fault system. To the south, most of the strain is accommodated by the San Andreas

Fault system while the southern Basin and Range is relatively inactive [START_REF] Kreemer | Geodetic constraints on areal changes in the Pacific North America plate boundary zone: What controls Basin and Range extension?[END_REF]. Although significant vertical deformation can occur during seismic events (Landers 1992, Northridge 1994 or Hector Mine 1999 earthquakes, red stars on Fig. 6), vertical motion observed in the area are nearly 10 times smaller than the horizontal velocities during interseismic periods [START_REF] Smith-Konter | Vertical crustal displacement due to interseismic deformation along San Andreas fault: Contraints from tide gauges[END_REF]. Consequently, we chose to analyse only horizontal motion.

We focus our study on the southern part of the San Andreas Fault system (SAFS) where highquality spatially dense GPS measurements are available. We use the CMM3 (Southern California Earthquake Center Crustal Motion Map Version 3.0, SCEC CMM3) velocity field as it was published by [START_REF] Kreemer | Geodetic constraints on areal changes in the Pacific North America plate boundary zone: What controls Basin and Range extension?[END_REF]. It is supposed to represent the interseismic motion that affects our region of interest. This means that all transient motions induced by the seismic events of Landers, Northridge and Hector Mine have been modelled and removed. These data are associated with relatively homogeneous uncertainties of 1.2 mm/yr in average.

A Lambert conformal conic projection is used to project the GPS velocity field on a Cartesian frame. To evaluate the effect of the choice of the domain, we analyse two overlapping areas shown in Fig. 6. We aim at checking that effective rigidity values remain invariant regardless of the chosen borders. The first area of interest, hereafter named Zone 1, is limited by a red dashed line on Fig. 6 and is identical to the one used by [START_REF] Chéry | Plate rigidity inversion in southern California using interseismic GPS velocity field[END_REF]. Then, a

Model parametrization and regularization coefficients

According to the synthetic experiments presented above, we choose a uniform grid spacing of 20 km. This configuration leads to meshes of 2284 elements.

At first, we attempt to evaluate the Tikhonov parameter . To do so, we analyse the trade-off between the normalized regularization member of the functional along the domain boundaries, , and the residual data misfit at all geodetic measurements within the domain (Fig. 7a). Each point of the curve represents an optimization for a given value of the regularization parameter . A decrease of corresponds to an increase of the regularization of the velocity field along the domain boundary, meaning that high gradient changes of are smoothed. This would confer some degree of smoothness to the solution. On the contrary, a reduction of regularization enables a better fit to high velocity gradient changes along the domain boundary. These particularly occur at the transition between highly deforming fault zones and rigid far fields. Nevertheless, this may induce undesirable velocity gradient variations where the velocity field is smooth. Hence, in order to find the appropriate balance between the regularity of the boundary conditions and velocity residual, we compare observed and modelled velocity distributions along the boundary (Fig. 7b-c-d). When the damping parameter is small (Fig. 7b), we allow the regularization term of the functional to be high. This in turn permits to better fit the observed velocities close to the boundaries and consequently within the domain. Nevertheless, the boundary solution may

show in some places a degree of sharpness that is not supported by any data. Increasing the damping parameter used in Fig. 7c increases the regularity of the boundary conditions while still fitting properly the data along the border. We observe that this is done without significantly altering the fit between modelled and observed velocities. Finally, increasing the damping parameter, which means that extremely smooth boundary conditions become admissible, leads to incompatibility between modelled and observed velocities along the border (Fig. 7d). From this analysis, we set the regularization parameter to for which the balance between the regularization of the boundary conditions and the fit to observed velocities within the domain appears to be optimal.

Results of the inversion

Considering the model geometry and parametrization previously described, we perform the inversion of the GPS velocities for the two selected zones (Zone 1 and Zone 2) of the Western United States.

Estimated relative rigidity and corresponding uncertainty distributions

The inversion of the interseismic velocities leads to the distribution of effective rigidity illustrated by Fig. 8a-b. In the case of Zone 1 (Fig. 8a), the lowest values of (1-1.5) are centered on the Mojave Desert, whereas slightly higher rigidities (1.5-4) are observed along the San Andreas Fault system and in the extreme South of the Eastern California Shear Zone.

However, lower values of rigidity (associated with higher deformation rates) are expected along the San Andreas Fault system rather than in the Mojave Desert. We expect that this artefact is likely due to an over-correction of the post-seismic motion of the seismic events of 1992 (Landers), 1994 (Northridge) and 1999 (Hector Mine) within the CMM3 velocity field [START_REF] Liu | Recovery of secular deformation field of Mojave shear zone in southern California from historical terrestrial and GPS measurements[END_REF] . When the GPS data are processed to only keep the interseismic velocity, the post-seismic answer of the earthquakes is estimated at its best. This artefact in our results could help identify the residual post-seismic motion left in the data. As for the high rigidities (>12), they are associated to the South Basin and Range and the South Sierra Nevada where no significant deformation needs to be accommodated. As an extension of Zone 1, the inversion in Zone 2 (Fig. 8b) produces similar rigidity distribution along the SAF and the extreme South of the ECSZ (1.5-4), with, again a surprisingly low rigidity (<1.5) in the Mojave Desert. However, one main difference can be underlined as a zone with rigidity ranging from 6 to 12 is found in the eastern part of the South Basin and Range.

As described in paragraph 3.4, we determine the uncertainties associated with our rigidity estimation which essentially result both from the local measurement density and the uncertainties associated with the data themselves. For each mesh element, we estimate the lower and upper admissible value for rigidity (Fig. 9a-b). First, along all the active fault systems, identified rigidity values are quite low as the amplitude between the upper and lower bounds are lower than 3. The reliability of our solution in deforming zones comes from the local high density of measurements and from large amplitude of the deformation. Conversely, when entering rigid zones, where only few measurements are available, the uncertainties increase very much reaching values that typically range from 2.5 to above 16 by several orders of magnitude. This is notably the case East of the ECSZ. Although the distribution of rigidity shown in Fig. 8a suggests an optimal value of 6-12, uncertainties in this area (Fig. 9)

indicate that a much larger rigidity value (higher than 16 by several orders of magnitude) is also valid. This can be noted in the inversion over the shifted domain (Fig. 8b).

Associated residual velocities

Alongside with the distribution of the rigidity, we evaluate the difference between GPS and the modelled velocities to produce the residual map (Fig. 10a-b). The fit between observed and modelled data is estimated using the normalized root mean square (

) (McCaffrey, 2005), ( 20
)
where and stand for the eastern and northern directions respectively, is the residual velocity, the data standard error and N the number of data. In addition to the , the weighted root mean square ( ) gives a measure of the a posteriori weighted scatter in the fits [START_REF] Mccaffrey | Block kinematics of the Pacific-North America plate boundary in the southwestern United States from inversion of GPS, seismological, and geologic data[END_REF],

For Zone 1 (Fig. 10a), we get a of 1.26, with a of 1.10 mm/yr that can be compared with the uncertainty of 1.20 mm/yr associated with the data. The highest residuals (>4.5 mm/yr) occur on the southern segments of the SAF, while intermediate residuals ( 2.5-4.5 mm/yr) are unevenly distributed between high and low data density zones.

A similar analysis for Zone 2 gives a of 1.25 (Fig. 10b) with a of 0.93 mm/yr.

The difference observed in the values of both zones can be explained by the data distribution.

Indeed, the second zone excludes some of the velocities that are poorly estimated by the optimization (notably on the Pacific plate) and includes few vectors that are better recovered.

DISCUSSION

Based on the hypothesis that interseismic strain mostly reflects rheological contrasts across the lithosphere, the solved inverse problem entirely depends on the quality of the CMM3 velocity field. Therefore, we first discuss the sensitivity of the model result with respect to the data. Then, we discuss our results (strain rate and rigidity distributions) in the light of the ones provided by previous studies on western US and California. We finally discuss the future use of our method for tectonic and geodynamic purposes.

Robustness of the inversion

In this study, we use the entire dataset of the Southern California Crustal Motion Map Version 3.0, involving 615 vectors for Zone 1 and 530 vectors for Zone 2. In order to evaluate the impact of data selection on rigidity distribution, we perform complementary inversions using identical parametrization, but removing GPS vectors whose residual norms are greater than a given threshold value. These residuals can be due to three different factors:

1) GPS uncertainty. Data uncertainties range from 0.16 mm/yr to 3.71 mm/yr for the horizontal components with a value of 1.20 mm/yr. The maximum data uncertainties are observed in the South of Mojave Desert, along the Los Angeles Bay and for a few isolated points in the Sierra Nevada and ECSZ.

2) Local motions. Besides interseismic plate motions, some sites may be affected by gravitational collapse, geothermal activity (e.g. [START_REF] Vasco | Geodetic Imaging: High Resolution Reservoir Monitoring Using Satellite Interferometry[END_REF] or the exploitation of aquifer systems (e.g. [START_REF] Galloway | Detection of aquifer system compaction and land subsidence using interferometric synthetic aperture radar, Antelope Valley, Mojave Desert, California[END_REF][START_REF] Hoffmann | Seasonal subsidence and rebound in Las Vegas Valley, Nevada, observed by synthetic aperture radar interferometry[END_REF]. Many different processes can locally obscure the GPS interseismic velocity component, such as unravelled postseismic motions.

3) Modelling. The optimization algorithm and the forward model can also be at the origin of the residual velocities. Indeed, a poor estimation of the velocity along the boundaries could be the reason why high residual velocities are observed at the junction between the fault and the boundaries of the domain for both synthetic and real data cases. Furthermore, our forward model includes several assumptions such as an absence of body forces. Also, the data are assumed free of post-seismic effects which can be inexact if all post-seismic effects due to the Landers, Northridge and Hector Mine earthquakes, for instance, have not been fully removed.

We choose to withdraw from 0.3% ( 6 mm/yr) up to 50% ( 1.3 mm/yr) of the data in order to analyse the stability of the solution of our inversion. The corresponding , and the correlation of the rigidity distribution with respect to the one obtained by the previously described inversion are gathered in Error! Reference source not found.. We choose experiment 1 as the reference solution to estimate the rigidity correlation. Removing up to 10% of the GPS measurements exhibiting significant residuals in our initial inversion (experiences 3 to 6 in Table 1) neither enhance the or the , nor significantly modify the rigidity distribution. But considering only 50% of the data (experience 7 in Table 1) improves the to 0.67 and the to 0.54 mm/yr and leads to a rigidity correlation of 0.831 with the main features preserved.

This approach echoes the strategy [START_REF] Meade | Block models of crustal motion in southern California constrained by GPS measurements[END_REF] developed to reduce the number of stations and therefore to minimize the uncertainty magnitude. While based on different quality criteria, they remove about 50% of the initial dataset (CMM3) to compute their inversions. Our results illustrate that selecting the data with the lowest residuals does not significantly influence the modelled rigidity (see correlation in Table 1). However, in areas where data density is poor, a reduction of 50% can lead to a completely different interpretation. Therefore, keeping the whole dataset seems preferable.

Strain rate: comparison with other approaches

Most of strain rate computations derived from GPS velocity measurements stand on a continuous approximation of a model velocity field. A simple way to compute the strain rate is to design a triangulation of the GPS points collection and then assume that the velocity field inside each triangle evolves linearly. However, this method generates a non-smooth strain map due to a linear interpolation of measured GPS velocities. This method can be adapted to areas covered by sparse GPS networks [START_REF] Masson | Seismic versus aseismic deformation in Iran inferred from earthquakes and geodetic data[END_REF] but generates erroneous strain rates when applied to dense networks such as the ones installed in California. In this case, a smooth approximation of the velocity field needs to be performed in order to avoid spurious strain rate modeling. Consequently, a suitable method must also account for high strain gradient occurring around fault zones. A large variety of mathematical approaches can be used to deduce a strain rate map, often leading to relatively large differences [START_REF] Feigl | Space Geodetic Measurement of Crustal Deformation in Central and Southern California , 1984-1992[END_REF]Mc Caffrey et al., 2005;[START_REF] Shen | Crustal deformation across and beyond the Los Angeles basin from geodetic measurements[END_REF][START_REF] Tape | Adjoint tomography of the southern California crust[END_REF].

Our optimal solution of rigidity distribution can be used for the determination of strain rates over the whole study area. But, we have shown above that our models systematically underestimate rigidity in very few deforming areas, typically far from the active fault systems.

This bias is partly counterbalanced by the information provided by the upper bound of the admissible rigidity values. These latter are very close to the optimal solution in deforming zones, while they suggest that a purely rigid behavior may be considered when the deformation is very small, even though a slight deformation remains admissible just considering geodetic measurements way to conform to geological considerations and block-model assumptions that state that, in most cases, far from the faults, the blocks are rigid. So, we used the upper bound rigidity distribution (Fig. 9b) for creating our strain rate map (Fig. 11a).

We compare in Fig. 11 our strain rate map (through the 2 nd invariant of the strain tensor) with the one obtained by a method originally proposed by [START_REF] Haines | A Procedure for Obtaining the Complete Horizontal Motions Within Zones Deformation of Strain Rate Data[END_REF] and later revised in the framework of the strain map global project [START_REF] Kreemer | A geodetic platemotion and Global Strain Rate Model[END_REF]. Although both methods depend on distinct assumptions, they produce similar intensities (> 64 nanostrain/yr) located near faulted areas along the SAF and the ECSZ. This overall similarity is probably due to the fact that both approaches are able to produce a low residual between the discrete and the continuous velocity fields. Using our strongest admissible rigidity solution leads to low strain estimates in weakly deforming areas that are similar to the ones obtained by the global strain map project. This can be noticed in the Great Valley between the SAF and the Sierra Nevada and along the Pacific coast.

A significant difference between the two strain rate maps can be found only on two limited areas: offshore the Pacific coast and east of the ECSZ. Because these two areas display low residuals (Fig. 10), we guess that our model is likely not able to locally estimate the strain rate precisely. This could be due on the one hand, to an improper estimate of the boundary conditions notably within the Pacific plate, and on the other hand, to a very low local data density. Indeed, whereas [START_REF] Kreemer | A geodetic platemotion and Global Strain Rate Model[END_REF] only interpolate the strain rate dataset to best fit the data, our solution aims at doing the same, but under the constraint of the stress equilibrium equation (Eq. 1). As demonstrated by the synthetic benchmarks presented in paragraph 3.4, evenly distributed data lead to a better estimation of the rigidity. Therefore, a future use of our methodology could be to invert interpolated GPS velocities (such as the ones provided by the Global Strain Rate Project) instead of the original GPS data to compute effective rigidity distribution at a continental scale.

Lastly, we compare the spatial distribution of our dilatational strain rate solution with the one obtained by [START_REF] Kreemer | A geodetic platemotion and Global Strain Rate Model[END_REF] (Fig. 12). We use the first invariant of the strain rate tensor (mean of its trace) as a first-order approximation of the dilatational strain rate.

Neither the strain compatibility approach used by [START_REF] Kreemer | A geodetic platemotion and Global Strain Rate Model[END_REF] nor our study, take vertical velocity measurements into account. Nevertheless, the plane stress formalism of our modelling leads to the prediction of vertical strain rates, which is not the case in [START_REF] Kreemer | A geodetic platemotion and Global Strain Rate Model[END_REF] analysis. Yet, recent analyses (e.g., [START_REF] Becker | Western US intermountain seismicity caused by changes in upper mantle flow[END_REF] suggest that the rate-change of vertical loading of the lithosphere may play a dominant role in defining the distributions of seismicity and therefore strain.

The study of [START_REF] Lowry | The role of crustal quartz in controlling Cordilleran deformation[END_REF] provides a map of the flexural elastic thickness ( ) for the entire western US. We assume that a linear relationship exists between the in-plane plate rigidity and its corresponding thickness (Chery, 2008 and present work).

Therefore, our map of is directly proportional to the distribution of shown in Fig. 8a.

Although such a linear relationship is valid only if elastic parameters do not vary with depth, it provides a simple way to estimate the effective elastic thickness for our modelling. For the purpose of comparison with Lowry and Pérez-Gussinyé ( 2011), we display their value of over Zone 1 (Fig. 13). Flexural and geodetic elastic thicknesses displayed in Fig. 13 show a very limited degree of agreement. For example, the flexural thickness map predicts a thick plate for most of the SAF, while a low geodetic elastic thickness is deduced using the interseismic velocity field. The only area suggesting some resemblance corresponds to the Basin and Range around the ECSZ and the SAF around the Salton Trough for which both methods display low elastic thickness. In order to find some justifications about the large discrepancies between and at least two lines of arguments could be investigated.

First, despite the formal similarity between flexural plate and shear plate theories [START_REF] Chéry | Plate rigidity inversion in southern California using interseismic GPS velocity field[END_REF], they may reflect two distinct lithospheric behaviours. For example, as stated by [START_REF] Thatcher | Temporal evolution of continental lithospheric strength in actively deforming regions[END_REF], plate flexure is the result of a long term loading over millions of years, implying that the strain rate in most of the lithosphere is close to zero. is a measure of stress that is supported dynamically over very long timescales by a lithosphere that is in a state of frictional failure and viscoelastic flow, meaning the strain rate is virtually zero.

However, given the shorter timescale of geodetic observation and the clear evidence for seismic release of significant elastic strain potential accumulated on century timescales, likely does predominantly reflects the elastic behavior of a thicker domain associated to interseismic deformation. Another difference may come from the lithospheric loading.

Vertical loads modify distinct components of the strain tensor. Indeed, those induce flexure and plate motions and therefore horizontal shear. Hence, distinct behaviours may emerge from these kinds of load.

In the brittle part of the crust, background seismicity is likely to reflect the loading of interseismic motion, therefore introducing an anelastic component into the analysed shear motion. Beneath the crust and especially under a shear zone like the San Andreas Fault system, the upper mantle presents a laterally variable and strong anisotropy [START_REF] Hartog | Depth-dependent mantle anisotropy below the San Andreas fault system: Apparent splitting parameters and waveforms[END_REF]. If such anisotropic behaviours occur at both crustal and mantle levels, flexural and horizontal loading may activate two different rheological systems that could result into significant differences in terms of effective elastic thickness.

A second way to investigate is to assume that flexural and geodetic thicknesses represent the same mechanical concept. However, they could be differently revealed by the data because of the formal differences between the two inverse problems. In the case of flexural thickness, the determination of is based on the correlation between topographic and gravimetric signal.

Among other factors, erosion can smooth or sharpen the topographic signal. Even if its influence can be accounted for in modelling approaches (e.g. [START_REF] Forsyth | Subsurface Loading and Estimates of the Flexural Rigidity of Continental Lithosphere[END_REF], the impact of erosion on the determination of seems difficult to quantify due to large uncertainty associated to past erosion. In addition, a geodynamical setting mostly involving shear motion may not be adapted at all for a flexural plate analysis because such a motion is not likely to produce neither topographic nor gravimetric signals. Last but not least, inverse theory of plate flexure requires that flexural thickness cannot be determined for resolutions smaller than the characteristic flexural wavelength [START_REF] Watts | Isostasy and flexure of the Lithosphere[END_REF]. This also may explain why a sharp rigidity variation across the SAF cannot be resolved by this method. Even if our methodology has never been used prior to [START_REF] Chéry | An integrated mechanical model of the San Andreas Fault in central and northern California[END_REF], the direct relation between shear strain and shear rigidity is likely to produce high resolution estimate of geodetic thickness for zones where the geodetic strain is well defined. Conversely, we acknowledge that our uncertainty analysis predicts inaccurate rigidity determination in zones of low strain-rate like the Sierra Nevada.

Also, lithospheric loads like body forces and basal stress coming from mantle motion can impact the strain-rate field and therefore altering the determination of the shear rigidity. The identification of the importance of such effects must be tackled by future studies.

In order to better understand the discrepancy between flexural and shear analysis, a tractable way would be to design a complete mechanical model of western US as it was done for example by [START_REF] Pollitz | Viscoelasticcycle model of interseismic deformation in the northwestern United States[END_REF]. Such a model could be used to predict synthetic topographic, gravimetric and deformation datasets obeying to momentum and constitutive equations. Then lithospheric deformation and compared to the rheological input of the forward model.

CONCLUSION

A global inversion strategy has been proposed for the identification of effective rigidity maps using GPS velocity fields under minimum a priori assumptions. Taking advantage of the selfadjoint nature of the governing equations, large dimensional problems coming from necessary high resolution distribution of the rigidity have been considered. Compared to the previous study carried out by [START_REF] Chéry | Plate rigidity inversion in southern California using interseismic GPS velocity field[END_REF], the results are now backed by uncertainty analysis which suggests that the effective rigidity can only be accurately determined in moderate or highly strained areas.

This is a high-resolution methodology which can be seen as a mechanical model to link shear rigidity to interseismic strain with no prior knowledge of fault locations. The main limitation of this approach relies to the plane stress hypothesis used in the forward model. Therefore, no strain variation occurs with depth for a given horizontal location over the plate. This behaviour is probably over simplified around active faults acting like screw dislocations as proposed by [START_REF] Savage | Geodetic determination of relative plate motion in central California[END_REF]. To complete what is presented here, the following directions can be considered:

1) The 2D effective rigidity model can be replaced by a 3D model of Western United

States including the effective elastic thickness as the main geophysical parameter.

Because this approach would include the full 3D strain rate tensor, it would provide a more realistic approximation of the plate behaviour of the lithosphere especially around faults.

2) The 2D approach can be used over wide areas, for instance at the continental plate scale, after a splitting in patches. This would permit to determine large scale rigidity maps in the framework of the global strain map project of [START_REF] Kreemer | A geodetic platemotion and Global Strain Rate Model[END_REF].

3) The strong spatial correlation between low rigidity areas and active fault zones also suggests that our methodology could be applied for deciphering active faults in tectonically poorly known areas.

similar solution where ground deformation is significant, and to a larger ratio between rigidity in non-deforming zones relatively to deforming areas. However, both solutions (using a 2-fold or 4-fold range of magnitude) fit the data globally in the same way. This confirms the fact that a very large range of rigidity is admissible far from the active fault systems. This statement is illustrated again by the lower and upper bound solutions (Fig. B.1b-c). The rigidity distribution is well constrained in significantly deforming zones, lower and upper bounds being close to the optimal solution, but badly constrained elsewhere. Typically, any value of rigidity higher than 10 is an admissible solution for rigid areas.

In conclusion, we see that almost all the significant ground deformation is captured using a 2orders of magnitude range for relative rigidity. Increasing this range allows either a better modelling of ground deformation induced by shallow surface creeping segments, or the assignment of higher rigidity values in non-deforming zones. But, in this latter case, the fit to the velocity field is very little improved and we know that the uncertainties associated with these rigidity values remain very high.

  translation moves Zone 1 by 100 km towards the Northeast to obtain the second region called Zone 2 (blue dashed lines on Fig. 6). Both areas include the central San Andreas Fault system (SAFS) segment, the Eastern California Shear Zone (ECSZ), the south Sierra Nevada (SN) to the North, the Mojave Desert (MD) in the centre, the Salton Sea (SS) and the south Basin and Range (SBR) to the East. The western part of Zone 1 contains a part of the Pacific Plate along the Californian coast whereas Zone 2 is directly bounded by the San Andreas fault to the West.

Fig. 1 .

 1 Fig. 1. Schematic representation of the optimization problem: the domain is meshed with elements of constant rigidity and submitted to Dirichlet boundary conditions along the boundary . The black arrows symbolise the geodetic measurements within the domain, the blue ones the Dirichlet conditions which are part of the optimization variables and the red ones are the solution of the model.

Fig. 2 .

 2 Fig. 2. Sketch of optimization algorithm applied to plate rigidity inversion. For each iteration we optimize both the rigidity within the domain and the velocity along the boundaries.

Fig. 3 .

 3 Fig. 3. Distribution of the synthetic velocities considering an evenly sampled domain (black arrows) and a randomly sampled domain whose density decreases with the distance to the fault (red arrows). The fault (green line) is a dextral strike slip fault locked during the interseismic motion.

Fig. 4 .

 4 Fig. 4. Results of the inversion for the synthetic case with a geodetic spatial sampling (black dots on (a), (c) and (d)) whose density is constant whatever the distance to the fault (black line at coordinate 0 along the Y-axis). (a) rigidity distribution determined by the optimization; (b) associated uncertainty; (c) norm of residual velocities; (d) misfit between predicted and theoretical ; (e) velocity measurements (black squares with uncertainty bars) and velocities predicted by our model (red circles) along the profile shown in white dotted line on (a); (f) rigidity values and their associated uncertainties along the same profile. The plotted values are estimated at the barycentre of the elements of the mesh in a 30-km wide bandwidth centred on the profile. The colour code is the same as in (a).

Fig. 5 .

 5 Fig. 5. Results of the inversion for the synthetic case that mimics a real geodetic spatial sampling (black dots on (a) (c) and (d)) whose density decreases with the distance to the fault (black line at coordinate 0 along the Y-axis). (a) rigidity distribution from the optimization; (b) associated uncertainty; (c) norm of residual velocities; (d) misfit between predicted and theoretical ; (e) velocity measurements (black squares with uncertainty bars) and velocities predicted by our model (red circles) along the South to North profile shown in white dotted line on (a); (f) rigidity values and their associated uncertainties along the same profile. The opaque rectangle over the northern termination of rigidity and velocity profile (e) and (f) highlights a zone where the density of measurements is low, leading to high rigidity uncertainties.

Fig. 6 .

 6 Fig. 6. Spatial distribution of the geodetic measurements on the SAF system used for the inversion. The black arrows show the velocity field in the North American reference frame. Our main domain of analysis is shown in the red rectangle (Zone 1 615 GPS velocities) while the blue one (Zone 2 530 GPS velocities) represents a translation of the area of interest. The red stars indicate the location of Landers (L), Northridge (N) and Hector Mine (HM) earthquakes.

Fig. 7

 7 Fig. 7 Pareto curve for different regularization parameters used in the optimization algorithm applied to the southern California. (a) Plot of the normalized velocity variations norm as a function of the residual data misfit as damping varies. (b), (c) and (d) represent the velocities observed within a 10-km distance of the domain boundary (green) and calculated on the boundaries (red) for different damping parameter .

Fig. 8 .

 8 Fig. 8. Distributions of the relative rigidity D represented for Zone 1 (a) and for Zone 2 (b). According to Eq. 23, the elastic thickness is proportional to . The GPS data are represented by black dots, the faults by green lines and the Landers (L), Northridge (N) and Hector Mines (HM) earthquakes by white stars. SAF, SN, ECSZ, SBR, MD and SS stand for San Andreas Fault, Sierra Nevada, Eastern California Shear Zone, South Basin and Range, Mojave Desert and Salton Sea respectively.

Fig. 9 :

 9 Fig. 9: (a) Lower and (b) upper bounds for the rigidity values around the optimal distribution displayed on Fig. 8a for Zone 1. In very few (or non-) deforming areas, reach valuesthat are several orders of magnitude higher than the optimal solution.

Fig. 10 .

 10 Fig. 10. Norm of residual velocities between GPS and modelled velocities associated with the rigidity distribution determined for a) Zone 1 ( =1.26) and b) Zone 2 ( =1.25).

Fig. 11 .

 11 Fig. 11. Distribution of the strain rate for Zone 1. Second invariant of the strain rate tensor from (a) the upper (strongest) admissible values of our rigidity optimization and (b) the global strain map of Kreemer et al. (2014).

Fig. 12 .

 12 Fig. 12. Distribution of the first invariant of the strain rate tensor for Zone 1 from (a) our optimization and (b) the global strain map of Kreemer et al. (2014).

Fig. 13 .

 13 Fig. 13. (a) Geodetic elastic thickness associated to our study; (b) flexural elastic thickness given by Lowry and Pérez-Gussinyé (2011).

Figure A. 1 :

 1 Figure A.1: Velocity and relative rigidity distributions for synthetic cases mimicking shallow creep on strike-slip faults. The locking depth is set to 2 km. The search interval for the relative rigidity has 2 orders (left) or 4 orders (right) of magnitude. (a and d) Theoretical velocities (green dots) and modelled velocities (black dots) at measurement locations along the profile (dotted black line) indicated in (c). (b and e) Theoretical (dashed curve) and modelled relative rigidity (colored dashes using same color palette as in (c) and (f) respectively) along the same profile. (c and f) Spatial distribution of relative rigidity determined by our inversion.

Figure A. 2 :

 2 Figure A.2: Same as Fig. A.1 but the locking depth is now set to 0 km.

Figure A. 2 :

 2 Figure A.2: (a) Relative rigidity distribution over the Parkfield segment of the SAF. Black squares indicate the location of the cities of San Juan Bautista (north) and Parkfield (south). White arrows are the velocities of the CMM3 database. Measured (red circles) and modelled (black circles) along-strike velocities along (b) A-B and (c) C-D profiles shown in (a).

Figure B1 :

 B1 Figure B1: Rigidity distribution over Southern California as determined using a 4-order of magnitude range of admissible relative rigidity (from 1 to 10000). (a) Optimal solution, (b) Difference with respect to the optimal solution using a 2-order magnitude range (Fig. 8a) (c) lower bound solution and (d) upper bound solution. To be compared with Figures 8 and 9.

  .

  The code does not incorporate the value of the plate thickness. Therefore, we cannot discriminate plate thickness and elastic properties of the lithosphere from the rigidity values. Hence, fixing an absolute value to the rigidity remains an open problem. Finally, the velocity field produced by the forward model depends on two input parameter sets: the velocity boundary conditions and the

	effective rigidity distribution	as	where	is the minimum value of
	over the domain. So, all distributions of	presented in this paper range from 1 to some
	maximum values.
	For a given spatial domain, we generate a uniform 2-D Delaunay mesh composed by

The Young modulus remains the only free mechanical parameter in this equation. Since the model is driven by a velocity condition, only . This means that any distribution of the form provides the same velocity field regardless the value of the constant C. For this reason, we define the non-dimensional triangles. In order to estimate the velocity field at geodetic measurement locations we use the academic 2-dimensional finite element code CAMEF.

Table 1 .

 1 List of the experiments, the / values of the residual velocities and the correlation of the rigidity distribution relatively to experiment 1.

	Experiment Zone 1 1 2 2 3 1 4 1 5 1 6 1 7 1	Grid Size 20 20 20 20 20 20 20	Number of GPS vectors 615 530 613 (99.7%) 600 (97.5%) 583 (94.7%) 553 (90%) 307 (50%)	of the residual velocities 1.26 1.25 1.25 1.21 1.16 1.10 0.67	(mm/yr) 1.10 0.93 1.08 1.04 1.00 0.93 0.54	Rigidity Correlation 1 0.831 0.995 0.986 0.987 0.958 0.831

ACKNOWLEDGEMENTS

We thank Prof. Riad Hassani from Nice-Sophia Antipolis University for making available to us his plane stress Finite Element code CAMEF. We also acknowledge Tony Lowry and an anonymous reviewer for their constructive remarks and suggestions. The PhD of S. Furst is supported by the Total Company and the LabEx NUMEV project (n° ANR-10-LABX-20) fund French National Research Agency (ANR).

Despite the difference in their estimation, both spatial distributions of the dilatation strain rate from [START_REF] Kreemer | A geodetic platemotion and Global Strain Rate Model[END_REF] and us are very similar. They notably highlight the compressive context of the SAF system along the central bend. The only noticeable difference can be found along the fault system located north of Los Angeles where vertical motion is known to occur along active thrust faults (e.g. Northridge or Compton faults).

Rigidity of the lithosphere and effective elastic thickness

In the following, we study the relation between in-plane rigidity associated with geodetic strain (this work) and the flexural rigidity deduced from gravity and topographic data analysis [START_REF] Audet | Dominant role of tectonic inheritance in supercontinent cycles[END_REF][START_REF] Lowry | The role of crustal quartz in controlling Cordilleran deformation[END_REF][START_REF] Tesauro | Ductile crustal flow in 265[END_REF]. In the case of a thin curved elastic plate, the relation between the bending moment and the flexural rigidity is given by: (

where is the vertical displacement of the plate and R(x) its local curvature radius (e.g. [START_REF] Turcotte | Geodynamics[END_REF]. Using Eq. 2, a horizontal force per unit area applied to a vertical section of the lithosphere can be defined as:

where is the stiffness of the lithosphere to horizontal strain. is equal to where is the shear modulus (Pa) and the plate thickness (m). Therefore, in-plane rigidity is expressed in N while the flexural rigidity is given in Nm, precluding a direct comparison between these two fields. In order to compare our relative rigidity map with the flexural rigidity deduced from gravity and topographic data analysis [START_REF] Audet | Dominant role of tectonic inheritance in supercontinent cycles[END_REF][START_REF] Lowry | The role of crustal quartz in controlling Cordilleran deformation[END_REF][START_REF] Tesauro | Ductile crustal flow in 265[END_REF], we use the elastic thickness associated to these two formalisms.

Appendix A) Application to shallow creeping faults

In this appendix, we briefly discuss the way our approach deals with very high strain rates and what this implies in term of prior search interval for the relative shear-rigidity values.

Three terms contribute to the spatial regularization (smoothing) of the modelled distribution of relative rigidity and, consequently, the modelled velocity field. The first one simply comes from mesh size and, therefore, from the spatial resolution of our model. Then, a Tikhonov regularization term can be activated in the cost function (Eq. 5). Finally, a spatial smoothing process is performed on the gradient of the functional (Fig. 2).

All these components of our approach are likely to limit the range of strain rates that can be properly modelled. In particular, active faults exhibiting very shallow creep are supposed to lead to local under-estimation of rigidity gradients.

In practice, the mesh size is a few kilometres (20 km in the case of this study over South California), the Tikhonov regularization term is not active, and the smoothing process over the gradient of the functional can have a very limited spatial extent. Despite the limitations of the effects of these regularization components, we see that the whole range of the expected strain rates (about 4 orders of magnitude) can hardly be captured by our approach.

Nevertheless, we demonstrate, hereafter, that (1) if the spatial sampling is high, then using a 4-order of magnitude search interval allows for a good fit to sharp velocity changes. ( 2) Using a 2-order of magnitude search interval is sufficient in most cases, small but significant local residues being present only in the case of extremely shallow creeping behaviour.

First, we run our inversion on synthetic cases that are similar to those presented in section 3, but for which the locking depth is much shallower: 2 km (Fig these experiments, the mesh size is set to 3km.

In the case of a 2-km locking depth, we find that exploring a 2-order of magnitude interval for rigidity is sufficient to perfectly model the velocity field (Fig. A.1). The only significant difference between the modelled and the theoretical values of the relative rigidity appears on non-significantly deforming zones at distances higher than about 20 km from the fault.

Indeed, in these areas, we know that changing the relative rigidity by several orders of magnitude has very small influence on the velocity prediction. So, again, we show that our approach underestimates rigidity far from the active faults. Exploring a 4-order of magnitude rigidity interval does not improve the fit to ground velocities. It only pushes a little bit away the limit where the real rigidity is underestimated.

If we consider, now, the extreme case of a fault experiencing creep up to the ground surface, then the strain rate across the fault is infinite. A similar analysis to the preceding case shows that significant residues remain close to the fault (Fig. A.2). Yet, their amplitudes are very small (no more than 2 mm/yr) except exactly on the fault itself. In such an extreme case, indeed, extending the search interval of rigidity from 2-to 4-order of magnitude allows for the proper modelling of the velocity field everywhere. Nevertheless, this better fit has been obtained at the cost of the regularity of the rigidity distribution which is slightly altered in non-deforming zones (upper and lower central areas).

segment is known to be experiencing very shallow creep (e.g. [START_REF] Rosen | Aseismic creep along the San Andreas Fault northwest of Parkfield, Ca measured by radar interferometry[END_REF]. We use a 10-km mesh size and run inversions varying the range of the admissible rigidity values and the smoothing factors for the gradient of the cost function. All lead to the same solution that is presented on Fig. A.3. Using a 4-order of magnitude search interval does not provide additional information that would not have been captured by a 2-order of magnitude search interval. This may be due to the 10-km resolution of our mesh, as well as the relatively low density of GPS measurements in the CMM3 database, notably with respect to other highresolution measurements like InSAR.

In order to illustrate the efficiency of our modelling, we compare the modelled velocity field with the CMM3 measurements along two profiles across the SAF system, one (A-B) just south of the city of San Juan Bautista (Fig Appendix B) Processing the SAF zone using a 4-orders of magnitude range for admissible rigidity

In the main body of this paper, we determine the relative rigidity distribution (both for synthetic and real cases) from the exploration of a 2-orders of magnitude range of rigidity values. Appendix A shows that, unless we need to deal with sharp velocity gradients that are typically found on very shallow creeping fault segments, this rigidity range is sufficient to properly fit the interseismic velocity field within their uncertainties. In this appendix, we

show that using a wider range of admissible rigidity value (4 orders of magnitude instead of 2) in the inversion process applied to Southern California leads to a very similar solution.

As noticed in Appendix A, the inversion becomes now longer and less steady.


Laurent Meersseman

A NOTE ON THE AUTOMORPHISM GROUP OF A COMPACT COMPLEX MANIFOLD

In this note, we give explicit examples of compact complex 3-folds which admit automorphisms that are isotopic to the identity through C ∞ -diffeomorphisms but not through biholomorphisms. These automorphisms play an important role in the construction of the Teichmüller stack of higher dimensional manifolds.

Introduction.

Let X be a compact complex manifold and M the underlying oriented C ∞ manifold. The automorphism group Aut(X) of X is a complex Lie group whose Lie algebra is the Lie algebra of holomorphic vector fields [START_REF] Bochner | Groups on analytic manifolds[END_REF]. We denote by Aut 0 (X) the connected component of the identity. Its elements are thus automorphisms f such that there exists a C ∞ -isotopy (1.1) t ∈ [0, 1] -→ f t ∈ Aut(X) with f 0 = Id and f 1 = f.

Note that Aut(X) has at most a countable number of connected components so the quotient Aut(X)/Aut 0 (X) is discrete.

Let Diff(M ) be the Fréchet Lie group of C ∞ -diffeomorphisms of M . It is tangent at the identity to the Lie algebra of C ∞ vector fields. Let Diff 0 (M ) be the connected component of the identity. Its elements are C ∞diffeomorphisms f such that there exists a C ∞ -isotopy (1.2) t ∈ [0, 1] -→ f t ∈ Diff(M ) with f 0 = Id and f 1 = f.

Note that the discrete group Diff(M )/Diff 0 (M ) is the well known mapping class group. Define now (1.3) Aut 1 (X) := Aut(X) ∩ Diff 0 (M ).

There are obvious inclusions of groups

(1.4) Aut 0 (X) ⊆ Aut 1 (X) ⊆ Aut(X)
In many examples, the first two groups are the same but differ from the third one (think of a complex torus). The purpose of this note is to describe an explicit family 1 of compact 3-folds X a,b such that Hence Aut 1 (X a,b ) has a connected components, and this number can be chosen arbitrarily large. Our main motivation, which is detailed in Section 2, comes from understanding the Teichmüller stack of M , that is the quotient stack of the set of complex operators on M modulo the action of Diff 0 (M ).

The construction of the manifolds as well as the computation of their automorphism groups are elementary. As often when looking at explicit examples, the crux of the matter was to find the idea that makes everything work. We asked several specialists but they did not know any such example. We tried several classical examples but it always failed. Finally, we came accross the good family when looking for deformations of Hopf surfaces over the projective line P 1 in connection with a different problem. The manifolds X a,b are such deformations with the following additional property. All the fibers are biholomorphic except for those that lie above 0 and above an ath root of unity. Every automorphism must preserve these special fibers so must project onto P 1 as a rotation of angle 2πk/a for some k. This explains the a connected components of the automorphism group. Finally, diffeomorphically there is no special fiber since X a,b is just a bundle. It is then not difficult to check that all rotations are allowed for diffeomorphisms so every automorphism is in Diff 0 (M ).

Motivations

Thanks to Newlander-Nirenberg Theorem [START_REF] Newlander | Complex analytic coordinates in almost complex manifolds[END_REF], a structure of a complex manifold X on M is equivalent to an integrable complex operator J on M , that is a C ∞ bundle operator J on the tangent bundle T M such that (2.1)

J 2 = -Id and [T 0,1 , T 0,1 ] ⊂ T 0,1
where T 0,1 is the subbundle of the complexified tangent bundle T M ⊗ C formed by the eigenvectors of J with eigenvalue -i.

It is easy to check that the complex manifolds X J := (M, J) and X J ′ are biholomorphic if and only if there exists a diffeomorphism f of M whose differential df commutes with J and J ′ . This defines an action of Diff(M ).

The Teichmüller space T (M ) of M is then defined as the quotient of the space I(M ) of integrable complex operators on M (inducing the orientation of M ) by Diff 0 (M ).

As such, this is a topological space whose quotient by the mapping class group is the moduli space of complex structures on M . If M is a surface, then this is the classical Teichmüller space and a complex manifold in a natural way. In higher dimensions, this is usually not even locally an analytic space, cf. [START_REF] Meersseman | The Teichmüller and Riemann Moduli Stacks[END_REF]Examples 11.3,11.6].

Kodaira-Spencer and Kuranishi classical deformation theory (see [START_REF] Kodaira | Complex Manifolds and Deformations of Complex Structures[END_REF] for a good introduction) provides each compact complex manifold X with an analytic space, its Kuranishi space, which encodes all the small deformations of its complex structure. This is not however a local moduli space, but must be thought of as the best approximation in the analytic category of a local moduli space. In particular, the same complex structure may be encoded in an infinite number of points.

In other words, there is a surjective map from the Kuranishi space K J of X J onto a neighborhood of J in T (M ) but which is in general far from being bijective (cf. [START_REF] Catanese | A Superficial Working Guide to Deformations and Moduli[END_REF] and [START_REF] Catanese | Topological methods in moduli theory[END_REF] where several results where equality holds are discussed). Especially, Aut 0 (X J ) acts2 on K J identifying equivalent complex structures [START_REF] Meersseman | Foliated Structure of The Kuranishi Space and Isomorphism of Deformation Families of Compact Complex Manifolds[END_REF], [8, §2.3]. The point here is that the dimension of Aut 0 (X J ) is only an upper semi-continuous function of J so may jump. When this occurs, the previous action is non-trivial and the Teichmüller space is not a local analytic space around J.

We are not finished yet. The quotient of K J by Aut 0 (X J ) is still not a priori a neighborhood of J in T (M ). The natural projection map may still have non-trivial discrete fibers. This phenomenon is thoroughly studied in [START_REF] Meersseman | The Teichmüller and Riemann Moduli Stacks[END_REF]. It is shown that the action of Diff 0 (M ) onto I(M ) is a foliation transversely modelled onto the stack [K J /Aut 0 (X J )] and that a holonomy groupoid can be defined for such an object. The holonomy morphisms describe exactly the non-trivial fibers (see in particular Remarks 10.4 and 10.8). Hence, holonomy measures the gap between the stacks T (M ) (locally at J) and [K J /Aut 0 (X J )]. In many cases however, there is no holonomy. This leads to the following problem Holonomy Problem. Give examples of Teichmüller spaces with nontrivial holonomy. This is the main motivation of this paper. Note that the isotropy group of a point J in I(M ) is exactly Aut 1 (X J ) and not Aut 0 (X J ). So letting (M, J) encoding one of the manifolds X a,b , our result says that this isotropy group is not connected. As a consequence, the Teichmüller space of M has non-trivial finite holonomy at J.

3.

The manifolds X a,b .

Let a and b be two nonnegative integers. Let λ be a non-zero complex number of modulus strictly less than one. For further use, we define the following two surfaces. Let X 0 be the Hopf surface defined as C 2 \ {(0, 0)} divided by the group generated by the contraction (z, w) → (λz, λw). Let X 1 be the Hopf surface defined as C 2 \{(0, 0)} divided by the group generated by the contraction (z, w) → (λz + w, λw). These two Hopf surfaces are not biholomorphic, cf. [START_REF] Kodaira | Complex Manifolds and Deformations of Complex Structures[END_REF].

We consider the vector bundle O(b)⊕O(a) → P 1 . Throughout the article, we make use of the charts

(3.1) (t, z 0 , w 0 ) ∈ C 3 and (s, z 1 , w 1 ) ∈ C 3
subject to the relations

(3.2) st = 1, z 1 = s b z 0 , w 1 = s a w 0 .
Let c > 0 and let σ be a non-zero holomorphic section of O(c). In accordance with (3.1) and (3.2), we represent it in local charts by two holomorphic maps σ 0 and σ 1 satisfying σ 1 (s) = s c σ 0 (t). Let W be O(b) ⊕ O(a) minus the zero section.

Lemma 3.1. Assume that b -a -c ≥ 0. Then, the holomorphic maps

(3.
3) (t, z 0 , w 0 ) → g 0 (t, z 0 , w 0 ) = (t, λz 0 + σ 0 (t)w 0 , λw 0 ) and

(3.4) (s, z 1 , w 1 ) → g 1 (s, z 1 , w 1 ) = (s, λz 1 + s b-a-c σ 1 (s)w 1 , λw 1 )
defines a biholomorphism g of W .

Proof. Just compute in the other chart

g 1 (s, z 1 , w 1 ) =(s, λz 1 + s b-a-c σ 1 (s)w 1 , λw 1 )
=(1/t, s b (λz 0 + σ 0 (t)w 0 ), s a (λw 0 )) so g 0 and g 1 glue in accordance with (3.2).

Consider now the group G = g . It acts freely and properly on W and fixes each fiber of W → P 1 . The quotient space W/G is thus a complex manifold. More precisely Proposition 3.2. The manifold W/G is a deformation of Hopf surfaces over P 1 . Moreover the fiber over t ∈ P 1 is biholomorphic to X 0 if t is a zero of σ, otherwise it is biholomorphic to X 1 .

In particular, W/G is compact.

Proof. We already observed that the bundle map W → P 1 descends as a holomorphic map π : W/G → P 1 . It is obviously a proper holomorphic submersion, hence it defines W/G as a deformation of complex manifolds parametrized by the projective line. The fiber over t is C 2 \ {(0, 0)} divided by the contracting map (z, w) → (λz + σ(t)w, λw). If t is a zero of σ, then this is exactly the Hopf surface X 0 . Otherwise, it is biholomorphic to X 1 , see [START_REF] Kodaira | Complex Manifolds and Deformations of Complex Structures[END_REF]. for t ∈ C.

For the rest of the paper, we assume that a is strictly greater than 3. Observe that the condition b ≥ 3a is nothing but b -a -c ≥ 0.

Computation of the automorphism groups.

We are in position to state and prove our main result.

(4.1) Aut 0 (X a,b ) ≃ α P 0 α | α ∈ C * , P ∈ C b-a [X] /G and (4.2) Aut(X a,b ) =Aut 1 (X a,b ) ≃G a × α P 0 α | α ∈ C * , P ∈ C b-a [X] /G
where G a is the group of a-th roots of unity, C b-a [X] is the space of complex polynomials with degree at most b -a and where the product in (4.2) is given by

(4.3) r, α P 0 α • r ′ , β Q 0 β = rr ′ , αβ αQ • r + βP • r 0 αβ
We note the immediate corollary Proof. Choose a fiber of X a,b → P 1 isomorphic to X 1 . Restrict f to it and compose with the projection onto the projective line. This gives a holomorphic map from X 1 to P 1 , hence a meromorphic function on X 1 . But the algebraic dimension of X 1 is zero, see [START_REF] Dabrowski | Moduli Spaces for Hopf Surfaces[END_REF], so this map is constant. In other words, f sends the π-fibers isomorphic to X 1 onto the π-fibers. By density of these fibers, f sends every π-fiber onto a π-fiber so descends as an automorphism h of P 1 .

Lemma 4.4. The automorphism h is a power of the rotation at 0 of angle 2π/a.

Proof. Note that f must send a π-fiber biholomorphic to X 0 onto a π-fiber biholomorphic to X 0 . Now the set of such fibers is the set of a-th roots of unity plus zero by Proposition 3.2 and (3.5). It follows from Lemma 4.3 that the automorphism h is an automorphism of the projective line which preserves this set. Since a > 3, it must preserve at least three different points of the unit circle, hence must preserve the unit circle. But this implies that zero is fixed. Schwarz Lemma shows now that it is a power of the rotation 3 at 0 of angle 2π/a.

Lift f as an automorphism F of the universal covering W of X a,b . We denote by (F 0 , F 1 ) its expression in the charts (3.1). Lemma 4.5. In the charts (3.1), the lifting F has the following form

(4.4) F 0 (t, z 0 , w 0 ) = (r k t, αz 0 + τ 0 (t)w 0 , αw 0 ) and (4.5) F 1 (s, z 1 , w 1 ) = (r -k s, r -kb (αz 1 + τ 1 (s)w 1 ), αw 1 )
where r = exp(2iπ/a), k is an integer, α a complex number and τ = (τ 0 , τ 1 ) is a section of O(b -a).

Proof. The first coordinate in (4.4) comes from Lemma 4.4. For the two other coordinates, recall from [START_REF] Wehler | Versal deformation of Hopf surfaces[END_REF] that the automorphism group of X 0 is GL 2 (C) (modulo quotient by the group generated by the contraction) and 3 As pointed out to us by F. Bosio, this is no more true for a = 3. Letting j = exp(2iπ/3), the map z → -(z -j)/(2jz + j 2 ) preserves the set but is not a rotation.

that of X 1 is the group of upper triangular matrices with both entries on the diagonal equal (modulo quotient by the group generated by the contraction).

Hence the general form of F 0 is (4.6) F 0 (t, z 0 , w 0 ) = (r k t, α 0 (t)z 0 + τ 0 (t)w 0 , α 0 (t)w 0 ) for α 0 and τ two holomorphic functions. But in the other chart, using the same more general form of (4.5), we must have

F 1 (s, z 1 , w 1 ) =(r -k s, r -kb (α 1 (s)z 1 + τ 1 (s)w 1 ), α 1 (s)w 1 ) =(1/(r k t), (r -k s) b (z 0 α 1 (1/t) + τ 1 (1/t)s a-b w 0 ), s a α 1 (1/t)w 0 )
which extends at s = 0 and glues with (4.6) if and only if α = (α 0 , α 1 ) is a constant and τ = (τ 0 , τ 1 ) is a section of O(b -a).

It remains to check whether these automorphisms really descend as automorphisms of X a,b , i.e. whether they commute with the contraction g of Lemma 3.1. We compute (3.5). A similar computation holds in the (s, z 1 , w 1 )-coordinates, so finally all these automorphisms descend. Lemma 4.6. An automorphism of X a,b is in the connected component of the identity if and only if it descends as the identity of P 1 .

g 0 • F 0 (t, z 0 , w 0 ) = (r k t, λαz 0 + λτ 0 (t)w 0 + ασ 0 (r k t)w 0 , λαw 0 ) = F 0 • g 0 (t, z 0 , w 0 ) since σ is G a -invariant, cf.
Proof. If an automorphism f of X a,b is in the connected component of the identity, then by Lemma 4.4 its projection h is isotopic to the identity through rotations of angle 2iπk/a. This is only possible if h is the identity. Conversely, if h is the identity, it is easy to see that in Lemma 4.5 we can move α to 1 and τ to the zero section and obtain a path of automorphisms from f to the identity. Lemma 4.7. Every automorphism of X a,b is isotopic to the identity through C ∞ -diffeomorphisms.

Proof. Let B a,b be the bundle over P 1 with fiber X 0 obtained by taking the quotient of W by the group generated by the λ-homothety in the fibers. Observe that X a,b can be deformed to B a,b putting a parameter ǫ ∈ C and considering the family

(4.7) W × C/ g -→ C
where the action is given by (we just write it down in the first chart):

(4.8) (t, z 0 , w 0 , ǫ) → g0 (t, z 0 , w 0 , ǫ) = (t, λz 0 + ǫσ 0 (t)w 0 , λw 0 , ǫ) Hence X a,b is C ∞ -diffeomorphic to B a,b .
More precisely, let X ǫ a,b be the fiber of the family (4.7) over ǫ. Then X 0 a,b = B a,b , and X 1 a,b = X a,b and by Ehresmann's Lemma there is an isotopy of C ∞ -diffeomorphisms φ ǫ from X ǫ a,b onto X 1 a,b with φ 1 equal to the identity. Let f be an automorphism of X a,b . It is easy to check, using Lemma 4.5, that f is still an automorphism of B a,b and of all the X ǫ a,b . We are saying that the map (4.9)

F (t, z 0 , w 0 , ǫ) = (F ǫ (t, z 0 , w 0 ), ǫ) = (f (t, z 0 , w 0 ), ǫ)
is an automorphism of the whole family which induces f on each fiber. The isotopy

φ ǫ • F ǫ • (φ ǫ ) -1 joins the automorphism f = F 1 of X a,b and φ 0 • F 0 • (φ 0 ) -1 through C ∞ -diffeomorphisms.
But now, at ǫ = 0, since B a,b is a holomorphic bundle, we may take any rotation at 0 as map h and thus construct a path of automorphisms g t of B a,b between F 0 and the identity. Combining the isotopy φ 0 • g t • (φ 0 ) -1 with the previous one, we obtain that the automorphism f of X a,b is isotopic to the identity through C ∞ -diffeomorphisms.

The proofs of Theorem 4.1 and Corollary 4.2 follow easily from the previous lemmas.

Final Comments.

The construction behaves well with respect to pull-backs of the base manifold P 1 . In particular, if g is a ramified covering from P 1 to P 1 , the pull-back manifold Y = g * X a,b will have a more complicated finite group Γ as quotient Aut 1 (Y)/Aut 0 (Y). In this way, we may construct examples with Γ being any cyclic or dihedral group, or with Γ being the tetrahedral, octahedral or icosahedral group. 4 There are several open questions left on this topic. First the manifolds X a,b are not even Kähler and it would be interesting to have a similar example with projective manifolds. Also, it would be interesting to have examples with Aut 0 (X) reduced to zero, but Aut 1 (X) not, especially examples of surfaces of general type, cf. [START_REF] Catanese | A Superficial Working Guide to Deformations and Moduli[END_REF], [START_REF] Catanese | Topological methods in moduli theory[END_REF] and the subsequent literature. Finally, the most exciting would be to find an example with Aut 1 (X)/Aut 0 (X) infinite since it would give a Teichmüller space with infinite holonomy at some point. As pointed out to us by S. Cantat, this cannot happen for Kähler manifolds, since the kernel of the action of Aut(X) on the cohomology contains Aut 0 (X) as a subgroup of finite index [START_REF] Lieberman | Compactness of the Chow scheme: applications to automorphisms and deformations of Kähler manifolds[END_REF]. But in the non-Kähler world everything is possible.

For all these additional questions (except perhaps for the first one), it seems that a really different type of examples is needed.

(1. 5 )

 5 Aut 0 (X a,b ) Aut 1 (X a,b ) = Aut(X a,b ) and (1.6) Aut 1 (X a,b )/Aut 0 (X a,b ) = Z aDate: April 10, 2017. 1991 Mathematics Subject Classification. 32M05,32G05 .1 It is indexed by two integers satisfying b ≥ 3a and a > 3.

1

 1

Definition 3 . 3 .

 33 Assume that c = 2a and that b ≥ 3a. We denote by X a,b the manifold W/G corresponding to the choice (3.5) σ 0 (t) = t a a-1 k=0 (t -exp(2iπk/a)).

Theorem 4 . 1 .

 41 The manifold X a,b satisfies

Corollary 4 . 2 .Lemma 4 . 3 .

 4243 The group Aut 1 (X a,b ) has a > 3 connected components and the quotient Aut 1 (X a,b )/Aut 0 (X a,b ) is isomorphic to the cyclic group Z a .Theorem 4.1 will be proved through a succession of Lemmas. Let f be an automorphism of X a,b . Then it respects π and descends as an automorphism h of P 1 .

This is not exactly an action, cf. [8, §2.3] for details.