# Plate rigidity inversion in southern California using interseismic GPS velocity field

Jean Chery, B. Mohammadi, Michel Peyret, Claire Joulain

► To cite this version:

Jean Chery, B. Mohammadi, Michel Peyret, Claire Joulain. Plate rigidity inversion in southern California using interseismic GPS velocity field. Geophysical Journal International, 2011, 187 (2), pp.783-796. 10.1111/j.1365-246X.2011.05192.x. hal-00669897

HAL Id: hal-00669897

https://hal.science/hal-00669897v1

Submitted on 11 Jun 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Plate rigidity inversion in southern California using interseismic GPS velocity field

J. Chéry,<fnref n="1" /> B. Mohammadi,<fnref n="2" /> M. Peyret<fnref n="1" /> and C. Joulain<fnref n="1" />

Accepted 2011 August 10. Received 2011 August 9; in original form 2010 July 23

## SUMMARY

This paper presents an inversion method using the interseismic velocity field to determine effective rigidity of the lithosphere. The method is based on the minimization of a cost function defined as the quadratic measure of the difference between measured and modelled velocity fields on a discrete set of points. The continuous mapping of the rigidity is fulfilled with a limited set of parameters and the forward solution is achieved using a plane stress finite element code. The computation of the cost function gradient in the parameters' space allows one to iteratively find the best parameters set through a suitable optimization algorithm.

We first design a benchmark including an abrupt rigidity variation that cannot be described by a continuous function. For such a case, we show that increasing the number of parameters is a way to accurately describe sharp variations of the rigidity map. Then, we use a dense GPS velocity field over the southwestern United States to estimate the corresponding rigidity variations for different spatial resolutions of the parameters' grid. We analyse the conceptual and practical difficulties associated with our methodology. Finally, rigidity maps obtained by our inversion method in southwestern United States and particularly across the San Andreas Fault System are reviewed and compared to current plate rigidity estimates and geophysical data over this area.

**Key words:** Inverse theory; Seismic cycle; Dynamics of lithosphere and mantle; Rheology; crust and lithosphere.

## 1 INTRODUCTION

The interseismic velocity field in actively deforming domains is traditionally modelled by a buried dislocation aiming to represent the long-term fault motion at depth (Savage & Burford 1973). Because it has been shown that most of these geodetic slip rates are fairly consistent with rates deduced from geological measurements (Reilinger et al. 2006), this approach is becoming paramount. However, this way to model interseismic strain is not unique and we adopt in this paper a different viewpoint. Indeed, we propose a mechanical relationship between the plate rigidity and the interseismic strain distribution that can be measured by various geodetic methods (continuous or campaign mode global positioning system (GPS), InSAR and leveling). This subject has been previously addressed for a 1-D horizontal simple shear situation (Chery 2008). Our purpose here is to show that the spatial distribution of the plate rigidity can be estimated in 2-D by inverting a dense velocity field data set with a suitable optimization method.

### 1.1 Flexural lithosphere rigidity

Large geological strain of the lithosphere mostly reflects irreversible stress–strain relations involving viscosity, plasticity and friction.

However, numerous topographic and gravimetric features of the oceanic and continental lithosphere have been interpreted like the flexure of a plate of effective elastic thickness (often called $T_e$) subjected to vertical and horizontal loads (Watts 2001). In oceanic domains, topographic profiles around oceanic volcanoes and bulges associated with plate curvature in subduction zones can be explained using simple flexural plate models (Fig. 1a). Loads can be distributed over the plate (case of a sedimentary basin) or point like (case of a volcano) and bending moment can be also added at plate ends (Weissel & Karner 1998). This formalism is quite appropriate to study the elastic behaviour of the lithosphere because plastic and frictional strains are generally small.

However, geological strain is often much higher than elastic strain due to plate flexure. For example, topography undulations in the central Indian Basin are associated with pervasive faulting that can be explained as an elasto-plastic buckling of the lithosphere submitted to horizontal forces (Gerbault 2000). Deciphering elastic strain in continental areas is still more challenging. Indeed, present-day topography results from the whole geological and geomorphologic history of the area that can largely affect the plate rigidity estimate. For this reason, the rigidity parameters obtained by adjusting an elastic model to best simulate topographic and gravimetric data must be termed as ‘effective’ as they are likely to incorporate the signature of non-elastic strain.

<footnote n="1">¹Université Montpellier 2, Géosciences Montpellier, CNRS UMR-5243, 34095 Montpellier, France. E-mail: jean.chery@gm.univ-montp2.fr</footnote>

<footnote n="2">²Université Montpellier 2, I3M, CNRS UMR-5149, 34095 Montpellier, France</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** (a) Vertical motion $u_z$ induced by a vertical load applied on one side of an elastic plate with a variable rigidity. (b) Horizontal motion $u_y$ induced by a horizontal load applied on the same plate. The asthenosphere resists passively to plate motion in both cases.</figcaption></figure>

Two methods have been used so far to compute effective rigidity of a bended lithospheric plate. When the loads can be geographically defined (such as for a volcano or a mountain belt), a 1-D or 2-D plate model can be used to invert its elastic properties. However, different studies over the same geographical areas may display large discrepancies of the inverted rigidities (Hetenyi *et al.* 2006). When loads appear to be distributed inside the lithosphere, it becomes more appropriate to use the correlation between topographic and gravimetric fields to find the plate rigidity explaining at best the observed signals (e.g. McNutt 1990; Audet & Bürgmann 2011).

## 1.2 Plate rigidity at different timescales

Lithospheric loads occurring at short timescales like earth tides reveal a purely elastic behaviour of the crust and the mantle. In contrast, permanent vertical loads like mountains or volcanoes lead to a complete relaxation of the underlying asthenosphere. Between these extrema, transient strain is observed, therefore revealing the time-dependent behaviour of the crust and the mantle (Thatcher & Pollitz 2008). In such a case, a relevant rheological model seems to be an elastic plate (the lithosphere) lying over a viscoelastic asthenosphere (Nur & Mavko 1974).

Transient viscoelastic behaviour occurs over timescales of 1–100 yr as attested by post-seismic motion (Ergintav *et al.* 2007) and motion rebound associated with large lakes artificially filling/infilling (Kaufmann & Amelung 1995). Inverted linear viscosities range from $10^{19}–10^{21}$ Pa s for the lower crust to $10^{18}–10^{19}$ Pa s for the upper mantle. At longer timescales (10–100 kyr), variation of loads induced by glaciers and large lakes induce a transient motion that can be measured geodetically (Milne *et al.* 2001) or by analysing abandoned shorelines (Bills *et al.* 1994). These latter authors found viscosities as low as $4 × 10^{17}$ Pa s at 40 km depth in the eastern Great Basin (western United States) and a viscosity of $2 × 10^{20}$ Pa s for the uppermost mantle. Concerning the loading timescale effect on the rheological response of the lithosphere, the studies above reveal that long-term loading results in a lower estimate of elastic plate thickness than observations of strain relaxation over much shorter times (Thatcher & Pollitz 2008). Although the measurements of these transient motions are useful to perform a joint inversion of plate rigidity and its underlying viscosity, the places where post-seismic or postglacial strain occurs are few and far between. Therefore, these observations do not allow for an extensive rigidity mapping on the continents.

## 1.3 Interseismic strain and shear plate rigidity

Contrary to post-seismic or postglacial loading that can be idealized like the response to a step function leading to a stress relaxation with

decaying velocities, interseismic strain appears like a kinematic loading associated with plate motion. If a measured geodetic strain is not affected by the coseismic or post-seismic motion of a previous large earthquake, it may therefore be representative of the slow loading phase of the lithosphere during the seismic cycle (Fig. 1b). Up to now, interseismic geodetic strain provided by GPS and InSAR methods have been mostly used to constrain kinematic models using thick rigid lithospheric blocks connected by slipping faults at depth, usually called block models (Meade & Hager 2005). Therefore, it is possible to invert a velocity field to obtain both blocks motion (translation and rotation on the sphere) and fault motion (slip rate and locking depth). In such a case, the velocity field across a strike-slip fault follows an arctangent variation associated with the solution of a screw dislocation in an elastic half-space (Chinnery 1963).

Mechanically speaking, the block model is made of infinitely rigid blocks connected by the locked part of the fault above the locking depth. We propose to modify the block model approach in two ways:

(1) We assume that the lithosphere rigidity can continuously vary across fault zones and rigid blocks. In this way, we define a continuous version of the block model with a laterally variable rigidity (as for flexural modelling) where the integrated shear rigidity of the lithosphere combines both the thickness and the intrinsic elastic parameters of the plate. We consider here a simple case for which interseismic shear only occurs on vertical planes. This leads to a strain setting for which the velocity field does not change with depth. Therefore, integrated rigidity is simply the product of the vertical average of the intrinsic shear rigidity in the plate with the plate thickness.

(2) The motion at the surface of a given portion of the lithosphere is driven by lateral forces applied to the boundary of the considered domain. This assumption is quite different to the one made by the block model that prescribes internal boundary conditions. This way of modelling the interseismic velocity field is close to the approach used in geomechanical modelling: the observed geological strain is modelled using a mechanical model with assumed rheology and boundary conditions (Chéry *et al.* 2001).

We still make some additional assumptions:

(1) No strong interaction occurs between the plate and the underlying mantle, as it occurs, for example, in subduction zones. Moreover, we assume that the viscous coupling with the asthenosphere is small during the interseismic phase. This can be justified as follows: mantle viscosities deduced from post-seismic motion and postglacial rebounds are on the order of $10^{18}–10^{19}$ Pa s as (see Subsection 1.2). If we assume a passive mantle–plate interaction in the zone of interest, strain rate at the base of the lithosphere has the
order of magnitude of the interseismic lithospheric strain, that is, $10^{-17}–10^{-14} s^{-1}$. The resulting interacting stress is therefore on the order of $10^1–10^5$ Pa, much smaller than the deviatoric stress acting in the lithosphere. According to this reasoning, we conjecture that the mechanical setting of a plate during interseismic loading is close to what occurs during the flexural rebound of the same plate submitted to a geological load. However, as these two mechanisms are not related to the same geophysical processes, we call *H* the thickness associated to the shear rigidity of the lithosphere to differentiate it from the flexural thickness (*T*<sub>e</sub>).

(2) We assume that the surface motion is equal to the motion at depth. We acknowledge that this simplification may be not true everywhere since local strain rotation may, for example, be associated with fault tips. However, we conjecture that this kind of strain pattern is not dominant and that it is therefore reasonable to model the plate deformation using a pure shear approach (i.e. main strain axes are horizontal and vertical). In our case, we model lithospheric strain using a 2-D plane stress approximation for which the relevant material parameter is the integrated shear rigidity, called *D* hereafter. Although this approach allows for the prediction of the vertical motion associated with horizontal strain, we did not compare it with geodetic vertical motion since tectonic motion in California is essentially horizontal.

Even with these simplifying assumptions, the inversion of the effective rigidity of the continental lithosphere based on the knowledge of the interseismic velocity field at the surface is not straightforward. Indeed, there is no direct relation between the geodetic strain observed at one location and the rigidity of this location. Rather, the geodetic strain at this location depends in an unknown manner on the whole rigidity distribution within a certain neighbourhood and on the imposed velocity field along the boundaries of the study domain. Moreover, the uniqueness of the solution is not warranted since two different rigidity distributions could lead (at least in principle) to the same velocity field. A last difficulty is linked to the relation between the number and spatial distribution of geodetic measurements and the optimal number of parameters needed to estimate the rigidity. In Section 2, we propose a way to deal with these difficulties using a global optimization approach.

The paper is organized as follows: we first present the methodology of inverse modelling with the chosen optimization method. A synthetic case is built to test the impact of rigidity parametrization on the solution quality. Then an application to southern California is done using the dense GPS velocity field of this region. Finally, the significance of the rigidity map of southern California is discussed in light of the geophysical data.

# 2 METHODOLOGY OF INVERSE MODELLING

## 2.1 Introduction

Because the physical properties of the Earth interior cannot be directly deduced from direct observation, inverse problems based on Earth's surface data are at the heart of many geophysical problems (Tarantola 1987). In the case of the inversion of geodetic data during interseismic strain, the parametrization of the problem is largely influenced by the nature of the forward model.

If the forward model is purely kinematical, it means that the predicted velocity field does not derive from a physical conservation equation. Therefore, the inverse problem is equivalent to data fitting. This approach is well illustrated by the construction of the strain

map using geodetic data (Tape *et al.* 2009), with the joint use of other data constraints like fault slip orientations and earthquake focal mechanisms (Kreemer *et al.* 2003). The choice of the forward model is largely dependent on the type of approximation function for the velocity field. Polynomial, splines or other continuous and derivable functions like wavelets can be used, thus allowing an easy computation of the strain map. The optimal solution is generally found using a least-squares minimization.

Another use of interseismic strain is the estimation of long-term fault slip rate and locking depth of crustal faults. In the case of an infinitely long single strike-slip fault, only these two parameters are searched (Savage & Burford 1973). However, most of intracontinental domains like the western United States are characterized by changing fault azimuth and dip angle. Also the expected fault slip rate along faults is expected to vary. A generalization of the Savage and Burford approach is therefore a model made of a collection of blocks separated by faults (locked during the interseismic period). The inverse problem is dedicated here to find the fault slip rate distribution at depth (beneath the locking depth) that fits at best the observed velocity. Motion of the blocks is simultaneously inverted (Meade & Hager 2005). As this approach seems to be well suited to invert fault slip rates around vast undeformed domains (Reilinger *et al.* 2006), the uniqueness of the inverse problem is not guaranteed when a collection of small blocks is used as it happens in California: a large trade-off between fault slip rate and locking depth prevents an unequivocal slip rate determination (D'Alessio *et al.* 2005).

We discuss in the following subsections the optimization scheme of the inverse problem, the forward model associated with the mechanical problem, the parametrization of the rigidity distribution, the choice of the cost function and the resulting numerical algorithm.

## 2.2 Global optimization and solution of boundary value problems

We aim to determine plate rigidity parameters that permit to best adjust a known surface velocity field inside a domain $\Omega$ bounded by a frontier $d\Omega$ (Fig. 2). The principle of this method is based on the minimization of a cost function $J$, which represents a measure of the discrepancy between the velocity field $v$ predicted by a forward model and a target velocity $v*$ associated with a discrete geodetic data set. The minimization of this cost function, which convexity is unknown, leads to a spatial distribution of rigidity, which should explain at best the velocity field within $\Omega$.

Since the cost function is not necessarily convex, it likely has multiple local minima. Consequently, a global optimization strategy has to be implemented. Several strategies may be chosen. Some algorithms carry out a direct search (derivation-free) in a multidimensional parameter space. The most widely used are simulated annealing algorithms (e.g. Kirkpatrick *et al.* 1983), genetic algorithms (e.g. Holland 1975), or more recently the neighbourhood algorithm (e.g. Sambridge 1999). They make it possible for the cost function to deteriorate under certain conditions to avoid any local optimum. However, such approaches generally become time-consuming when the dimension of the parameter space increases (typically more than 20). An alternative approach consists in using gradient-based methods within a global search strategy. This is accomplished by using, sequentially or in parallel, various different start points for a standard gradient descending algorithm (e.g. Ugray *et al.* 2007).

In this study, we use a global optimizer of this latter type, which combines a standard gradient-based method within a second
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Spatial view of the optimization problem corresponding to three layers. The data layer is made of velocity measurements included in the domain $\Omega$ (black arrows). The parameter layer defines rigidity values over a grid (green dots). The finite element layer corresponds to a mesh having variable rigidity $D(x,y)$ from one element to another with boundary conditions defined on $d\Omega$. This formalism allows for the prediction of velocities at exact data points location.</figcaption></figure>

minimization problem whose objective is to provide an initial start that belongs to the attraction basin of the global minimum of the cost function (Mohammadi & Saiac 2003). The minimization is driven by a global optimization algorithm. We briefly explain how we introduce global search features into a classical, typically gradient based, local minimization algorithm. Details on the mathematical background and applications of this approach to academic and industrial problems can be found in Mohammadi & Saiac (2003), Ivorra (2006), Ivorra & Mohammadi (2006) and Mohammadi & Pironneau (2009).

Classical gradient-based minimization algorithms can be seen as discrete forms of initial value problem for a first-order dynamical system. For instance, $p$ being a vector of parameters to adjust starting from its initial value $p^0$, an iterative minimizing process is defined by ($p^0$ given)

$$p^{n+1} = p^n - \rho M^{-1}(p^n) \frac{dJ(p^n)}{dp}, \qquad (1)$$

where $\rho$ is a scalar (either fixed or optimized at each step of the iterative process) and $M(p)$ may be the Identity matrix in the case of the steepest descent method or the Hessian matrix of $J$ (respectively, an approximate of the Hessian matrix) in the case of the Newton (respectively, quasi-Newton) method. Eq. (1) is a discrete form of the following first-order differential problem:

$$M(p)p' = -\frac{dJ}{dp}, \quad p(t=0) = p^0, \qquad (2)$$

where $p'$ represents the first derivative of $p$. Let suppose the global minimum is unique. Such algorithms converge towards the global minimum if the initial condition belongs to its attraction basin. Otherwise, the minimizing sequence of parameter $p'$ leads to a local minimum. In that sense, the problem of global minimization with a gradient-based algorithm becomes the prescription of an initial condition for the mentioned initial value problem in the suitable attraction basin. This objective leads to consider this global optimization as a boundary value problem. Let us suppose that the value of infimum $J_m$ of the functional is known, but not where it is reached. Then, we can define the initial condition as an additional

unknown $p(0) = v$ that should minimize the misfit $h(u)$ between the local minimum and the expected global minimum. Hence, our global optimization problem can be written as follows:

$$M(p)p' = -\frac{dJ}{dp}, \quad p(t=0) = v, \quad v = \arg\min_{u}(h(u)). \qquad (3)$$

Now we have two nested minimization problems. If one uses a steepest descent method for the outer minimization problem, this leads to two coupled dynamical systems.

$$M(p)p' = -\frac{dJ}{dp}, \quad p(t=0) = v, \qquad (4)$$

$$v' = -\frac{dh}{dv}, \quad v(t=0) = v^0 \qquad (5)$$

which we solve using the algorithm described in Mohammadi & Saiac (2003).

## 2.3 Forward modelling and basic assumptions in our case

We consider a 2-D deformable plate in the framework of linear and isotropic elasticity. Velocity boundary conditions (also called Dirichlet conditions) are applied all along its boundary $d\Omega$ as shown in Fig. 2. Because the plate traction at the surface corresponds to a free normal traction, a plane stress assumption is used. We do not consider the vertical plate motion associated with its internal strain. The plate has a thickness $H$ and the stress tensor $\sigma$ can be written like

$$\sigma_{ij} = E \left[ \frac{\mu}{1 - \mu^2} \varepsilon_{kk} \delta_{ij} + \frac{1}{1 + \mu} \varepsilon_{ij} \right], \qquad (6)$$

where $\varepsilon$ is the strain tensor and where $i, j$ and $k = 1, 2$. We assume in the following that the Poisson's ratio $\mu$ is constant and equal to 0.25, which is close to its seismological value. The only free mechanical parameter is then the Young's modulus $E(x,y)$. Because the model is driven by a velocity condition, only the relative variation of the Young's modulus matters for strain computation, that is to say that a distribution $C \times E(x,y)$ provides the same velocity field $v$ whatever the value of the constant $C$. For this reason, we use in the following an adimensional rigidity distribution $D(x,y)$ defined by $E(x,y)/E_{min}$ where $E_{min}$ is the minimum value of $E$ found on the domain $\Omega$. Consequently, a non-unit scale is used throughout the paper. The link between the relative rigidity $D$ and its geophysical interpretation in terms of effective lithospheric thickness is treated in the discussion section.

This direct problem is solved with the finite element (FE) method (Zienkiewicz et al. 2005), using the academic code CAMEF (courtesy of Riad Hassani). The rigidity value of each element of the FE mesh is set depending on the model parameters estimated from an independent set of points, called thereafter 'control points'.

## 2.4 Parametrization of the rigidity distribution

Mapping a continuous function based on control points can be done by a variety of interpolation or approximation methods. Since we do not have a priori information on the shape and the wavelengths contained in the rigidity distribution, we choose to define the continuous rigidity distribution $D$ as the interpolation of some scalar values $p_i$ (our model parameters) located at np control points $P_i$,
using a characteristic length scale $d_0$ (eq. 7).

$$D(x, y, d_0) = \frac{\sum_{i=1}^{np} p_i \cdot \exp(-d_i^2/d_0^2)}{\sum_{i=1}^{np} \exp(-d_i^2/d_0^2)} \quad (7)$$

with

$$d_i = \sqrt{(x - x_i)^2 + (y - y_i)^2}, \quad (8)$$

where $(x, y)$ are the coordinates of the current point and $(x_i, y_i)$ the coordinates of the ith control point. Knowing the barycentric coordinates of each element, the above formulae allow for the computation of element rigidity on the whole FE mesh. In this work, we choose to set a square parameter grid, the length scale $d_0$ being equal to the distance between two adjacent control points. In this way, a coarse grid of parameters on the domain $\Omega$ decreases the control point numbers $np$ in the optimization process. A denser parameter grid allows for a better description of the rigidity distribution as $d_0$ decreases. One may note that no relation is given between the parameter grid length scale and the mesh resolution. However, the element size should be smaller than $d_0$ to capture the rigidity variation defined by eqs (7) and (8).

## 2.5 Choice of the cost function

In a continuous framework, a cost function $J$ is a measure of the distance between two continuous fields $v$ and $v*$. In the case where $v$ is a solution provided by the forward analysis and $v*$ is a target solution, $J$ can be defined as a residual estimate with the following $L^2$ integral:

$$J = \frac{1}{|\Omega|} \int_{\Omega} (v - v^*)^2 d\Omega. \quad (9)$$

In practical applications in geodynamics, $v*$ is given by a discrete collection of velocity vectors in places where geodetic measurements are acquired. This leads to redefine $J$ as a discrete sum of velocity discrepancies weighted by a user-defined spatial data density.

$$J = \sum_{i=1}^{n} (v_i - v_i^*)^2 \cdot w_i^2 / \sum_{i=1}^{n} w_i^2. \quad (10)$$

Density weighting coefficients $w_i$ are given to balance the spatial heterogeneities of the GPS measurements. We set the weight on each velocity vector proportional to the minimum distance between this observation and the closest observation on the domain. We increase this way the influence of the areas where GPS benchmarks are few and far between.

If we include the covariance matrix of the data $C$ in the cost function, $J$ can be rewritten in a discrete formulation as

$$J = (v - v^*)^T W^{-1} C^{-1} W (v - v^*). \quad (11)$$

This equation is equivalent to eq. (10) if $C$ is taken as an identity matrix. Using this notation, $W$ is now a diagonal matrix made of the weighting coefficients for all GPS measurements.

## 2.6 Numerical algorithm

To summarize this section, the inverse problem applied to geodetic data is based on a global minimization algorithm. From a known plate geometry, boundary conditions and a given initial rigidity

distribution, a plane stress FE model is used to calculate the cost function $J$ associated with a rigidity distribution which is controlled by a parameter grid. The cost function is minimized to find the model that best predicts the measured velocities. The loop ends when the cost function is found to have decreased enough (this is based on a user choice) or when the variations of the parameters are found negligible. Of course, in both cases there is no guarantee that the solution is a global minimum. The flow of calculation is shown in Fig. 3. In the following sections, we will present the square root of the cost function $(J_{min})^{1/2}$ expressed in mm yr$^{-1}$ to provide a simple comparison with the velocity data.

## 3 SYNTHETIC CASE: A HARD INCLUSION IN A SOFT MATERIAL

Before addressing the problem of inversion of the lithosphere rigidity using real geodetic data, we test the capability of the optimization algorithm to retrieve a given rigidity distribution $D^*$ from its associated velocity field $v^*$. More specifically, we focus on the influence of the spatial resolution of our models.

We attempt to invert a discontinuous rigidity distribution $D^*$ made up of two domains having constant rigidities (Fig. 4). The rigidity of the polygonal inclusion is 10 times higher than that of the surrounding medium. The model is submitted to simple shear on its lateral sides at a rate of 20 mm yr$^{-1}$. The lower and upper sides are let free. The target velocity field $v^*$ is computed with the FE code CAMEF. The 196 nodes of the FE grid are used to define the target velocity field $v^*$. The inverse processing aims to retrieve the rigidity parameters on a regular grid of control points. In this case, it is important to note that all sides of the forward model are constrained to have the velocity corresponding to the desired solution $v^*$. Only the inner nodes of the mesh are treated as unknown. At each iteration, the rigidity value within each mesh element is computed according to eqs (7) and (8).

Due to the discontinuity of the target solution, the parametric rigidity field cannot be described adequately unless an infinitely fine grid is used. Therefore, we focus here on the relation between the fineness of the grid (i.e. the number of model parameters) and the accuracy of the inverted rigidity solution.

The coarse grid (3 × 3) does not really allow for the detection of the rigid body location (Figs 5a and 6). In particular, the rigidity

<figure><img src="image_4.png" /><figcaption>Figure 3. Sketch of optimization algorithm applied to plate rigidity inversion. This is the local version of the used search method. The context of global optimization is described in Section 1.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Benchmark geometry, relative rigidity distribution (D), boundary conditions and target velocity field $v^*$. The inner body is 10 times more rigid than the surrounding material, but the obtained velocity field is scale-invariant with respect to the rigidity distribution. Horizontal spatial units are in metres. The colour scale has been chosen to highlight the spatial variation of the relative rigidity in this objective distribution as well as in the inverted solutions shown in Fig. 5.</figcaption></figure>

contrast is largely underestimated. The medium grid (8 × 8) correctly estimates the location and the diameter of the rigid inclusion (Fig. 5b). The rigidity ratio is underestimated by a factor of 2. The

fine grid resolution (15 × 15) permits the capture of the contour of the hard material (Fig. 5c). The residual velocity $v - v^*$ associated with the three experiments reveal interesting features. First, the coarse grid has already a relatively small maximum residual (1.4 mm yr$^{-1}$) with respect to the velocity applied on the lateral sides (20 mm yr$^{-1}$). This is likely due to the fact that the difference between the target solution $v^*$ and the solution associated with a constant rigidity is already modest. Indeed, all these solutions are constrained by the same boundary conditions that limit the observed differences. However, it remains that the residual velocity dramatically decreases when the grid resolution increases (Table 1), ($J_{min}$)$^{1/2}$ for the fine grid being as low as 0.12 mm yr$^{-1}$ (Fig. 7c).

Interestingly, the location of zones of high residue is correlated with the areas where the rigidity is not well retrieved. Therefore, the spatial distribution of $v - v^*$ (that can always be mapped even if $D^*$ is unknown) provides information about the spatial distribution of the error made on the rigidity map.

In conclusion of this benchmarking, the benefit of having a fine grid in the zone of discontinuities is clearly visible as the inverted solution is greatly improved with respect to the true rigidity distribution. However, refining the parameter grid may appear to be useless and costly in regions where the rigidity remains constant. A multiscale approach for the grid definition would probably be a good strategy to get an optimal model parametrization.

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Results from (a) relative rigidity estimate D using a coarse grid (3 × 3), (b) same with a medium grid (8 × 8) and (c) same using a fine grid (15 × 15). Parameter locations are given by the circles. The colour scale is the same as in Fig. 4 where the objective rigidity distribution is shown.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Relative rigidity inverted for the three parameter grids (Fig. 5), along a profile corresponding to $Y = 50$ m.</figcaption></figure>

<table><caption>**Table 1.** Benchmark cases description.</caption><thead><tr><td>Case</td><td>Number of parameters</td><td>(J<sub>min</sub>)<sup>1/2</sup> (mm yr<sup>-1</sup>)</td></tr></thead><tbody><tr><td>1a</td><td>3×3</td><td>0.82</td></tr><tr><td>1b</td><td>8×8</td><td>0.26</td></tr><tr><td>1c</td><td>15×15</td><td>0.12</td></tr></tbody></table>
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Residual velocities associated with (a) low resolution, (b) medium resolution and (c) high-resolution inversion grid.</figcaption></figure>

## 4 RIGIDITY INVERSION IN SOUTHERN CALIFORNIA

The southern part of the San Andreas Fault System (SAFS) appears appropriate to study the relationship between the mechanical rigidity and the interseismic velocity field. Indeed, this zone is sampled by a large amount of high-quality GPS measurements published in a common reference frame (Shen et al. 2003). Moreover, the tectonics of the SAFS mostly involves horizontal motion that may be well described with a plane stress formulation.

The velocity data set used in our study is mostly based on the SCEC Crustal Motion Map version 3.0 (Shen et al. 2003), referred hereafter as CMM3 model as it was published by (Kreemer & Hammond 2007). The transient motion induced by the Landers–Hector Mine event is supposed to have been removed (Kreemer, personal communication) and we assume that these velocity data only represent interseismic motion associated with the regular loading between the Pacific and the North American plates. The $1\sigma$ uncertainties associated with these data are quite homogeneous and display an average of 1 mm yr$^{-1}$. For this reason, we do not include these uncertainties in the cost function computation and we set the matrix C to identity (eq. 11). The measurements are distributed heterogeneously (Fig. 8). Indeed, the density measurement is generally high close to active faults where high strain gradient are expected, but GPS sites are few and far between in weakly deformed areas like the Sierra Nevada or the Mojave Desert.

### 4.1 Mechanical problem and boundary conditions

As for the synthetic case, we use the forward model described in Subsection 2.3 based on the hypotheses given in Subsection 1.3 to model the strain associated with interseismic loading. Our chief assumption is that the observed velocity field is associated with the strain of a variable rigidity plate. A second assumption is that the choice of the velocity reference frame should not change the optimal parameter estimation. Indeed, if one believes that the rigidity inverted from the velocity field is an intrinsic plate property, the chosen reference frame should not affect the minimization process and its final result. The reader will find in Appendix the corresponding demonstration for the reference frame independence and its domain of validity.

To solve the mechanical problem using a plane stress FE approximation, we project the geographically defined GPS velocity field on a Cartesian frame. We use a Lambert conformal conic projection with a map projection centre set to 242° longitude and 34° latitude and the two standard parallels set to 32° and 36° latitude. We limit our study domain to a rectangular area, which includes the easternmost part of the Pacific Plate along the Californian coast to the

west, the Basin and Range and the Mojave Desert to the east and the central San Andreas Fault (SAF) segment and the south Sierra Nevada to the north. It is worthwhile to note that to the north of this selected area, the SAF loading cannot be assimilated to an interseismic loading since the fault is creeping up to the San Francisco Bay. To the south, the model ends near the Mexican border as the CMM3 data set is restricted to USA. The selected area contains 620 velocity vectors from the total data set.

Because GPS data do not directly sample the bounds of the rectangular area that we define (Fig. 8), we need to define the boundary conditions along this contour. To do so, we consider GPS data located within a band having a user-defined width (Fig. 8). We approximate separately each velocity component on the profile using a local linear regression of the velocity field in the strip. We set the bandwidth to 20 km in zones densely sampled with GPS measurements. However, the bandwidth is as large as 100 km in areas where sites are few and far between.

<figure><img src="image_9.png" /><figcaption>**Figure 8.** Spatial distribution of the geodetic measurements on the SAFS used for the inversion. The velocity field is shown (grey and red arrows) in the North American reference frame. The selected area in the rectangle (dashed lines) corresponds to our domain of analysis. The velocity distribution along the boundary conditions is shown with blue arrows. These boundary conditions are computed using a subset of the velocity field marked with red arrows (see text). Active faults are drawn in green (Jennings 1992). Letters L, H and N indicate the location of the 1993 Landers, the 1994 Northridge and the 1999 Hector Mine earthquakes, respectively. SC refers to Santa Catalina Island.</figcaption></figure>
<table><caption>**Table 2.** Summary of southern California rigidity inversions.</caption><thead><tr><th>Case</th><th>Number of parameters</th><th>(J<sub>min</sub>)<sup>1/2</sup><br>(mm yr<sup>-1</sup>)</th><th>Remarks</th></tr></thead><tbody><tr><td>2a</td><td>–</td><td>4.87</td><td>Constant rigidity (no optimization)</td></tr><tr><td>2b</td><td>6*7</td><td>2.60</td><td>Coarse grid</td></tr><tr><td>2c</td><td>12*15</td><td>2.13</td><td>Fine grid</td></tr></tbody></table>

## 4.2 Inversion results

Using previously described geometry, interpolated boundary conditions and GPS data, we perform three experiments summarized in Table 2. Cases 2b and 2c are stopped after 100 iterations allowing the cost function to reach a minimum.

### 4.2.1 Constant rigidity experiment

It is worthwhile to recall (see Subsection 2.3) that only the rigidity variation matters for the strain prediction. This would not be the case if stress would be used as a control data. Therefore, the predicted velocity field remains unchanged for any constant rigidity. We run an experiment using the boundary conditions previously described with a constant rigidity (case 2a) to evaluate the cost function without rigidity optimization. In such a case, the forward model leads to a velocity field $v$ that is markedly different from the data velocity field $v*$ with differences exceeding 6 mm yr$^{-1}$ in some places. The residual velocity field $v - v*$ is close to zero near the domain border for which the model velocity field is fully enforced by the boundary conditions that are derived from the data. The square root of the cost function is 4.87 mm yr$^{-1}$ for this experiment that is significantly higher than the $1\sigma$ uncertainty of ~1 mm yr$^{-1}$ associated with the data.

The largest residual velocities (>6 mm yr$^{-1}$) are located (1) between the coast and the SAF in southern California with an S-E orientation, (2) in central California between the SAF and the Sierra Nevada and (3) in the Mojave Desert with an N-NW orientation. This discrepancy is due to the strong distortion of the velocity field across the SAF zone (Fig. 8) that cannot be reproduced by a constant rigidity model.

<figure><img src="image_10.png" /><figcaption>**Figure 9.** (a) Relative rigidity for a coarse parameter grid of case 2b and (b) associated residual velocity field between GPS and modelled velocity.</figcaption></figure>

### 4.2.2 A coarse resolution experiment

The first optimization experiment (case 3b) is run with a 6*7 parameters rectangular grid orientated along the rectangular domain's borders. We limit the variation of the rigidity to a 1:20 ratio as we did for the synthetic benchmarks. Estimated rigidity values cover the whole parameter range [1:20] (Fig. 9a). Lowest values [1–5] are associated with the SAF zone but also with the Mojave Desert. Intermediate values occur in the eastern California shear zone ECSZ. High rigidities [10–20] appear in the Great Basin-Sierra Nevada, the Mojave Block and along the Pacific Plate. The optimization process leads the square root of the cost function to decrease down to 2.60 mm yr$^{-1}$ that is markedly less than the constant rigidity result, but still higher than the data uncertainty. Most of the residual velocities are below 5 mm yr$^{-1}$, the highest residues being concentrated in highly deformed areas around the SAF and in the Mojave Desert (Fig. 9b).

### 4.2.3 High-resolution experiment

The case 2c corresponds to a refined grid of parameters (12 × 15) still aligned with the borders orientations. The rigidity pattern provided by this experiment (Fig. 10a) is broadly consistent with the one obtained with the coarse inversion. However, the transition between low and high-rigidity areas is here sharper, suggesting that the coarse grid does not contain the high spatial frequencies required to correctly predict the observed velocity field. All zones associated with high strain and active tectonics (SAF zone, Mojave Block and eastern California shear zone) are associated with a low rigidity [1–5]. The Sierra Nevada, the Basin and Range and the Pacific Plate display high [10–20] rigidity. Surprisingly, a low-rigidity area emerges near the Santa Catalina Island in a zone where the plate could be considered as rigid. This pattern, that is not visible in the coarse inversion, may pose the problem of the relationship between the data set (spatial density and precision), the choice of the parameter grid and the quality of the inversion. This high-resolution grid leads the square root of the cost function to reach 2.13 mm yr$^{-1}$ (Table 2) that corresponds to 18 per cent of improvement with respect to case 2b. Residual velocities are not homogeneously
<figure><img src="image_11.png" /><figcaption>**Figure 10.** (a) Distribution of relative rigidity for a fine parameter grid (case 3c) and (b) associated residual velocity field between GPS and predicted velocity.</figcaption></figure>

distributed: they are relatively high (>2 mm yr⁻¹) in some strained zones along the SAF (Fig. 10b) and in the southwestern Sierra Nevada. They remain modest in most of the other areas such as the eastern hinge of the Sierra Nevada and the eastern California shear zone, the Basin and Range and in southernmost California.

### 4.2.4 *Experiments comparison along profiles*

Due to the major impact of the SAF zone on the interseismic strain in California, the relationship between the velocity field and the rigidity is well visible on a cross-section perpendicular to the fault. We show here how different inversions match the fault-parallel velocity field on two profiles (Figs 8 and 11). The north profile intersects the SAF just north of the Big Bend and also crosses the southern Sierra Nevada and the southern California Shear Zone. The south profile crosses the SAF north of the Salton Sea and ends in the south Basin and Range.

For a constant rigidity (case 2a), the lithospheric plate is homogeneous and the lateral variations of computed velocity field are mostly driven by the boundary conditions. The modelled velocity field along the two profiles is quite smooth and does not follow the high velocity gradients visible in the data. Therefore, most of the data are not fitted by the model within the errors bar.

The experiment with a coarse grid (case 2b) improves the agreement between the model and the data for the two profiles. However, the model is closer to the data for the south profile for which most of the data points agree with the model within the error bars (Fig. 11a). This better fit is probably associated with a smoother velocity gradient in this area with respect to the northern profile. Both profiles indicate a minimum of rigidity in the zone where the velocity gradient is maximum.

The high-resolution experiment (case 2c) provides a fair model-data agreement for both profiles. Most of the velocity values on the GPS points are explained within the error bars. Two notable exceptions are (1) the northeast part of the southern profile where the predicted velocities are ~2 mm lower than the observed values (see also Fig. 10b) and a group of points located 20–30 km northeast to the SAF trace on the northern profile. The high-resolution

experiment more adequately adjusts the gradients shown by the data than the coarse resolution case. The two minima of rigidity revealed by the inversion match well the location of the SAF to the west and the location of the ECSZ to the east (Figs 11e and f).

## 5 DISCUSSION

The general optimization process and the forward plane stress model described in Section 3 have been used to invert the rigidity field for two kinds of mechanical problems. The first type of benchmark includes a sharp rigidity discontinuity that cannot be described by the continuous function associated with the parametrization. In such a case, we show that the target rigidity $D^*$ associated with the velocity field $v^*$ can be approached by minimizing the cost function. A way to better describe the rigidity field $D$ is to expand the number of control points (i.e. the number of rigidity parameters to estimate), therefore generating higher spatial frequencies of the rigidity map near the imposed discontinuity. A second group of experiments aims to invert the effective rigidity of southern California using the horizontal GPS velocity field as the target $v^*$. In this case, the real effective rigidity is unknown and the boundary conditions on the model sides are not precisely known.

We analyse first the conceptual and practical difficulties associated with our methodology. Then, we compare the result of our inversion in southern California with the geophysical data that may help to constrain the lithosphere rigidity in this area.

### 5.1 Critical analysis of the inversion method

#### 5.1.1 *Definition of the cost function*

For a given distribution of rigidity, the most intuitive approach for estimating the value of the functional cost consists in evaluating the misfits right on the locations of the geodetic measurements. However, these data sets are often irregularly distributed in space. Some areas are poorly sampled whereas others probably exhibit redundant measurements. This spatial variability has to be taken into account when assessing the reliability of our optimal models of rigidity
<figure><img src="image_12.png" /><figcaption>**Figure 11.** Velocity (mm yr⁻¹) and relative rigidity associated with north profile (right-hand column) and south profile (left-hand column). (a) and (b) fault-parallel velocity model (black lines), data and errors in red; (c) and (d) fault-perpendicular velocity model; (e) and (f) normalized rigidity along the profile for coarse and fine parameter grids. Dashed lines indicate the intersection of the profile with the fault systems (EL = Elsinore Fault; SJ = San Jacinto Fault; SAF = San Andreas Fault, see Fig. 8 for profile location).</figcaption></figure>

distribution. Consequently, it sounds reasonable to integrate these geodetic measurements into a spatially continuous representation (Kreemer et al. 2003; Tape et al. 2009). This latter objective can be performed via the use of different basis functions (e.g. polynomial, exponential, splines and wavelets). This way, it would be possible to generate pseudo-data on a regular data grid especially in areas that may be lacking real measurements.

### 5.1.2 Choice of the boundary conditions

Boundary conditions (i.e. the distribution of velocities along the border of the study domain) impose heavy constraints on the minimization process. In this first approach, we define them simply by interpolating the existing velocity field nearby. However, it often happens that data may lack along the model boundary. This fact is illustrated in southern California both to the east, where very few measurements exist in the Mojave Desert, and to the west, where we attempt to extend our domain of analysis offshore. Therefore, the boundary conditions can result from a rough interpolation that may impact the rigidity estimate. To overcome this limitation and to relax this postulate on a *priori* boundary conditions, a coherent optimization approach should deal with rigidity distribution and boundary conditions all together.

### 5.1.3 Rigidity parametrization

Synthetic cases presented in this paper suggest that rigidity tends to very large values in areas that exhibit little internal deformation. Thus, its amplitude ranges from a given minimum to infinity if some zones are not significantly deformed. This semi-open varia-

tion domain is not adapted to any numerical search. Consequently, a wiser rigidity parametrization could use the compliance of the material (the rigidity inverse) instead of the rigidity. The search interval for the compliance would range from a very small value (a quasi-rigid body) to a finite value (maximum compliance). This use of compliance could insure a greater stability of the inversion process.

Another improvement of the parametrization could be to define a parameter grid adapted to the variable data density and to the spatial variation of the rigidity (or the compliance). In the previous sections, we tested various spatially uniform grids of control points by changing the sampling step. Although these tests led to similar results, it seems obvious that the spatial distribution of these control points has to roughly match the map of the local gradient of deformation. The higher is this latter, the denser should be the distribution of control points and conversely. The presence of a few points in a place where high spatial frequencies affect the velocity field may not be sufficient. Consequently, an enhancement of our approach consists in estimating the best locations of our control points, leading to an irregular parameter grid that optimally matches the velocity field. This can be performed within a hierarchical multilevel analysis framework. The use of the residual velocity map to define the optimal locations of the control points could be a way to define a comprehensive algorithm.

### 5.1.4 Uniqueness and error of the solution

The synthetic case and the natural case that we study reveal that the resulting rigidity map strongly depends on the grid choice. Also, in the case of the SAFS, the error and the spatial distribution of the
data must have an impact on the uncertainty of the estimated rigidity parameters. In this paper, we just compare the residual velocity map with the data along two profiles and we show a progressive convergence of the model velocity towards the data GPS velocities within their uncertainties. However, it remains to understand the relationship between the distribution of the velocity residues and the accuracy of parameters estimation. This major concern is beyond the scope of the paper and will be addressed in a future work.

Another important point is that the best solution $D*$ may not perfectly fit the target field $v*$ even with an optimal parameter grid design. Indeed, there is no guarantee that $v*$ corresponds to a plane stress solution even if the data are errorless. First, the 2-D plane stress model could be too rough to describe some 3-D effect occurring in the lithosphere. Also, other forces may apply to the lithosphere in addition to the lateral constraints, such as the body forces due to the topography, the basal drag from the mantle or any kind of transient motion. For example, it is notable that two spots of high velocity residual on Fig. 10b are located in the areas of 1992 Landers and 1994 Northridge earthquakes. These two large events certainly affect the local velocity field that probably exhibits some transient effects that are not taken into account in our approach. In this respect, our method could allow for the detection of a departure of the velocity field from a plane stress solution. This interesting property has been recently used to interpret an anomalous velocity field in Morocco within a convergent area (Perouse et al. 2010).

## 5.2 Geophysical interpretation of plate rigidity variations in southern California

### 5.2.1 Comparison with other rigidity estimations

The inversion of effective rigidity associated with southern California velocity field suggests that an accurate rigidity model can be computed from a dense, high-precision GPS velocity field (CMM3 model). Stable features emerge from our experiments that can be compared with other rigidity proxies.

Variations in effective flexural thickness in the United States have been estimated at large scale (Bechtel et al. 1990). Using coherence properties of topographic and gravimetric properties in the frequency domain, this study displays high flexural thickness in central and north United States then decreasing towards small values of 4–8 km in southwestern United States close to our zone of study. However, the authors note that fault bounded blocks like the Sierra Nevada may have a greater rigidity than the surrounding areas. Another study at a smaller scale in western United States (Lowry et al. 2000) seems to support this conjecture with displaying higher flexural thickness in the Sierra Nevada (~15 km) than in the Basin and Range (5 km). Direct plate gravity modelling of individual tectonic provinces such as the Sierra Nevada has been used to infer local rigidities (Kennelly & Chase 1989). These authors conclude that the Sierra Nevada has a high rigidity due to its large isostatic residual anomaly. In summary, in agreement with our results, inversions of flexural thickness based on gravity and topography provide some clues that some rigidity contrast may exist between the Sierra Nevada–Great Valley unit and the surrounding area. However, it must be recalled that this kind of estimation method assumes a constant rigidity over the window of estimation even if the true rigidity varies significantly over short spatial scales (Lowry et al. 2000). Therefore, the comparison between these estimates and the inversion provided here cannot be used to validate our approach.

Another estimation of the effective elastic thickness based on geodetic velocity field in southern California has been computed using simple 1-D models (Chery 2008). Despite the rough approximations of the governing equation of elastic plate made in this work, it emerges common features with the inversion presented here. A major contrast is found between a low rigidity of the SAF zone (Carrizo area and Salton Sea area) and a high rigidity of the Pacific Plate, the Sierra Nevada and the Mojave Block. As for the 2-D inversion, the rigidity contrast between low and high rigidity reaches 10 and more, illustrating that highly rigid zones such as the unit Great Valley–Sierra Nevada are not presently strained (Dixon et al. 2000). However, we must recognize that this comparison mostly indicates that the 1-D analysis is a *posteriori* validated by the 2-D analysis.

### 5.2.2 Relative rigidity and geophysical setting of southern California

Another way to infer the realness of our relative rigidity map $D$ is to infer the relative rigidity variations from related geophysical features. If we consider the mechanical lithosphere as a stress guide for which the inner part is elastic, the integrated elastic shear rigidity $D$ is proportional to the elastic plate thickness $H$ times its intrinsic elastic properties (see eq. 6). The former term is therefore mostly dependent on the seismological properties of the lithosphere like Young's modulus $E$. Seismological data indicate that $E$ varies laterally and with depth. If we consider a vertical profile in the lithosphere, $D$ is therefore proportional to $\bar{E}H$ where $\bar{E}$ is the average Young's modulus over the thickness $H$. Horizontal variations of the shear modulus (which is proportional to Young's modulus) have been invoked to explain asymmetric patterns of geodetic velocity fields (Le Pichon et al. 2005). However, it has been noted that lateral shear modulus contrast does not exceed a ratio of 2–3, the smallest shear modulus values generally corresponding to the fault zones. At scales larger than 10 km, the relative $V_p$ and $V_s$ variations within the crust at a given depth do not exceed 10 per cent except for shallow layers associated with sedimentary basins (Lin et al. 2007). It is therefore unlikely that intrinsic elastic properties are responsible for large lateral variations of effective rigidity. Rather, we think that mostly the effective elastic thickness $H$ varies over the continents. For the specific case of a lithosphere submitted to a shear load, it has been proposed that high heat flow areas correspond to smaller elastic thickness than the one associated with low heat flow provinces (Chery 2008). If this is true, an inverse correlation should emerge between inverted rigidity and the heat flow data of this area (Sass et al. 1994; Sass et al. 1997).

Therefore, we compare our map of rigidity distribution with the heat flow of southern and central California (Fig. 12). Being aware of the large variability of heat flow data due to perturbations of lithospheric heat flow by hydrologic or volcanic heat transfer, some trends seem to emerge however. In central California, low heat flow values (<60 mW m⁻²) are found in the Sierra Nevada–Great Valley and correlate well with high rigidities. High values (>80 mW m⁻²) occur to the east in the eastern California shear zone and to the west around the central SAF (~70 mW m⁻²). They are places where strain is localized and consequently, places where we predict low effective rigidity. Nevertheless, such a correlation vanishes in the south Basin and Range, which is characterized by a heat flow higher than 80 mW m⁻² but displays a high rigidity.

Overall, the correlation between heat flow and rigidity maps is rather, but not entirely, good. This finding has been observed at
<figure><img src="image_13.png" /><figcaption>**Figure 12.** Heat flow map produced from the Western United States Geothermal Database (http://smu.edu/geothermal/heatflow/heatflow.htm) (Blackwell & Richards 2004). Only data whose quality is high or medium have been selected. White lines are the contour lines of our relative rigidity model presented in Fig. 10(a).</figcaption></figure>

world scale (Audet & Bürgmann 2011). Indeed, these authors found a weak correlation coefficient (−0.32) between flexural thickness estimates and heat flow. It is worthwhile to recall that different factors prevent the prediction of lithosphere rigidity from heat flow measurements. In addition to shallow heat flow perturbations, large variations of radiogenic heat occur in the crust, therefore modifying the temperature distribution at depth in the lower crust and in the upper mantle (Sass et al. 1994). Also, rheological variations in the crust and the mantle due to chemical variations are likely to change the temperature of the seismic/aseismic transition (Williams 1996).

To get rid of these crustal effects, we compared our rigidity map with the seismic velocity models at the base of the lithosphere (in the upper mantle at depth ~100 km) obtained by Rayleigh wave tomography (Yang & Forsyth 2006), ambient noise and earthquake tomography (Yang et al. 2008) and surface wave tomography (Pollitz & Snoke 2010). Throughout the crust, the Sierra-Nevada exhibits high-velocity anomalies. At greater depths (~100 km), a positive anomaly stands in southern Great Valley, which has been explained as a detachment of cold lithosphere by mantle delamination. On the contrary, a low-velocity anomaly appears in the eastern Sierra-Nevada and along Walker Lane interpreted as upwelling of the asthenosphere. The same pattern exists beneath Salton Trough probably associated with the extension that affects this area. As stated at the scale of the whole western United States, these models are globally correlated with the heat flow map shown in Fig. 12. Therefore, we reach the same conclusions, namely: (1) in the Sierra Nevada-Great Valley, our estimation of high rigidity correlates with high seismic velocities and low heat flow, (2) low-velocity anomalies are found along the eastern California shear zone and the Salton Trough in agreement with our estimation of low rigidity, whereas (3) the southern Basin and Range does not exhibit either high-velocity anomalies, or low heat flow, as one could have expected from our high-rigidity estimate. This latter persistent discrepancy may be explained by the lack of dense geodetic measurements in this area located close to the limits of our study domain. An inversion of GPS velocities at a larger scale (whole western United States) will perhaps lead to a slightly different relative rigidity distribution in this area.

In addition to the heat flow versus rigidity correlation, most of low-rigidity areas are associated with active faults (Fig. 12). Strikingly, high-rigidity zones are not crossed by active faults. Although the rigidity map that we present here is only based on the knowledge of the interseismic geodetic strain with assumed locked faults, a causal link may explain the correlation. Because actively faulted zones display a high heat flow of ~80 mW m⁻², the depth of 350 °C isotherm which is usually associated with the maximum seismicity depth in the continental crust is estimated to 10–15 km (Sibson 1982). If the strength of the lithosphere is limited to the seismogenic crust in these high heat flow areas, the effective elastic thickness should be close here to 10 km. Therefore, the stress accumulation of these weak areas eventually leads the fault to failure in line with long-term mechanical modelling of active faulting in California (Bird & Kong 1994; Chéry et al. 2001). We must point out that our explanation does not apply to eastern Mojave area where high rigidity and the lack of faulting are associated with a high heat flow.

To summarize, the comparison between our inverted rigidity map of the SAF zone and the geophysical knowledge of the area suggests that a large variation of integrated rigidity has its root in horizontal thermal state variations within the lithosphere as proposed for example for northwestern America (Hyndman et al. 2009). In southern California, the location of low-rigidity areas is consistent with observed high heat flow, active faulting and the geodynamical setting of the SAF zone (Teyssier & Tikoff 1998).

Our interpretation is at odds with the mechanical understanding of interseismic strain published in the literature. Indeed, interseismic strain in fault zones is often modelled using thick elastic blocks separated by slipping faults beneath a locking depth (e.g. Savage & Burford 1973; Meade & Hager 2005) or by viscoelastic cycle deformation (Savage & Prescott 1978). Despite the fact that block models can explain data equally well, we point out that using a fault driven at depth by a dislocation is virtually equivalent to assuming a null strength on the slipping fault. If one adopts this point of view, it is then possible to interpret the locking depth of the block model as the base of a laterally variable elastic plate, the locked fault being a low-rigidity area (see Chery (2008) for a discussion). However, additional interseismic strain sources like steady slip at depth and viscoelasticity have been modelled together with effective rigidity variations (Pollitz et al. 2010). These authors conclude that much (but not all) the strain rate must be attributed to effective rigidity variation, a result that ultimately supports the approach taken in this paper.

# 6 CONCLUSION

This work is a first attempt to determine the effective rigidity of the continental lithosphere from 2-D interseismic velocity field. Both synthetic and natural cases suggest that a global optimization scheme aiming to minimize a cost function converges towards a unique solution. However, our experiments also demonstrate that the choice of the parameter space is central to minimize efficiently the cost function and to capture the strain gradient associated with velocity data. We acknowledge that some aspects of the method need improvement to strengthen our methodology. We aim to address the following aspects in a near future:

(1) The formal calculation of the resolution of the method, that is the *a posteriori* error on the parameters as a function of the data distribution and its *a priori* error;

(2) The development of a multiscale grid for the parameter space;
(3) The optimization of the boundary conditions together with the elastic rigidity. This way, the rigidity map would only depend on the choice of the geographical limits of the studied domain and the chosen parametrization;

(4) The formulation of the problem with a plane stress approach defined on the sphere to avoid a data projection onto a Cartesian frame. Ultimately, a full 3-D forward model would account for the whole equilibrium of the mechanical problem including the shape of the lithosphere and its associated elastic variations.

Beyond the inversion of effective elastic properties of the lithosphere, this new methodology could be also useful to estimate the stress rate associated with the observed geodetic strain using a consistent framework. Indeed, the optimal rigidity solution has direct strain rate and stress rate counterparts throughout the model. Therefore, the projection of the stress rate tensor over the active faults should give a prediction of stress accumulation over the considered interseismic period.

## ACKNOWLEDGMENTS

We thank Riad Hassani for providing his plane stress Finite Element code CAMEF and for a thorough reading of an early version of the manuscript. The paper was improved by the constructive comments of two anonymous reviewers and the Associate Editor.

## REFERENCES

Audet, P. & Bürgmann, R., 2011. Dominant role of tectonic inheritance in supercontinent cycles, *Nat. Geosci.*, **4**, 184–187, doi:10.1038/NGEO1080.

Bechtel, T.D., Forsyth, D.W., Sharpton, V.L. & Grieve, R.A.F., 1990. Variations in effective elastic thickness of the North American lithosphere, *Nature*, **343**, 636–638.

Bills, B.G., Currey, D.R. & Marshall, G.A., 1994. Viscosity estimates for the crust and upper mantle from patterns of lacustrine shoreline deformation in the Eastern Great Basin, *J. geophys. Res.*, **99**, 22 029–22 086.

Bird, P. & Kong, X., 1994. Computer simulations of California tectonics confirm very low strength of major faults, *Bull. geol. Soc. Am.*, **106**, 159–174.

Blackwell, D.D. & Richards, M., 2004. Calibration of the AAPG geothermal survey of North America BHT data base, *AAPG Annual Meeting*, Dallas, TX, Poster session, paper 87616.

Chery, J., 2008. Geodetic strain across the San Andreas Fault reflects elastic plate thickness variations (rather than fault slip rate), *Earth planet. Sci. Lett.*, **269**, doi:10.1016/j.epsl.2008.01.046.

Chéry, J., Zoback, M.D. & Hassani, R., 2001. An integrated mechanical model of the San Andreas Fault in central and northern California, *J. geophys. Res.*, **106**(B10), 22 051–22 066.

Chinnery, M.A., 1963. The stress changes that accompany strike-slip faulting, *Bull. seism. Soc. Am.*, **5**, 921–932.

D'Alessio, M.A., Johanson, I.A. & Burgmann, R., 2005. Slicing up the San Francisco Bay area: block kinematics and fault slip rates from GPS-derived surface velocities, *J. geophys. Res.*, **110**(B06403), doi:10.1029/2004JB003496.

Dixon, T., Miller, M., Farina, F., Wang, H. & Johnson, D., 2000. Present-day motion of the Sierra Nevada block and some tectonic implications for the Basin and Range province, North American cordillera, *Tectonics*, **19**(1), 1–24.

Ergintav, S., Dogan, U., Gerstenecker, C., Cakmak, R., Belgen, A., Demirel, H., Aydin, C. & Reilinger, R., 2007. A snapshot (2003–2005) of the 3D postseismic deformation for the 1999 Mw = 7.4 Izmit earthquake in the Marmara Region, Turkey, by first results of joint gravity and GPS monitoring, *J. Geodyn.*, **44**, 1–18, doi:10.1016/j.jog.2006.12.005.

Gerbault, M., 2000. At what stress level is the central Indian Ocean lithosphere buckling? *Earth planet. Sci. Lett.*, **178**(3–4), 165–181.

Hetenyi, H., Cattin, R., Vergne, J. & Nabelek, J.L., 2006. The effective elastic thickness of the India Plate from receiver function imaging, gravity anomalies and thermo-mechanical modelling, *Geophys. J. Int.*, **167**, 1106–1118.

Holland, J.H., 1975. *Adaptation in Natural and Artificial Systems*, University of Michigan Press, Ann Arbor, MI.

Hyndman, R., Currie, C.A., Mazzotti, S. & Frederiksen, A., 2009. Temperature control of continental lithosphere elastic thickness, Te vs Vs, *Earth planet. Sci. Lett.*, **277**, 539–548.

Ivorra, B., 2006. Semi-deterministic global optimization and industrial applications, *PhD thesis*, Université de Montpellier II.

Ivorra, B., Universite 2006. Semi-deterministic vs. genetic algorithms for global optimization of multichannel optical filters, *Int. J. Comput. Sci. Eng.*, **2**(3), 170–178.

Jennings, C.W., 1992. *Preliminary Fault Activity Map of California*, California Division of Mines and Geology, Sacramento, Open-File Report 92-03.

Kaufmann, G. & Amelung, F., 1995. Reservoir-induced deformation and continental rheology in vicinity of Lake Mead, Nevada, *J. geophys. Res.*, **105**, 16 341–16 358.

Kennelly, P.J. & Chase, C.G., 1989. Flexure and isostatic residual gravity of the Sierra Nevada, *J. geophys. Res.*, **94**, 1759–1764.

Kirkpatrick, S., Gelatt, Jr, C.D. & Vecchi, M.P., 1983. Optimization by simulated annealing, *Science*, **220**(4598), 671–680.

Kreemer, C. & Hammond, W.C., 2007. Geodetic constraints on areal changes in the Pacific–North America plate boundary zone: what controls basin and range extension? *Geology*, **35**(10), 943–946.

Kreemer, C., Holt, W. & Haines, J., 2003. An integrated global model of present-day plate motions and plate boundary deformation, *Geophys. J. Int.*, **154**, 8–34.

Le Pichon, X., Kreemer, C. & Chamot-Rooke, N., 2005. Asymmetry in elastic properties and the evolution of large continental strike-slip faults, *J. geophys. Res.*, **110**(B03405), doi:10.1029/2004/B003343.

Lin, G., Shearer, P.M., Hauksson, E. & Thurber, C.H., 2007. A three-dimensional crustal seismic velocity model for southern California from a composite event method, *J. geophys. Res.*, **112**(B11306), doi:10.1029/2007JB004977.

Lowry, A.R., Ribe, N. & Smith, R.B., 2000. Dynamic elevation of the Cordillera, western United States, *J. geophys. Res.*, **105**, 23 371–23 390.

McNutt, M., 1990. Flexure reveals great depth, *Nature*, **343**, 596–597.

Meade, B.J. & Hager, B.H., 2005. Block models of crustal motion in southern California constrained by GPS measurements, *J. geophys. Res.*, **110**(B03403), doi:10.1029/2004JB003209.

Milne, G.A., Davis, J.L., Mitrovica, J.X., Scherneck, H.-G., Johansson, J.M., Vermeer, M. & Koivula, H., 2001. Space-geodetic constraints on glacial isostatic adjustment in Fennoscandia, *Science*, **291**, 2381–2385.

Mohammadi, B. & Pironneau, O., 2009. *Applied Shape Optimization for Fluids*, Oxford University Press, Oxford.

Mohammadi, B. & Saiac, J.-H., 2003. *Pratique De La Simulation Numérique*, Dunod, Paris.

Nur, A. & Mavko, G., 1974. Postseismic viscoelastic rebound, *Science*, **183**, 204–206.

Perouse, E., Vernant, P., Chery, J. & Reilinger, R., 2010. Active surface deformation and sublithospheric processes in the western Mediterranean constrained by numerical models, *Geology*, **38**, 823–826.

Pollitz, F.F. & Snoke, J.A., 2010. Raleigh-wave phase-velocity maps and three-dimensional shear velocity structure of the western US from local non-plane surface wave tomography, *Geophys. J. Int.*, **180**, 1153–1169.

Pollitz, F.F., McCrory, P., Wilson, D., Svarc, J., Puskas, C. & Smith, R.B., 2010. Viscoelastic-cycle model of interseismic deformation in the north-western United States, *Geophys. J. Int.*, **181**, 665–696.

Reilinger, R. et al., 2006. GPS constraints on continental deformation in the Africa-Arabia-Eurasia continental collision zone and implications for the dynamics of plate interactions, *J. geophys. Res.*, **111**, doi:10.1029/2005JB004051.

Sambridge, M., 1999. Geophysical inversion with a neighbourhood algorithm – I. Searching a parameter space, *Geophys. J. Int.*, **138**, 479–494.
Sass, J.H., Lachenbruch, A.H., Galanis, S.P., Jr., Morgan, P., Priest, S.S., Moses, T.H., Jr. & Munroe, R.J., 1994. Thermal regime of the southern Basin and Range Province: 1. Heat flow data from Arizona and the Mojave desert of California and Nevada, *J. geophys. Res.*, **99**, 22 093–22 119.

Sass, J.H., Williams, C.F., Lachenbruch, A.H., Galanis, S.P., Jr. & Grubb, F.V., 1997. Thermal regime of the San Andreas fault near Parkfield, California, *J. geophys. Res.*, **102**, 27 575–27 585.

Savage, J.C. & Burford, R.O., 1973. Geodetic determination of relative plate motion in central California, *J. geophys. Res.*, **78**, 832–845.

Savage, J.C. & Prescott, W.H., 1978. Asthenosphere readjustment and the earthquake cycle, *J. geophys. Res.*, **83**, 3369–3376.

Shen, Z.K. et al., 2003. The SCEC Crustal Motion Map, Version 3.0, Available at: http://epicenter.usc.edu/cmm3 (last accessed 2011 September).

Sibson, R.H., 1982. Fault zone model, heat flow, and the depth distribution of earthquakes in the continental crust of the United States, *Bull. seism. Soc. Am.*, **72**, 151–163.

Tape, C., Muse, P., Simons, M., Dong, D. & Webb, F., 2009. Multiscale estimation of GPS velocity fields, *Geophys. J. Int.*, **179**(2), 945–971.

Tarantola, A., 1987. *Inverse Problem Theory, Methods for Data Fitting and Model Parameter Estimation*, Elsevier, Amsterdam.

Teyssier, C. & Tikoff, B., 1998. Strike-slip partitioned transpression of the San Andreas fault system: a lithospheric-scale approach, in *Continental Transpressional and Transtensional Tectonics*, Vol. 135, pp. 143–158, eds Holdsworth, R.E., Strachan, R.A. & Dewey, J.F., Geological Society, London.

Thatcher, W. & Pollitz, F.F., 2008. Temporal evolution of continental lithospheric strength in actively deforming regions, *GSA Today* **18**(4–5), doi:10.1130/GSAT01804-5A.1.

Ugray, Z., Lasdon, L., Plummer, J., Glover, F., Kelly, J. & Marti, R., 2007. Scatter search and local NLP solvers: a multistart framework for global optimization, *INFORMS J. Comp.* **19**(3), 328–340.

Watts, A.B., 2001. *Isostasy and Flexure of the Lithosphere*, Cambridge University Press, Cambridge.

Weissel, J.K. & Karner, G.D., 1998. Flexural uplift of rift flanks due to mechanical unloading of the lithosphere during extension, *J. geophys. Res.*, **94**, 13 919–13 950.

Williams, C.F., 1996. Temperature and the seismic/aseismic transition: observations from the 1992 Landers earthquake, *Geophys. Res. Lett.*, **23**, 2029–2032.

Yang, Y. & Forsyth, D.W., 2006. Rayleigh wave phase velocities, small-scale convection, and azimuthal anisotropy, beneath southern California, *J. geophys. Res.*, **111**, B07306, doi:10.1029/2005JB004180.

Yang, Y., Ritzwoller, M.H., Lin, F.C., Moschetti, M.P. & Shapiro, N.M., 2008. Structure of the crust and uppermost mantle beneath the western United States revealed by ambient noise and earthquake tomography, *J. geophys. Res.*, **113**, B12310, doi:10.1029/2008JB005833.

Zienkiewicz, O.C., Taylor, R.L. & Zhu, J.Z., 2005. *The Finite Element Method: Its Basis and Fundamentals*, Elsevier Butterworth-Heinemann.

## APPENDIX: VELOCITY REFERENCE FRAME INVARIANCE

To demonstrate this assumption in the context of the minimization of the cost function, one needs to prove that the gradient of *J* does

not change if the velocity fields *v* and *v*\* are replaced, respectively, by rotated and translated fields *v*<sub>R</sub> and *v*<sub>R</sub>\*. Let us start with a simple calculus for the gradient of a functional based on least-square differences of two velocity fields as given in eq. (11). In its discrete formulation, and considering a covariance matrix equal to the Identity matrix, the gradient of this functional for a rotated and translated velocity field in a Cartesian coordinate system is given by

$$ 
\begin{align}
\nabla J_R &= \operatorname{Grad}[(v_R - v_R^*)^T W^T W (v_R - v_R^*)] \nonumber \\
&= (\nabla v_R^T W^T W (v_R - v_R^*) + (v_R - v_R^*)^T W^T W \nabla v_R) \nonumber \\
&= (\nabla v^T R^T W^T W R (v - v^*) + (v - v^*)^T R^T W^T W R \nabla v) \nonumber \\
&= (\nabla v^T W^T W (v - v^*) + (v - v^*)^T W^T W \nabla v) \nonumber \\
&= \nabla J \tag{A1}
\end{align}
$$

with

$$ v_R = R(v + T), $$

$$ v_R^* = R(v^* + T), $$

$$ R^T R = I, $$

where *R* denotes a diagonal matrix of identical 2-D-rotation matrices, *T* a translation vector and *W* the weighting diagonal matrix, uniform for the east- and north-velocity component of one single sample. This proves analytically that the gradient of functional based on least-square differences of two velocity fields with respect to the rigidity parameters is the same whatever rotated and translated velocity fields being used.

If we include the covariance matrix of the data *C* in the cost function, *J* can be rewritten in a discrete formulation as eq. (11).

$$ J = (v - v^*)^T W^T C^{-1} W (v - v^*). \tag{A2} $$

In such a case, the above analysis holds if $R^T W^T C^{-1} W R = W^T C^{-1} W$. Indeed, the gradient in its matrix formulation is given by

$$ 
\begin{align}
\nabla J_R &= \operatorname{Grad}[(v_R - v_R^*)^T W^T C^{-1} W (v_R - v_R^*)] \nonumber \\
&= (\nabla v_R^T W^T C^{-1} W (v_R - v_R^*) + (v_R - v_R^*)^T W^T C^{-1} W \nabla v_R) \nonumber \\
&= (\nabla v^T R^T W^T C^{-1} W R (v - v^*) + (v - v^*)^T R^T W^T C^{-1} W R \nabla v) \nonumber \\
&\neq (\nabla v^T W^T C^{-1} W (v - v^*) + (v - v^*)^T W^T C^{-1} W \nabla v). \tag{A3}
\end{align}
$$

The equality of the last line of eq. (A3) would require that no correlation exists both between the uncertainties of the east and north component for each site, and between the uncertainties associated with two different sites. Thus, the inverse matrix of covariance is not diagonal. Consequently, the change of velocity reference frame is likely to induce some changes in the gradient of our functional to minimize. Hence, we have no guarantee to converge to the same optimal model even though our tests suggest that this effect is likely to be negligible.