Joint estimate of the rupture area and slip distribution of the 2009 L ’ Aquila earthquake by a Bayesian inversion of GPS data

S U M M A R Y Usually, when inverting geodetic data to estimate the slip distributions on a fault, the area is made large enough to more than cover the rupture zone, with regularization producing regions of large slip with very small slip over the rest of the surface. We have developed a new inverse method which assumes that nonzero slip is confined to a rectangular region, and which jointly estimates, using Bayesian methods, the boundaries of this region as well as the slip distribution within it, using a smoothing parameter also determined as part of the inversion. Synthetic tests show that our method can successfully image deeper slip regions not resolved by previous methods, and does not produce spurious regions of nonzero slip. We apply our method to coseismic displacements measured by GPS for the 2009 L’Aquila earthquake, first determining the orientation of the fault assuming a simplified model with uniform slip, and then determining probability density functions for the location, length, and width of the rupture area and for the slip distribution. The standard deviation of slip is about 10 cm and describes a normal-faulting earthquake with a maximum slip of 88 ± 11 cm and seismic moment of 3.32+0.30 −0.29 × 1018 N m.


I N T RO D U C T I O N
Surface displacements and gravity anomalies observed by means of different geodetic techniques (GPS, DInSar, and GRACE and GOCE satellite missions) have been used to study earthquakes.The information that we can obtain by inversion of these observations depends on the earthquake and Earth models and the model parameters that we decide to estimate.In addition to oversimplified earthquake models, like point-like seismic sources for estimating the focal mechanism, magnitude and centroid location (Cambiotti & Sabadini 2013) or rectangular finite faults with uniform slip providing also a first indication of the dimensions of the rupture area (Anzidei et al. 2009;Devoti et al. 2012;Cambiotti & Sabadini 2012;Wang et al. 2012), the scientific community has devoted particular attention to heterogeneous slip distributions over the fault surface in order to obtain a detailed description of the rupture and to determine the stress drop within the gouge of the fault (Serpelloni et al. 2012;Wang et al. 2015), as well as to improve the match of the modelling with observations.
Our understanding of earthquakes indicates that they can be described as a rupture of a specific area of a new or already existing fault (i.e. the area unlocked by the earthquake) and that this rupture mainly consists of a slip across the fault (or tangential dislocation).According to arguments from continuum mechanics, the slip distribution must be continuous over the fault (Dahlen & Tromp 1998).
At present, in order to make the inverse problem overdetermined and to avoid unrealistic slip distributions, geodetic data inversions rely on prior constraints on the roughness of the slip distribution.There exist several objective or subjective methods to quantify the degree of smoothing from information from observations (Fukuda & Johnson 2008).Nevertheless, none of these methods has considered the rupture area where the slip distribution is confined as an unknown which must be determined as part of the inversion.Indeed, the extension of the finite fault is always fixed a priori and chosen in a subjective way on the basis of where the slip is expected to be localized and the rupture area is inferred a posteriori by inspection of the spatial distribution of the estimated slip.With the aim of minimizing the impact of this choice on the results (Zhou et al. 2014), the extension of the finite fault is often chosen large enough to cover the expected region of largest slip and part of the surrounding region where the slip goes smoothly to zero.
In this work we will present a new inverse method that incorporates prior information on the slip roughness and on the requirement that the slip must take place in a specific rupture area.Within a Bayesian framework, where the full knowledge about the unknowns can be described by probability distributions, we will objectively estimate the extension of the rupture area, jointly with the slip distribution and its degree of smoothing, from GPS data.In contrast to the classical approach, this method does not need to make a specific choice about the extension of the finite fault and, for this purpose, we will rely on the fully Bayesian approach proposed by Fukuda & Johnson (2008) starting from the Akaike Bayesian information criterion (Akaike 1980;Yabuki & Matsu'ura 1992).
As a case study, we will apply this inverse method to the 2009 L'Aquila earthquake.This earthquake has been widely studied with the aid of geodetic (GPS and DInSAR) and strong motion data in a two-step approach where, first, the fault geometry is determined under the assumption of uniform slip and, then, used to estimate the slip distribution (Walters et al. 2009;Atzori et al. 2009;Cheloni et al. 2010;Serpelloni et al. 2012).In the first step, this approach assumes that the slip is uniform over a rectangular rupture area and geodetic data are inverted for estimating two parameters describing the uniform slip (its along strike and updip components) and seven parameters describing the extension of the rectangular rupture area and the planar geometry of the fault (the length and width of the rectangle, the depth of its top edge, the dip and strike angles and the geographic coordinates of a chosen point of the rectangle, e.g. its centre).These parameters are then used in the second step to define the planar fault over which the slip distribution is estimated from the same data set used in the first step.Nevertheless, the extension of the finite fault for this second inversion is subjectively chosen and taken larger than that estimated in the first step so that it could contain the region of largest slip and part of the surrounding region where the slip goes smoothly to zero.It is worth mentioning that this procedure to determine the fault plane is typical for geodetic data inversions.In waveform inversions, instead, the fault plane is constrained from the nucleation point position, centroid parameters (location, strike, dip), including sometimes aftershock distributions, and varied until it delivers the best waveform.
Our inverse method can be regarded as an improvement of the two step approach where, once the fault geometry has been established in the first step under the assumption of uniform slip, the slip distribution is estimated from observations, jointly with the rupture area where the slip takes place.In this respect, our method makes objective also the estimate of how large and where the rupture area is located and, in this way, avoids making the a priori choice of a sufficiently large finite fault over which the slip distribution is inverted.For the sake of simplicity, we will keep the assumption of planar fault surfaces and rectangular rupture areas and we will adopt flat homogeneous elastic Earth models in order to simulate surface displacement caused by earthquakes (Okada 1992), rather than relying on complex, spherically symmetric Earth models as in Cambiotti et al. (2011) and Cambiotti & Sabadini (2015).

E A RT H Q UA K E M O D E L
Let us consider a flat Earth model and assume that the fault surface is planar and that we know its dip angle and line of strike (i.e. the intersection between the planar fault and the Earth surface).In order to identify the points of this planar fault, as depicted in Fig. 1, we introduce the along strike, x 1 , and updip, x 2 , coordinates with respect to an origin taken on the line of strike.
We assume also that the rupture area reaches the Earth surface and is rectangular.In this respect, it can be specified by its along strike, L, and along dip, W, lengths and by the along strike coordinate, x 1 = l, of the mid-point.We put these three parameters in the following 3-D array where A is the space of all possible rupture areas reaching the Earth surface.It is possible to consider also rupture areas which do not reach the Earth surface by including a fourth parameter indicating the updip coordinate of the top edge of the rupture area.For the sake of simplicity and in view of the application to the 2009 L'Aquila earthquake, we do not consider this possibility here.Following Yabuki & Matsu'ura (1992), we parametrize the slip distribution using bicubic splines defined over a regular grid of the rectangular rupture area, with N 1 nodes along strike and N 2 nodes along dip.The bicubic splines of those nodes adjacent to and on the edges of the rupture area are defined in such a way that, at the left, right and bottom edges, the slip is zero and that, at the top edge, the second derivative of the slip with respect to the updip coordinate is 994 G. Cambiotti et al. zero.The latter constraint allows the slip to reach the Earth surface and, at the same time, avoids anomalous results at this edge.We provide details about the regular grid and the bicubic splines in Supporting Information Appendix A.
We decompose the slip distribution, u, into its along strike, u 1 , and updip, u 2 , components (2) where x1 and x2 are the along strike and updip unit vectors, and x is the vector position over the fault surface Then, we express each slip component as follows: where n is the bicubic spline of the nth node, u k,n is the slip coefficient associated to this node and to the along strike (k = 1) or updip (k = 2) component, and A is the rupture area in the fault reference system We note that, by definition of the modified bicubic splines, the slip distribution goes to zero at the left, right and bottom edges of the rupture area and, by definition of rupture area, it is set to zero outside.Furthermore, we note that the updip coordinates x 2 = −W and 0 identify the bottom and the top edge of the rupture area.According to eq. ( 4), the slip is parametrized by a linear combination of expansion (basis) functions, which are represented by overlapping bicubic splines, defined over N = N 1 N 2 nodes.We collect the slip coefficients u k,n into the following (2 N)-dimensional array where 2 N is the total number of coefficients, S is the space of all possible coefficients and s k is given by The number of nodes of the regular grid is established in such a way that the grid steps do not exceed a spatial resolution specified both along strike, δl 0 , and along dip, δw 0 , that is, with the cautionary remark that the along strike, N 1 , and updip, N 2 , numbers of nodes are never chosen smaller than 4 in order to define properly the (modified) bicubic splines.

I N F O R M AT I O N F RO M O B S E RVAT I O N S
Within the framework of the dislocation theory and linear elastic Earth models, coseismic surface displacements depend linearly on the slip coefficients and non-linearly on the rupture area (Okada 1992;Yabuki & Matsu'ura 1992).The relation between data and model parameters thus takes the following form ỹ = y(a, s) + e (10) where ỹ and e are the M-dimensional arrays of the observed surface displacements and errors, and y is the M-dimensional array of modelled surface displacements y(a, s) = K(a) s.
Here, M is the number of data and K is the M × (2N)-matrix that linearly relates the modelled surface displacements, y, to the slip coefficients, s, and depends nonlinearly on the rupture area, a.The elements of this matrix, K iq , can be computed from the elastostatic Green functions for tangential dislocations (e.g.Okada 1992) Here, the index i identifies the component of the surface displacement at a specific GPS site (i.e. a datum), the index q is given by q = (k − 1) N + j according to the arrangement of slip coefficients introduced in eqs ( 6) and ( 7), j is the modified bicubic spline centred at the jth node of the regular grid, and G k i is the elastostatic Green function for the along strike (k = 1) and updip (k = 2) point-like unit slip.
The array of errors, e, contains observational errors from data analysis and modelling errors that result, for example, from nonplanar fault geometries and heterogeneities in the Earth structure, which are not taken into account by our earthquake and Earth models.Due to difficulties in assessing modelling errors and following Fukuda & Johnson (2008), we use a simple representation of observational and modelling errors and assume that they obey a Gaussian distribution with zero mean and covariance matrix where E is the covariance matrix of observational errors and σ ∈ (0, ∞) is a hyperparameter that must be estimated as part of the inversion.The latter can be interpreted as follows: there are modelling errors when σ > 1, they are negligible when σ = 1 and the earthquake model is over-parametrized (i.e. it fits the observations too well) when σ < 1.Alternatively, we can say that observational errors have been underestimated during the data analysis when σ > 1, correctly estimated when σ = 1 and overestimated when σ < 1.We note that, by definition, this assumption on the structure of observational and modelling errors guarantees the statistical agreement between observed and modelled surface displacements, on average and within one-sigma error.The statistical agreement is obtained by rescaling the observational errors with the hyperparameter σ and, in the most common case of σ > 1, this means that there are modelling errors or that observational errors have been underestimated during the data analysis.
From eqs (10), ( 11) and ( 13), the probability density function (PDF) for the data given the earthquake model parameters and the hyperparameter σ , can be written as where X y is the misfit between observed and modelled surface displacements Bayesian inversion of GPS data 995 and • and T stand for the determinant and the transpose, respectively.Given the data ỹ, eq. ( 14) depends on the earthquake model parameters (s and a) and of the hyperparameter σ and, thus, provides the information on the parameters, supplied by the observations (Yabuki & Matsu'ura 1992).

I N F O R M AT I O N F RO M P R I O R C O N S T R A I N T S
It is common to introduce as prior information that the slip distribution is smooth to some degree, to make the inverse problem overdetermined and avoid implausible results.In this respect, we consider the measure of slip roughness proposed by Yabuki & Matsu'ura (1992) and defined as follows: This involves the integration over the rupture area, A, of the second order derivative of the along strike, u 1 , and updip, u 2 , components of the slip with respect to the along strike, x 1 , and updip, x 2 , coordinates.
Although this kind of constraint on the slip roughness is arbitrary (alternative definitions of slip roughness have been considered (Zhou et al. 2014)), its implementation is mandatory to invert static earthquake data recorded by different geodetic techniques.Only the inclusion of kinematic earthquake data, together with the physics of the rupture process describing plausible time evolutions of the slip distribution, allows to remove this kind of constraint (Minson et al. 2013).
Substituting eq. ( 4) into eq.( 16), we obtain or, equivalently, Here B is a symmetric positive definite N × N-matrix, whose (n, n )th element can be obtained analytically performing the following integral and A is the following block (2 N) × (2N)-matrix We note that the matrices B and A depend on the rectangular rupture area, a.They also depend on the parametrization of the slip distribution and, particularly, on the along strike and along dip spatial resolutions used to define the regular grid, eq. ( 9).With this measure of slip roughness, the prior PDF for the slip coefficients reads where β ∈ (0, ∞) is a hyperparameter that weights information from the constraint on the slip roughness; we made use of the relationship A = B 2 .Eq. ( 21) can be seen as the prior PDF for the slip coefficients given the rupture area, a, and the hyperparameter β.Within this framework we are not including any positivity constraints on the slip distribution, meaning that we are not introducing any prior information about the slip direction.Although it would be possible to include this kind of prior information, it is not necessary for the case of the 2009 L'Aquila earthquake because of the objective estimate of the smoothing parameter β, as discussed in Sections 6 and 7.
In addition to this prior PDF, we may introduce also the prior PDF for the rupture area and the hyperparameters, say f (a, σ, β).Combining this PDF with the PDF for the slip given the rupture area and the hyperparameter β, we thus obtain the prior PDF for all earthquake model parameters and hyperparameters As far as we have no knowledge about these parameters, but for the domain in which they take their values, that is,

P O S T E R I O R P D F
According to the Bayes' theorem (Tarantola 2005;Fukuda & Johnson 2008), the PDF for the data and the prior PDF, eqs ( 14) and ( 23), can be combined in order to obtain the posterior PDF for all the earthquake model parameters and hyperparameters given the data where c 0 is a constant of normalization.Eq. ( 24) takes into account information from observations and prior constraints, and constitutes the solution of the inverse problem.It can be used to calculate any estimator, as optimal parameters, means, standard deviations and confidence intervals, as well as to investigate the marginal PDFs for selected earthquake model parameters or hyperparameters.From now on, for the sake of brevity, we will not explicit any longer the dependence on the data.We will thus write χ y (a, s) and p(a, s, σ, β) instead of χ y (a, s, ỹ), eq. ( 15), and p(a, s, σ, e | ỹ), eq. ( 24), respectively.
By considering that eq. ( 24) depends on s only via the argument of the exponential function and that the latter is a bilinear form in s, this equation consists of a linear inverse problem for the slip coefficients given the rupture area and the hyperparameters.Following Yabuki & Matsu'ura (1992), the conditional PDF for the slip given the rupture area and the hyperparameters can be cast in the form of a Gaussian distribution with mean s and covariance matrix σ S where On the basis of this result, it is straightforward to marginalize the posterior PDF with respect to the slip s.We obtain where p(a, σ, β) is the marginal PDF for the rupture area and the hyperparameters

The fully Bayesian approach
As discussed in Fukuda & Johnson (2008), the approach proposed by Yabuki & Matsu'ura (1992) is based on the Akaike Bayesian information criterion (ABIC) and does not account for the uncertainties of the hyperparameters, which is appropriate only if the marginal PDF for the hyperparameters has a dominant and sharp peak at its maximum.For this reason and with the purpose of including positivity constraints on the slip distribution, Fukuda & Johnson (2008) proposed a fully Bayesian approach which mainly consists in estimating the mean slip distribution and its standard deviation by integration of the posterior PDF over the slip coefficients and hyperparameters in order to account also for the uncertainties of the hyperparameters.The fully Bayesian approach, once applied to the present case without positivity constraints, allows us to account for the uncertainties in the rupture area, in addition to those of the hyperparameters, and therefore we will rely on it.
In this perspective and following Fukuda & Johnson (2008), we marginalize p(a, σ, β) with respect to the hyperparameters and obtain where p(a) is the marginal PDF for the rupture area and p(σ, β | a) is the conditional PDF for the hyperparameters given the rupture area Furthermore, we define the expectation, E(•), and covariance, C( • , •), operators as follows where h and g are two functions of the earthquake model parameters and hyperparameters, and E(•| a) is the expectation operator given the rupture area We will use the expectation and covariance operators to estimate the earthquake model parameters, the hyperparameters and their uncertainties.As discussed in Supporting Information Appendix B, and for some specific functions of the earthquake model parameters, the integration over the slip coefficients can be performed analytically thanks to the fact that the conditional PDF for the slip takes the form of a Gaussian distribution, eq. ( 25).Furthermore, by substituting the hyperparameter β with the new hyperparameter γ given by eq. ( 28), it is possible to perform analytically the integration over the hyperparameter σ .On the other hand, the integrations over the rupture area a and the new hyperparameter γ must be performed numerically and we do it using the Romberg algorithm (Quarteroni et al. 2000).
In particular, as described in detail in Supporting Information Appendix B, the mean and covariance of the along strike, u 1 , and updip, u 2 , components of the slip at a given point x of the planar fault are given by ) and require numerical integration over only the new hyperparameter γ .For the sake of simplicity, when we will discuss the mean slip amplitude and its variance, we will consider the following approximation.Given the definition of the slip amplitude, the approximation of its expectation and variance read For the case study of the 2009 L'Aquila earthquake, we consider the surface displacements obtained from the time-series analysis of 50 continuous survey and 27 campaign GPS stations by Serpelloni et al. (2012), Fig. 2 (black arrows).Considering that each GPS station records the three (east, north and vertical) components of the surface displacement, but for the MSNI GPS station for which there is no the vertical component, the number of data is M = 77 × 3 − 1 = 230.In this case, the covariance matrix of the observational errors, E, is diagonal because the GPS time series analysis has been obtained assuming that there is no noise correlation among GPS stations and the different components of the recorded surface displacements.As already obtained by Serpelloni et al. (2012), the fault surface estimated from the inversion of these GPS data under the assumption of uniform slip is characterized by dip and strike angles of 50.5 • and 129.4 • , respectively, and the mid-point of the top edge of the extended finite fault (i.e. a point of the line of strike) is at (13.4714 • E, 42.3609 • N).These values of the dip and strike angles and of the mid-point of the top edge of the fault are kept fixed and used as exact prior information in the following GPS data inversion.We choose the mid-point of the top edge of the extended finite fault area as origin of the fault reference system introduced in Section 2.
For computational reasons, we take a finite rupture area space where the along strike coordinate of the mid-point, l, varies from −5 to 5 km, and the length, L, and the width, W, vary from 10 to 40 (40) We sample this space over a regular grid of 65 nodes for each parameters (i.e. 2 n + 1 with n = 6 according to the discretization required by the Romberg algorithm), which totals 65 3 = 274 625 nodes.Furthermore, as in Serpelloni et al. (2012), we fix at 2 km the spatial resolution of the regular grid used for the parametrization of the slip distribution, both along strike and along dip.

Rupture area
Fig. 3 shows the marginal PDFs for the along strike coordinate of the mid-point, p(l), the length, p(L) and the width, p(W), of the rectangular rupture area and their comparison with Gaussian distributions with the same means and standard deviations of these marginal PDFs.The estimates of the means, standard deviations and 95 per cent confidence intervals are reported in each panels.We note that the marginal PDFs are bell shaped and have only one maximum at l = 0.47 km, L = 18.44 km. and W = 18.13 km.In this respect, we note that the restriction to a finite rupture area space, eq. ( 40), does not affect our conclusion because the probability of those rupture areas that we are excluding is negligible.Furthermore, the Gaussian distribution is a good approximation for l and L, while it is not for W due to the asymmetry of this parameter at larger values.This asymmetry is due to the fact that the resolving power of surface coseismic displacements on the slip distribution decreases with increasing depth and to the poor coverage of GPS network at the Earth surface above the deep portion of the fault (from Fig. 2, we note that there is no coverage southwestwards the ROIO GPS station for about 25 km, up to the PSCA GPS station; see also the discussion in Supporting Information Appendix C).These results show that the GPS data for the 2009 L'Aquila earthquakes resolve the rupture area with a precision (95 per cent confidence intervals) of a few kilometres along strike (l = 0.61 +1.50  −1.38 km and L = 18.93 +3.59 −2.74 km), and several kilometres along dip (W = 20.72 +14.22  −4.02 km).Particularly, the midpoint of the top edge of the rupture area is at (13.4771 • E, 42.3574 • ).We note that our estimates of the length and width are greater than those obtained by Serpelloni et al. (2012) under the assumption of uniform slip, which are L = 14.3 +3.8  −2.4 km and W = 15.8 +7.9 −5.5 km: this is due to our definition of length and width referring to the whole region of the fault unlocked by the earthquake, even with small slip, which is necessarily larger than that obtained under the assumption of uniform slip, which is controlled by high slip values.

Slip distribution
Fig. 4 shows the mean slip distribution and its standard deviation in the fault and geographic reference frames, which we obtain using eqs (37) and (39).We note that the slip distribution describes a normal earthquake with average rake angle of −103.3 •+4.2 −4.8 and is well localized around the region of largest slip.Furthermore, it is characterized by a single (global) maximum of 88.2 ± 10.8 cm at depth of 6.2 km.The seismic moment is M S = 3.32 +0.30  −0.29 × 10 18 N m and the moment magnitude is M W = 6.281 +0.025 −0.026 assuming a shear modulus of 26 GPa as in Serpelloni et al. (2012) (derived from Di Luzio et al. 2009).Left-lateral and right-lateral strike slips of about 20 cm characterize the shallowest few kilometres of the fault northwestwards and southeastwards the PAGA GPS station, respectively.The standard deviation is of the order of 10 cm in the mean rupture area, with a maximum of 17 cm at 12.5 km depth.Outside the rupture area it decreases to zero in less than a few kilometres along strike and in a few tens of kilometres downdip.This difference in the dependence of the standard deviation with respect to the along strike and updip coordinates shows how, according to the fully Bayesian method, uncertainties on the dimensions of the rupture area (especially those on the width Fig. 3) are translated to uncertainties on the slip distribution.Furthermore, we note that the standard deviation increases in the shallowest portion of the fault.This increase can be explained considering that the shallower the slip is, the more localized is the deformation at the Earth surface and, so, the more dense must be the GPS network to resolve it.Only in the very proximity of GPS stations, the shallowest slip can be resolved, as shown by the local minimum of the standard deviation attributable to the nearby PAGA GPS station, of about 3 cm.Deep and intermediate slip, instead, deforms larger regions of the Earth surface and, thus, can be resolved by more distant GPS stations.
The slip and, in particular, the standard deviation that we have estimated according to the fully Bayesian approach are not exactly zero outside the mean rupture area because they have been obtained by means of the expectation and covariance operators, eqs ( 35) and ( 36), which take into account the slip distribution for all possible rupture areas weighted by the probability of the latter, p(a).This implies that we estimate (almost) zero slip and standard deviation only at those points of the fault surface which do not belong to likelihood rupture areas, that is, sufficiently far from the mean rupture area.
The one asperity model for the synthetic test T1 discussed in Supporting Information Appendix C, Fig. S1(A), closely resembles the mean slip distribution estimated from the original data, Fig. 4(a), and it is well recovered by the GPS network.More complex slip distributions, like those characterized by two and four asperities for the synthetic tests T2 and T3, Supporting Information Fig. S1(B,C), can also be recovered by the GPS network, although the ability of resolving multiple slip maxima deteriorates with increasing depth.In light of this, we asses that GPS data for the 2009 L'Aquila earthquakes do not indicate the presence of one shallow and one deep asperities as instead suggested by broad band ground motion seismometers (Cirella et al. 2012;Gallovič et al. 2015).Since Ameri et al. (2012) and Tinti et al. (2014) have shown clear evidence of two asperities, separated in time, in the seismic data, it is likely that this difference reflects aseismic afterslip that affects the GPS offsets over their 24 hr sampling interval.

Observed and modelled surface displacements
According to the fully Bayesian approach, it is not meaningful to consider modelled surface displacements for specific (e.g.mean or optimal) rupture areas and slip distributions.Rather, similarly to the expectation of the slip distribution, we must consider the expectation of modelled surface displacements, that is, with E(s | a) given by eq.(B.19a).Furthermore, according to the assumption made on the structure of observational and modelling errors, we must estimate also the hyperparameter σ in order to rescale the covariance matrix from the GPS data analysis, eq. ( 13).
Particularly, from eqs (B.10) and (B.11), we obtain E(σ ) = 9.39 and C(σ , σ ) = 0.85 and, then, the estimate (mean and standard deviation) of this hyperparameter, σ = 9.39 ± 0.92.In light of this, the standard deviations of the east, north and vertical components of the observed surface displacements must be increased by a factor of √ σ = 3.06 ± 0.96.Fig. 2 shows the comparison between observed, ŷ, and modelled, E(y), surface displacements and their residuals.The root mean square (RMS) errors are 4.3 4.5 10.1 mm for the east, north and vertical components of the surface displacements, respectively, and, as implied by the fully Bayesian method, observed and modelled displacements are in overall agreement within one-sigma error.This is shown in Figs 2(b,d), where we compare the residual between observed and modelled surface displacements with the observational errors rescaled by the square root of the hyperparameter σ .Particularly, we note that, in the near field, the agreement between observed and modelled surface displacements is better in the hanging wall than in the foot wall, both for the horizontal and vertical components (the PAGA GPS station is just in the hanging wall).Furthermore, the agreement is better in the far field than in the near field, although this is mainly due to the fact that far field surface displacements are smaller than their uncertainties.
The weighted residual square sum (WRSS), defined as follows: yields WRSS/M = 8.5, with M = 230 being the number of data.This definition of WRSS constitutes an indication of the quality of the fit between observed and modelled surface displacements with respect to the original observational errors from GPS data analysis, that is, without scaling the data covariance matrix by the hyperparameter σ .We note that Serpelloni et al. (2012) obtained a value of WRSS (7.1 M) smaller than that here obtained (8.5 M) and this is mainly attributable to the different (subjective and objective) methods used to quantify the slip roughness.Indeed, as we will discuss in Section 7, the subjective method implemented by these authors selects a slip distribution that is rougher than ours, thus allowing a better fit to the observations.

THE A P R I O R I C H O I C E O F T H E E X T E N S I O N O F T H E F I N I T E FAU LT
In order to gain insight into the effects of the a priori choice of the extension of the finite fault on the estimate of the slip distribution, we now invert the GPS data for different choices of the extension of the finite fault.This means that, instead of eq. ( 24), we will consider the posterior PDF given by the conditional PDF for the slip and the hyperparameters given the rupture area and we will make different choices of the rupture area parameters.
The comparison between the results obtained in this way and those discussed in Section 6, will show the advantages of our method and enlighten the consequences of fixing the extension of the finite fault.This comparison is also discussed in Supporting Information Appendix C on the basis of synthetic tests.Figs 5(a) and (b) show the slip distribution estimated assuming the same finite fault as chosen by Serpelloni et al. (2012) and characterized by l = 0 km, L = 30 km and W = 24 km.We note that the slip distribution resembles that shown in Fig. 4(a), except for the fact that it is slightly smoother.Particularly, it goes to zero more slowly towards the borders and the maximum slip is 85.3 ± 8.5 cm.Furthermore, we note that this choice of the finite fault yields smaller standard deviations in the region of largest slip and in the surrounding few kilometres than those obtained according to our method (Fig. 4), and larger standard deviations elsewhere.It is worth mentioning also that this a priori choice of the extension of the finite fault sets the standard deviation arbitrarily to zero at depths greater than 24 km, where our method predicts instead small, but nonzero, standard deviations.By definition, indeed, the standard deviation, as well as the slip, is zero outside the chosen finite fault according to the classical method.
Figs 5(c) and (d) show the slip distribution estimated assuming a larger finite fault characterized by l = 0 km, L = 80 km and W = 80 km.We note that the slip distribution is smoother than that obtained for the smaller rupture area shown in Fig. 5(a) and (b), and the maximum slip is 78.1 ± 5.8 cm.This is due to the fact that the larger finite fault requires a higher degree of smoothing to stabilize the slip distribution.The largeness of the chosen finite fault also yields a few local maxima far from the region of largest slip, which can be regarded as the attempt of the slip distribution to improve the match of the modelling with observations, even at locations where unlocking is not expected on physical grounds.As specific example of this behaviour, we put our attention on the local maximum close to the line of strike and at about the along strike coordinate x 1 = 25 km, which improves the fitting to the surface displacement measured by the nearby BOMI GPS station.Similar to the case of the smaller finite fault, the standard deviations in the region of largest slip and in the surrounding few kilometres are smaller than those obtained by our method (Fig. 4), while they  43).Results are shown for the 30 × 24 km 2 and 80 × 80 km 2 rupture areas discussed in the main text (top and bottom panels, respectively).The red rectangles indicate the 30 × 24 km 2 rupture area, which is the same as chosen by Serpelloni et al. (2012).The black arrows and circles in panels (a) and (b) indicate the mean slip and standard deviation and the contour lines are given every 10 and 2 cm for the mean slip distribution and standard deviation, respectively.are larger elsewhere.The further increase of the dimension of the finite fault would make things even worse, up to the scenario where the uncertainties of the slip distribution becomes larger and larger going far away from the region of largest slip.
These results show that the pattern and amplitude of the slip distribution and, especially, of the standard deviation depend on the a priori choice of the extension of the finite fault.Our inverse method, instead, does not have this dependence and, objectively estimating the rupture area from information from observations, is able to provide the rupture area where the slip actually takes place and translates uncertainties on the dimensions of the rupture area to uncertainties on the slip distribution.
By comparing the slip distribution shown in Fig. 5(a) with the slip distribution obtained by Serpelloni et al. (2012) we can also understand the impact of objective, rather than subjective, methods for estimating the slip roughness.Indeed, both results have been obtained using the same data set and finite fault, and, thus, they differs only by the way in which the prior constraint on the slip roughness has been implemented.The fully Bayesian method estimates a slip distribution smoother than that obtained by Serpelloni et al. (2012) who choose the value of the smoothing near to the location on the trade-off curve where only a small decrease in the misfit between observed and modelled surface displacements is gained by adding more slip roughness (Harris & Segall 1987;Du et al. 1992).Our slip distribution, indeed, is characterized by only one (global) maximum of 85.2 ± 8.5 cm, while that obtained by Serpelloni et al. (2012) has a global maximum of more than about 100 cm and local maxima away from the region of largest slip.In light of the ability of our method to resolve multiple maxima (see Supporting Information Appendix C), these local maxima should be considered as artefacts rather than being an information provided by the GPS data.A similar remark holds also for the case of the slip distribution obtained by Cheloni et al. (2010) from GPS data, which is indeed characterized by four maxima and, then, rougher than ours.Finally, we note that the objective method does not need positivity constraints on the slip distribution in order to estimate realistic rake angles.
Another example of the difficulties encountered in fixing a priori the extension of the finite fault is offered by the slip distribution estimated through the single value decomposition (SVD) inversion of DInSAR data by Atzori et al. (2009) where, as the same authors say, patches with null slip or very high uncertainty at the border have been removed, without giving further details.In contrast to these previous works, the slip distribution obtained by Walters et al. (2009) inverting DinSAR and body-waves is smooth and closely resembles our results.Nevertheless, we note that this agreement is mainly due to the fact that these authors fix a priori a very small finite fault of 20 × 18.5 km 2 , which is almost the mean rupture area that we have objectively estimated from GPS data, about 19 × 21 km 2 .In this respect, we can say that their choice was right, but not inferred from observations.These results show that our methodology, building on the objective estimation of the extension of the rupture area from observations, overcomes the need of a somehow subjective choice of the extension of the finite fault pertaining to the second step of two step approach for retrieving the heterogeneous slip distribution once the planar fault geometry has been established.This advantage, at the same time, is achieved without introducing any bias on the estimate of the slip distribution.

C O N C L U S I O N S
We have presented a new inverse theory that incorporates prior information both on the slip roughness and on the evidence that slip must take place in the region unlocked by the earthquake, that is, the rupture area.This theory, based on the fully Bayesian approach proposed by Fukuda & Johnson (2008), objectively estimates the roughness of the slip distribution, as well as the location and dimensions of the rupture area, from information from observations.For the sake of simplicity, we have assumed that the fault surface is planar and that the rupture is rectangular and extends up to the Earth surface.
Our objective estimation of the rupture area from information from GPS data constitutes a step ahead of previous inverse methods aiming at estimating the heterogeneous slip distribution over a finite fault.Indeed, independently from the use of subjective or objective methods for constraining the slip roughness (Fukuda & Johnson 2008), all previous inverse methods fix a priori the exten-sion of the finite fault.This means that they implicitly assume that the slip, as well as its standard deviation, is zero outside the chosen finite fault.In contrast, our method consider the probability of each rupture area and estimate the mean slip distribution by integration, weighted according to these probabilities, over the whole space of possible rupture areas.In this context, our new inverse method does not make a unique choice of the extension of the finite fault and predicts nonzero slip and standard deviation only at those points belonging to likelihood rupture areas.
We have applied this new inverse method to the 2009 L'Aquila earthquake, by employing the same GPS data of Serpelloni et al. (2012), as well as their estimates of the dip angle and of line of strike.From the analysis of the posterior PDF, we conclude that GPS data are able to constrain the rupture area with a precision (95 per cent confidence intervals) of a few kilometres along strike and several kilometres along dip (see Fig. 3).The estimated slip distribution clearly describes a normal earthquake of magnitude M W = 6.281 +0.025 −0.026 with average dip angle of −103.3 •+4.2 −4.8 .It is characterized by only one (global) maximum of 88.2 ± 10.8 cm at depth of 6.2 km.The standard deviation is of the order of 10 cm within the mean rupture area and decreases to zero in a few kilometres along strike and a few tens of kilometres downdip where the slip goes to zero.In contrast, we have shown that the standard deviation of the slip distribution estimated fixing a priori the extension of the finite fault depend on this a priori information.
As shown by the synthetic tests discussed in Supporting Information Appendix C, our fully Bayesian method does not introduce any bias in the resolving single-and multi-asperity slip distributions.Instead, it provides results which closely resemble those obtained according to the classical approach when the a priori choice of the extension of the finite fault is optimal, that is to say that it coincides exactly with the rupture area of the synthetic slip distribution.In realistic situations, where this information is not available in advance, our fully Bayesian method performs even better than the classical approach in resolving the synthetic slip distribution, especially when there exist deep asperities.
We have also shown that the objective estimation of the slip roughness makes the slip distribution smoother than those obtained by using subjective methods.This avoids the appearance of spurious regions of nonzero slip away from the region of largest slip that, in view of the ability of our objective method to resolve multi asperities, should be considered as artefacts of subjective methods rather than being an information provided by GPS data.Furthermore, the objective estimate of the slip roughness for the case of the 2009 L'Aquila earthquake avoids the inclusion of positivity constraints on the slip distribution, which have been implemented in other works, as in Serpelloni et al. (2012).We have also shown that the subjective choice of a too large finite fault may lead to the appearance of artefacts, such as unrealistic local maxima away from the region of largest slip, apart from the subjective or objective method used to estimate the slip roughness.In this respect, our new inverse method, jointly estimating the slip distribution and the rupture area where the slip takes place, solves this issue.Overall, despite the computational effort in estimating the slip distribution for different rupture areas, our new inverse method allows to obtain more realistic slip distributions and a more reliable estimate of its uncertainties compared to previous works.

A C K N O W L E D G E M E N T S
This research was partially supported by the CAS/CAFEA international partnership program for creative research teams (No.

Figure 1 .
Figure 1.A cartoon depicting the projection on the Earth surface of the along strike and updip axes of the fault reference system and of the rectangular rupture area.The along strike and updip coordinates of the vertexes and of the mid-points of the top and bottom edges are indicated.

)
The expectations of the slip coefficients, E(u k,n | a), and of their products, E(u k,n u k ,n | a), given the rupture area are elements of the expectations E(s | a) and E(s s T | a), which are given by eq.(B.19

Figure 2 .
Figure 2.Comparison between observed and modelled (a) horizontal and (c) vertical components of the surface coseismic displacements (black and grey arrows, respectively).Panels (b) and (d) show the residual (i.e. the difference) between observed and modelled surface coseismic displacements.The thin circles represent the standard deviation of the observed horizontal and vertical displacements, after having been rescaled by the square root of the expectation of the hyperparameter σ , that is, by E(σ ) = 3.06.In this way we account for the structure of both observational and modelling errors estimated a posteriori within the framework of the fully Bayesian approach.The dashed dot line and rectangle indicate the line of strike taken fromSerpelloni et al. (2012) and the estimated mean rectangular rupture area.

Figure 3 .
Figure 3. Marginal PDFs for (a) the along strike coordinate of the mid-point, p(l), (b) the length, p(L), and (c) the width, p(W), of the rectangular rupture area and their approximation with Gaussian distributions (solid and dashed lines, respectively).The grey regions indicate the 95 per cent confidence intervals.In each panel, we also report the extremes of the 95 per cent confidence intervals, the means and standard deviations of the marginal PDFs.

Figure 4 .
Figure 4. (a,c) Mean slip distribution and (b,d) standard deviation in the fault and geographic reference frames (top and bottom panels, respectively) obtained by the posteriori PDF, eq.(24).The grey arrows and circles in panel (a) indicate the mean slip and standard deviation and are given every 2 km, and the contour lines are given every 10 and 2 cm for the mean slip distribution and standard deviation, respectively.The red rectangle indicates the mean rupture area, and the grey small circles and labels in panels c and d show the locations and names of the GPS stations.

Figure 5 .
Figure 5. (a,c) Mean slip distribution and (b,d) standard deviation in the fault reference frame obtained by the conditional PDF given the rupture area, eq.(43).Results are shown for the 30 × 24 km 2 and 80 × 80 km 2 rupture areas discussed in the main text (top and bottom panels, respectively).The red rectangles indicate the 30 × 24 km 2 rupture area, which is the same as chosen bySerpelloni et al. (2012).The black arrows and circles in panels (a) and (b) indicate the mean slip and standard deviation and the contour lines are given every 10 and 2 cm for the mean slip distribution and standard deviation, respectively.