An integral method to estimate the moment accumulation rate on the Creeping Section of the San Andreas Fault

SUMMARY Moment accumulation rate (also referred to as moment deﬁcit rate) is a fundamental quantity for evaluating seismic hazard. The conventional approach for evaluating moment accumulation rate of creeping faults is to invert for the slip distribution from geodetic measurements, although even with perfect data these slip-rate inversions are non-unique. In this study, we show that the slip-rate versus depth inversion is not needed because moment accumulation rate can be estimated directly from surface geodetic data. We propose an integral approach that uses dense geodetic observations from Interferometric Synthetic Aperture Radar (InSAR) and the Global Positioning System (GPS) to constrain the moment accumulation rate. The moment accumulation rate is related to the integral of the product of the along-strike velocity and the distance from the fault. We demonstrate our methods by studying the Creeping Section of the San Andreas fault observed by GPS and radar interferometry onboard the ERS and ALOS satellites. Along-strike variation of the moment accumulation rate is derived in order to investigate the degree of partial locking of the Creeping Section. The central Creeping Segment has a moment accumulation rate of 0.25–3.1 × 10 15 Nm yr − 1 km − 1 . The upper and lower bounds of the moment accumulation rates are derived based on the statistics of the noise. Our best-ﬁtting model indicates that the central portion of the Creeping Section is accumulating seismic moment at rates that are about 5 per cent to 23 per cent of the fully locked Carrizo segment that will eventually be released seismically. A cumulative moment budget calculation with the historical earthquake catalogue ( M > 5.5) since 1857 shows that the net moment deﬁcit at present is equivalent to a M w 6.3 earthquake.


I N T RO D U C T I O N
Recent observations of unusually large slip near the shallow toe of the 2011 Tohoku, Japan earthquake (Fujiwara et al. 2011) have shaken the conventional view of earthquake and tsunami hazards on faults that are potentially creeping. There are two hypotheses that attempt to explain the unusually large slip near the trench: first, a recent study by Noda & Lapusta (2013) suggests that even fully creeping faults may participate in a large rupture in the presence of 'dynamic weakening', which involves heating and thermal expansion of pore fluids during rapid shear. Alternatively, large slip near the trench might be due to small but significant moment accumulation on the shallowest portion of the subduction zones during the * Now at: Department of Earth and Space Sciences, University of Washington, 4000 15th Avenue NE, Seattle, WA 98195, USA.
interseismic period (i.e. partial locking on the plate interface), which is difficult to recover due to the lack of current offshore geodetic measurements. Earthquakes in the past may not have completely filled the 'seismic gap' in the shallowest portion of the plate interface, resulting in moment accumulating beyond one seismic cycle (Lorito et al. 2011). The estimation of the moment accumulation rate along the plate interface will provide insights on the fraction of the seismic and aseismic slip on the fault interface and the origin of large earthquakes along partially creeping faults.
One place to study the moment accumulation rate along creeping faults on the continents is at the Creeping Section of the San Andreas fault (SAF). It is an approximately 150 km long segment extending from Parkfield to San Juan Bautista in central California ( Fig. 1; Bakun et al. 2005). The central Creeping Section is historically void of large earthquakes (M > 7) though significant earthquakes (e.g. the 2003 M6.5 San Simeon earthquake and the 2004 M6.0 Parkfield earthquake) have occurred near this region. If the central section is completely unlocked, it could act as a barrier to a future through-going rupture. While the centre section creeps aseismically at a rate comparable to the long-term fault slip rate (Sieh & Jahns 1984;Argus & Gordon 2001;Titus et al. 2006;Moore & Rymer 2007;Tong et al. 2014), the creep rate is lower towards the ends of the Creeping Section and thus some seismic moment may be accumulating (Tong et al. 2013;Field et al. 2014;Lienkaemper et al. 2014;Maurer & Johnson 2014).
For the entire San Andreas Fault System, the moment accumulation along a given fault should be balanced by the moment release from all the past earthquakes (Reid 1910;Savage & Burford 1973;Thatcher & Rundle 1979). Therefore, the moment released during the next big earthquake can be estimated as the current moment accumulation rate times the time since the last major event. The common approach to deriving the moment accumulation rate is through constructing elastic block models constrained by geodetic measurements (Smith & Sandwell 2003;Hammond & Thatcher 2005;McCaffrey 2005;Meade & Hager 2005;d'Alessio et al. 2006;Bird 2009;Smith-Konter et al. 2011). These models are used to estimate slip rates along faults, which are multiplied by the depth of the seismogenic zone and the shear modulus to arrive at moment accumulation rate. However, a correction must be applied in the case of creeping faults. For example, surface creep rates are estimated for all fault segments using various techniques such as GPS, InSAR, and alignment arrays when developing the Uniform California Earthquake Rupture Forecast Version 3 (UCERF3) earthquake hazard model (Field et al. 2014). In the UCERF3 study, when surface creep rates are less than 50 per cent of the long-term fault slip rates, a mechanical model based upon stress equilibrium (Savage & Lisowski 1993) is used to account for shallow creep rate variations with depth. Alternatively, when the surface creep rates are 100 per cent of the slip rates, the faults are assumed to be completely unlocked so no moment accumulates. When the creep rates reach from 50 to 100 per cent of the long-term slip rates, Field et al. (2014) employed an empirical approach based on the relationship of the creep rate and the fraction of locking. Field et al. (2014) also pointed out that the moment reduction caused by shallow creep could be more accurately solved for the areas having sufficiently dense geodetic measurements.
We propose a novel approach to estimate the moment accumulation rate along strike-slip faults during the interseismic period of the seismic cycle, constrained by dense near-fault geodetic observations, such as the GPS array or radar interferometry data. We use dense geodetic data on the central portion of the Creeping Section of the SAF to directly integrate surface velocities to estimate the moment accumulation rate and therefore bypass the non-unique inversion for the depth variations in slip rate. Finally, we compare the results with those derived by the moment-bounding inverse method.

Integral method
The moment accumulation rate can be recovered by studying the surface velocity profile of a long strike-slip fault. A long strike-slip fault that is creeping from the surface to great depth will have a step-like velocity profile (Fig. 2). By taking the difference between the observed surface velocity and the step-like velocity, we obtain the residual surface velocity. If some portion of the fault is locked, there will be a residual surface velocity v r after removing the steplike velocity profile, which is denoted by the back-slip rate (Savage & Burford 1973) versus depth function s. Assuming a 2-D dislocation in an elastic half-space, the fault-parallel residual velocity as a function of distance x from the fault is v r (x): where z is the depth, s(z) is the back-slip (also referred as slip deficit) rate as a function of depth, and D m is the maximum slip deficit depth. The back-slip rate function s (z) for the simple 2-D model can be used to calculate the moment accumulation rate M per unit fault length L as where we assumed a constant shear modulus μ. The shear modulus μ is set to 30 GPa throughout our study. As shown in Fig. 2, slip rate-versus-depth functions are very different among the boxcar slip distribution, the stress-free crack model, and a tapered crack model yet they produce nearly identical velocity profiles (Segall 2010). Our ability to differentiate these slip rate-depth profile depends on the accuracy and resolution of the measurements. Indeed the inverse problem of solving for the slip rate-versus-depth profile from surface velocity measurements is the classical problem used to illustrate non-uniqueness. The inverse problem can be stabilized by making assumptions about the smoothing and the non-negativity of the model solution (Parker 1994). Nonetheless, there is still an infinite set of solutions to the proposed problem (1) that satisfy the observations. Downloaded from https://academic.oup.com/gji/article-abstract/203/1/48/577323 by guest on 15 March 2020 Rather than solving for the slip rate variations with depth, we are interested in the moment accumulation rate (sometimes also called moment deficit or moment deficit rate). It can be shown that the moment deficit rate is directly related to the integral of the surface velocity times the distance from the fault in a 2-D dislocation model (Appendix) where W is the maximum distance from the fault to perform this integral. When the finite width W goes toward infinity, the integral approaches the true moment accumulation rate asymptotically. In eq. (3) the locking depth D m disappears, which is simply a manifestation that the shape of the velocity profile reflects the depth of locking. For example, a fault with a deeper locking depth will have a greater moment accumulation rate because the integral area enclosed by its velocity profile and the step-like function is larger (Fig. 2). We are enlightened by this theoretical calculation, which demonstrates that the moment can be obtained directly from geodetic measurement v r . It is thus unnecessary to solve the slip rate-versusdepth profile in an inverse problem to calculate moment accumulation rate. To illustrate this idea, the moment accumulation rate from three different slip-rate models are shown in Fig. 2 (shaded grey). It can be proved that various slip-rate models that predict the same velocity profiles will have the identical moment accumulation rate (Appendix).
As shown in eq. (3) the moment accumulation rate is directly linked to the integral of the surface velocity field. A spatially continuous velocity field with high precision and high resolution is generally not readily available because the surface velocities are usually observed at discrete locations by ground-based instruments (e.g. GPS sites). Interpolation of the discrete GPS velocity vectors can induce unreliable artefacts due to the sensitivity of the spatial distribution of the data and the interpolation scheme. InSAR can provide a continuous velocity field but the noise of the InSAR observables, especially at long wavelengths, could cause significant bias to the estimation. This is particularly true when estimating the contributions to the moment from velocity measurements far from the fault where the moment arm is large. Small errors in these far-field velocity measurements will result in a large error in the moment rate. Our solution to this problem is to divide the moment accumulation rate into two components depending on wavelength, taking advantage of the strength of the respective GPS and InSAR techniques. The continuous velocity field from InSAR is high-pass filtered so the integrated velocity far from the fault is zero. Far-field contributions to the moment, based on the 3-D earthquake cycle model, are added to the near-field data to establish the total moment accumulation rate (Smith & Sandwell 2006;Tong et al. 2013Tong et al. , 2014.
We propose a simple method to integrate the high-pass filtered InSAR-derived interseismic velocities to estimate the shortwavelength moment accumulation rate. This high-pass filtering approach can be represented by where ν r is the filtered residual velocity data, H is the high-pass filter, and l is the low-pass filter. Specifically, we used a Gaussian 2σ 2 in this study, where σ defines the filter width. The high-pass filtering effectively reduces the long wave-length noise in the InSAR data contaminated by the delay effects of the atmosphere and/or the ionosphere.
The high-pass filtered moment accumulation rate is denoted as: The integral width W should be determined by the filter wavelength σ and the observations. We must not let W go to infinity when calculating M L in that ν r will go to zero at a certain distance away from the fault. The asymptotic property of eq. (3) is preserved in the low-pass filtered component of the moment and this longwavelength component should be combined with eq. (5) to obtain the final estimate. In Appendix, we have used analytic formulations to demonstrate the modulation effect of the filtering on the slip rate-versus-depth function.

Moment-bounding inversion
Another method for obtaining moment accumulation rate is to solve a non-negative least squares (NNLS) problem with an equality constraint to the moment (Johnson et al. 1994;Maurer & Johnson 2014). This NNLS problem can be expressed as where G is the design matrix derived from an elastic dislocation model (i.e. Okada 1985), s is unknown distribution of the back-slip rate, and v is the observable, that is residual velocities at the surface.
is the diagonal matrix containing the standard deviation of each observable. R represents the first derivative of the slip rate distribution and λ is its weight. a is a vector that represents the discretized form of the integral of s which is the moment accumulation rate functional. M is the desired moment rate. The goal here is to find a complete set of the functional that fit the observables within uncertainties. The uncertainties of the observables provide statistical significance of the measurements so a range of alternative solutions can be evaluated.
We implemented a new scheme of the moment-bounding inversion method. The inversion matrix, according to eq. (6), is written as: ⎡ where G ALOS , G ERS and G GPS are the corresponding elastic responses to unit slip s of a rectangular fault patch calculated within an elastic half-space. The dimensions of the fault patch are ∼10 km in the strike direction and 3 km, increasing to 9 km, in the down-dip direction (Section 4.3). For GPS, the elastic response contains the horizontal displacement in east and north direction. For InSAR, the elastic responses are converted into Line-Of-Sight (LOS) direction using horizontal displacement and information on the radar look geometry. v ALOS , v ERS and v GPS are the residual velocities from ALOS, ERS and GPS respectively. −1 AL O S , −1 E RS and −1 G P S are diagonal matrices representing the inverse of the standard deviations of the measurements. In order to integrate short-wavelength deformation from InSAR with GPS, we designed and applied a high-pass filter H to both the velocity field and the elastic response function for InSAR data. The parameter p is a weight to restrict the moment accumulation rate, which is set to be large enough so that the last row of eq. (7) is always satisfied for all the possible solutions of s.
To obtain the final estimate, we perturb the value of M (eqs 6 and 7) to solve a set of NNLS problems. The bounds of M range from zero to a value corresponding to a fully locked fault. The interval of the perturbation of M is usually one tenth of the upper bound of M. For along-strike variations in the 3-D problem, we need to adjust function a for each set of NNLS problems to specify which fault segment to enforce the constraint. The faults are divided into 21 fault segments (see Section 4.3) so we need to solve for 21 sets of NNLS problems. Overall, the moment-bounding inversion calls for solving hundreds of NNLS problems under the same observational constraints. Before inversion, data outliers must be carefully removed in order to satisfy the presumed probability density function (Parker & Song 2005) because the bounds of the moment accumulation rate determined by the statistics are sensitive to the noise structure of the velocity data. We found that (see Section 3 for details) the presence of data outliers would bias the chi-squared distribution of the observables (Parker & Song 2005), resulting in unrealistic estimation of the bounds.

DATA ANALYSIS
The GPS velocity data used in this study is adopted from Rolandone et al. (2008;Fig. 3). 225 horizontal GPS velocity vectors span our study area (Fig. 3a). Some of the GPS sites are within 1 km of the active creeping zone and are prone to local effects (addressed in Section 4.2). Data coverage is relatively dense near Parkfield where significant post-seismic motion has occurred since the 2004 event. GPS data outliers, defined as the data points with residuals that are greater than 1 standard deviation, were removed by constructing a preliminary best-fitting model. To reduce the observed GPS velocities to the residual velocities, we first constructed a 3-D elastic half-space model by assuming no moment accumulation along the entire extent of the creeping faults and removed this forward model from the observation. The fault model extends from the surface to 9999 km in depth and is approximately 200 km in length, closely (i.e. within 500 m) following the fault geometry derived from the quaternary fault trace mapped by the U.S. Geological Survey (Jennings 1975). The fault plane is set to be dipping vertically, consistent with the locations of the microseismicity (Waldhauser & Schaff 2008). In the half-space model, we set the long-term slip rate for the Creeping section to be 34 mm yr −1 based on a joint-inversion constrained by various geodetic and geologic data (Tong et al. 2014). To the north, the slip rate is partitioned to 17 mm yr −1 along the Creeping Section and 17 mm yr −1 along the Paicines segment of the Calaveras fault (Tong et al. 2014).
The predicted GPS velocity from a fully creeping fault is shown in Fig. 3(b) and the residual velocities are shown in Fig. 3(c). The along-strike variations of the residual velocities are readily discernable, which can be attributed to the transition from fault locking to creeping. The Parkfield section has the largest residual velocities, illuminating deformation due to partial locking of the fault. There is a significant amount of residual velocity (∼5 mm yr −1 ) at the central Creeping Section. The residual velocities to the north of the San Juan Bautista segment, where the fault returns to a locking state, are larger than the central Creeping Section but smaller  ALOS/PALSAR data are along ascending orbits and ERS data are along descending orbits. The dots with three different colours highlight the steps in the SAR image alignment process. The black lines show the high-quality interferograms that are used in deriving the InSAR velocity data. The grey lines show the low-quality interferograms that are contaminated by atmosphere noise or unwrapping errors. The 'master' images are marked by solid lines and the 'surrogate' images are marked by dashed lines  than the Parkfield region. We used this GPS velocity field in two ways: first, we treated it as an independent data set to validate the InSAR-derived result along the central Creeping Section. Second we incorporated it directly into the moment-bounding inversion to analyse the along-strike variation of the moment accumulation rate for the entire region.
The InSAR data were acquired by both the C-band radar onboard the ERS satellite from the European Space Agency (ESA) and the L-band radar onboard the ALOS satellite from the Japan Aerospace and Exploration Agency (JAXA). The two satellites acquired SAR data at two different time periods: ERS-1 and 2 spans from 1992 to 2009 and ALOS spans from 2007 to 2011. Three ascending tracks from ALOS (track 219, 220 and 221) and two descending tracks from ERS (track 27 and 256) cover the entire Creeping Section from the locked region south of Parkfield to the north of the Calaveras-Paicines fault. The baseline versus time information, along with the InSAR pair selections, for the InSAR data used in this study is provided in Fig. 4. C-band ERS data suffer from temporal decorrelation so fewer interferograms are usable compared to the L-band ALOS data. Methods related to the InSAR data processing and integration with GPS observations are discussed in Tong et al. (2013).
We used the InSAR stacking method to estimate the secular component of the interseismic velocities (Tong et al. 2013). The velocity field measured by InSAR is in the look direction of the radar and thus it is called the line-of-sight (LOS) velocity. We combined the LOS velocities from the descending ERS and ascending ALOS satellites (Wright et al. 2004) to distinguish horizontal from vertical crustal deformation.
We applied a spatial Gaussian high-pass filter with a 40 km full wavelength to the InSAR LOS velocities (Tong et al. 2013). The In-SAR velocities have spatial resolution of approximately 0.5 km. We designed a down-sampling algorithm based on the second invari-ant of the strain rate field to reduce redundant InSAR data points. We first determined seven categories of sampling rate based on the second invariant of the strain rate tensor estimated by simply taking the spatial gradient of the interpolated GPS velocity field. The highest velocity gradient region is near the San Andreas fault, so the InSAR data near the fault are sampled at a spatial resolution of 0.01 • in longitude by 0.008 • in latitude. The lowest gradient region is further away from the active fault and the data in these regions are sampled at a spatial resolution of 0.08 • in longitude by 0.064 • in latitude. The regions with moderate gradient have sampling rates that lie between these two end-member scenarios. Before each decimation step, we applied a medium low-pass filter to the data to avoid aliasing effects. After decimation, there are 8730 data points from ALOS and 4990 data points from ERS. Finally, analogous to the GPS comparison described above, a 3-D elastic half-space model, assuming no moment accumulation throughout the entire Creeping Section, was removed from the raw observed velocities to reduce these to the residual InSAR velocities. Note that the elastic model was filtered the same way as the raw InSAR velocities before taking the difference. Figs 5(a) and (d) show the filtered residual velocity from InSAR that were used in the next step of analysis. This velocity field reflects the amount of locking on the creeping faults.
There is great deal of information embedded in this velocity map. Significant residual velocities concentrate near the faults, indicating the amount of partial locking. The broadness and the magnitude of the residual velocity field reveals clearly the amount of moment accumulation rate. These velocities range from −3 to 3 mm yr −1 in radar LOS direction, which is indeed a subtle signal successfully imaged by radar on board the satellites. The moment deficit rate results based on Fig. 5 will be discussed in the next section.
We note that there are complexities in the residual velocity field, which cannot be explained by the model (Figs 5c and f). These  (Lindsey et al. 2014). We can make inferences on the direction of the deformation since we have two look directions from both the ascending and descending orbits. By comparing consistent patterns in the residual velocities fields (Figs 5c and f), we infer that there are areas of subsidence (shown as red) around the Parkfield region as well as at the junction of the San Andreas fault and the Calaveras-Paicines fault. More accurate vertical geodetic measurements are needed to investigate these issues.
We first focus on the velocity profile taken across the central Creeping Section (profile marked by a dash line in Fig. 5). This transect extends 20 km in width and 40 km in length, centred at (−120.927 • , 36.338 • ). The region covered by this transect is marked by a box in Figs 1 and 5. We projected the high-pass filtered residual data points within the box onto the profile that is perpendicular to fault strike (orientation of 50 • clockwise from north). Next we decomposed the ALOS and ERS LOS velocities into faultparallel and vertical velocities by assuming there are no significant variations in the fault-perpendicular movement over the aperture of the InSAR data coverage. The assumption on negligible faultperpendicular deformation needs clarification. The contractional strain of the central SAF varies spatially. It has been estimated that the fault-perpendicular velocity variation is approximately 1.7 mm yr −1 for a GPS network spanning a 60-80 km wide distance across the central SAF (Rolandone et al. 2008). This variation, however, should be much smaller for the near-fault region (∼20 km) imaged by InSAR (Fig. 5).
Next we smoothed the fault-parallel velocity data by binning the data at irregular intervals. The intervals were chosen to keep dense sampling near the fault and gradually decreases to sparser sampling away from the fault. We estimated a mean and a standard deviation by averaging each bin. Assuming that each point is an independent observable that contains a deterministic signal and a random noise, the averaging reduces the variance of the observables by a factor of √ n, n being the number of data point in each bin. The estimation of the standard deviation (i.e. variance) is as important as the mean because they determine the bounds of the moment rate functional (Section 4.3).
The reduced fault-parallel residual velocity profile from InSAR across the central Creeping Section is shown in Fig. 6. There are 44 data points that cover 20 km on either side of the fault. The residual velocity variation (∼1-5 mm yr −1 ) is distributed from 0 to 10 km distance away from the Creeping Segment. The residual velocity profile recovers off-fault deformation associated with the moment deficit of the main creeping fault. This data was used in both the integral method and the moment-bounding inverse method.
The geodetic observations, shown in Figs 3 and 5, were projected into Cartesian coordinates. We examined the residuals in the best- Figure 6. Filtered residual velocity from InSAR (circles with error bars) and the model predictions (solid and dashed line). This profile is a transect crossing the central creeping section (Fig. 5). The uncertainty of the velocity profile is derived from the data statistics (see text for details). The solid black line corresponds to the best-fitting model and the dashed lines correspond to the model that enforces a bound to the moment. These moment bounds with 50 per cent confidence interval are derived from Fig. 7. The blue data points are the outliers that were removed during modelling. fitting model to detect data outliers. The outliers were defined as those residuals that are greater than 1.5 standard deviations for InSAR data and for GPS velocities respectively. InSAR data points that are more than 20 km away from the fault were also removed because these data are high-pass filtered at 40 km wavelength. The reduced numbers of observables are 3577 for ALOS, 2636 for ERS, and 213 for GPS.

R E S U LT S
We applied both the integral method and the moment-bounding inverse method using the geodetic data over the Creeping Section along the San Andreas fault. First, we used the analytical formulation of eq. (5) to derive the moment accumulation rate M L and its bounds for a 40 km wide transect crossing the centre portion of the Creeping Section. Next we applied the moment-bounding inverse method to analyse both the central Creeping Section and the along-strike variations of the moment accumulation rate from Parkfield to San Juan Bautista. Finally, the results obtained by these two different methods were compared against each other and with published results.

Integral method on the Central Creeping Section
According to eq. (5), the short wavelength component of the moment accumulation rate can be approximated by the integral of the product of the distance and the associated surface velocity averaged over a finite width 2W. We applied eq. (5) to the InSAR 2-D profile of data shown in Fig. 6. In Fig. 6, there are data outliers between −5 and −10 km distances but there are no known quaternary faults near this area (Jennings 1975). We interpret these outliers to be due to the noise in the InSAR data rather than a new structure accommodating part of the plate motion by creep. This argument is supported by the fact that there is no discernable velocity gradient shown by the GPS velocity data 5-10 km west of the San Andrea fault (Fig. 3).
We set up the NNLS problem (eq. 6) to solve for the best-fitting slip-rate model with the original data points. The residuals from the best-fitting model are carefully examined. Next, we removed outliers whose residuals are greater than 1.5 times the standard moment accumulation rate moment accumulation rate deviations. There are also large outliers within 1 km from the fault zone, possibly due to local effects. Linear interpolation is used to fill in the gap caused by the outlier removal. We used the cleaned data set to derive the moment accumulation rate estimates. Fig. 7 shows the moment accumulation rate estimates M L and the uncertainties as a function of the integration width W. The moment accumulation rate increases with W from 0 to 10 km, and then it decreases gently from 10 to 20 km. The behaviour of the moment rate functional exhibits the modulation effects determined by the filter that was applied to the velocity field (Appendix A). The moment accumulation rate M L will gradually diminish as W continues to increase. Therefore, we chose an optimal W based on the width of the filter. We took the value at W = 20 km, 2.7 × 10 15 Nm yr −1 km −1 , as the final estimation of M L . We estimated the bounds of moment accumulation rate M L by adding and subtracting one standard deviation from the mean. Because of the relatively large noise present in this limited data set, the upper bound of M L is 5.97 × 10 15 Nm yr −1 km −1 and the lower bound is close to zero.
It should be noted that the moment accumulation calculated by this method only represents the short-wavelength component. It is a lower limit of the true M L . However, if we assume that there is no moment accumulation at long-wavelengths, the estimate of M L can be deemed the true M L . We have investigated this assumption by comparing the estimation obtained by simple integration to the bestfitting value obtained by the 3-D elastic dislocation model inverted from both the GPS and InSAR velocity data. The moment accumulation rate from the integral method is 2.7 × 10 15 Nm yr −1 km −1 and a similar value of 2.5-3.4 × 10 15 Nm yr −1 km −1 is obtained by the elastic dislocation model at the same location of the fault (see Section 4.3). These two estimates are surprisingly close to each other. As shown by previous studies, the elastic dislocation model constrained by the GPS velocity data is sensitive to the moment deficit in the upper crust (Rolandone et al. 2008;Ryder & Bürgmann 2008). The filtered InSAR velocity is relatively sensitive to the moment deficit at a shallower depth. The fact that the estimations from the two different approaches are similar to each other suggests that there is a moment deficit in the upper crust of the Creeping Section.

Moment-bounding inversion on the Central Creeping Section
We present the moment-bounding inversion method in 2-D by applying it to the InSAR velocity profile shown in Fig. 6. We used eq. (6) to calculate the upper and lower bounds of the moment accumulation rate. We also calculated the measure of the misfit by employing the χ 2 statistics for each of the plausible moment rates. For an underdetermined NNLS problem, the expected χ 2 misfit for the best-fitting model should be n = o−f, and the variance 2n, where o is the number of the observables and f is the degree of freedom in the system (Parker & Song 2005). We chose f to be 2 since the simplest 2-D dislocation model is a boxcar function and it has two degrees of freedom (the slip rate and the depth of the creeping zone). For example, there are 37 data points after removing outliers (Fig. 6), thus n = 35. Fig. 8 shows the probability of the misfit due to random chance as a function of the moment rate per unit length of the fault. For a probability of 90 per cent, the chi-squared statistics gives the expected χ 2 0 = 46.06. The intersection of the probability curve and the 90 per cent probability line defines the bounds (Johnson et al. 1994).
As shown in Fig. 8, the current InSAR measurements and its noise characteristics restrict the moment accumulation rate between 0.15 and 5.6 × 10 15 Nm yr −1 km −1 with a probability of 90 per cent. The trough of the misfit curve indicates the best-fitting moment accumulation rate, which is 2.75 × 10 15 Nm yr −1 km −1 . For the best-fitting model, we obtained a minimum chi-squared misfit of 18.8, which is acceptable at a probability of 99.43 per cent. The bounds are moderately sensitive to the smoothing weights used in the inversion. If the smoothing weights change from 1 to 10, the bounds are reduced to between 0.3 and 5.1 × 10 15 Nm yr −1 km −1 .
To check if the model is in agreement with the actual data, we show the predicted velocities from the lower bound, the best fit, and the upper bound models in Fig. 6. The bounded models closely envelop the InSAR velocities and their standard deviations.
We validated the InSAR-derived velocity model with the GPS observations. In order to make comparisons with the GPS data in a straightforward way, we combined the prediction from the best-fitting model (solid line in Fig. 6) with a step-like velocity profile. The bounds were derived in the same way (dashed lines in Fig. 6). The reconstructed velocity field reflects the moment accumulation rate in the seismogenic zone. Fig. 9 highlights the velocity profile crossing the central Creeping Section. 14 out of 15 GPS measurements (except EADE) fall within the bounds of the model. The residual of the GPS from the best-fitting model has a RMS of 2.7 mm yr −1 with a mean of 0.6 mm yr −1 . Indeed, we found it remarkable that the upper and lower bounds deduced from InSAR measurements matches the GPS velocities closely (Fig. 9). This good agreement supports our assumption, which is that the moment accumulation rate M L measured by InSAR can represent the true M L .

Moment-bounding inversion and the along-strike variations
We established the along-strike variation of the moment accumulation rate using the moment-bounding inversion method (eq. 7). The inversion method can be summarized by eq. (7) and the details of the method are described in Section 2 of this paper. Similar to the 2-D problem (Section 4.2) we established the moment rate functional in 3-D by altering the value of M and solving a set of the constrained NNLS problems. Then we determined the best-fitting model and the upper and lower bounds of the moment by constructing Downloaded from https://academic.oup.com/gji/article-abstract/203/1/48/577323 by guest on 15 March 2020 a relationship of the Chi-squared misfits and the moment rate M (Fig. 8). The variation for the 3-D problem is that the bounds only apply to subsegments of the fault model specified by the vector a. When establishing the bounds for certain subsegments, the slip rates on other subsegments are not constrained. After each estimation, we modify vector a subsequently to constrain the moment accumulation rate at other locations along the fault strike and repeat the above procedures. The moment rate function is solved at 10 km intervals in the along-fault distance and is illustrated in Fig. 10. The creeping faults were divided into 21 subsegments along strike, precisely following the surface fault geometry. The fault model is approximately 200 km in length. The rectangular fault patches are approximately 10 km in length. We divided each subsegment into two subpatches according to depth. The width of the fault patches increase with depth so they are 3 km in the shallow layer and 9 km in the deep layer. The fault model extends from the surface to 12 km in depth. The complexity of the 3-D problem lies in inferring the bounds of the moment accumulation rate because we must determine the expected chi-squared value in a different way. The expected chisquared value should be determined by the number of observational data and the degrees of the freedom in the model. The number of observables varies for each subsegment and the sensitivity of the observables to unit slip on a subsegment decays with distance. The total number of observables no longer determines the expected chisquared value. To search for the observables that can be attributed to the slip on a specific fault patch, we found a maximum value of the predicted displacement due to unit slip on the fault subsegment, and compared the predicted displacement of the surrounding points to the maximum value. We deemed the data points whose kernel function are greater than 10 per cent of the maximum value to be the local observables that can be related to certain fault patch.
The best-fitting model is capable of reproducing the details in the observed velocity field. In Fig. 3(c), the predictions (blue vectors) match the observations (black vectors) closely, especially in the central Creeping Section where the detected deformations are small. In  Figs 5(b) and (e), our model predictions are in excellent agreement with the filtered residual velocities (Figs 5a and d) for both the ALOS and the ERS data. The InSAR observations agree with the GPS in recovering subtle but significant deviations from the freely slipping behaviour in the central Creeping Section. The misfits, after removing the best-fitting model (shown in Figs 5c and f), are sufficiently reduced. The degree of partial locking varies along-strike from Parkfield to San Juan Bautista, which is discernable from both Figs 3 and 5.
The along-strike variations of the moment accumulation rate indicate the location of the partially locked patch during the interseismic period. As shown in Fig. 10(b), the best-fitting model suggests that the moment accumulation rate reaches a minimum (0.6 × 10 15 Nm yr −1 km −1 ) near Monarch Peak. To the south, there is increased partial locking from Middle Mountain (4.1 × 10 15 Nm yr −1 km −1 moment rate) to Parkfield (7.6 × 10 15 Nm yr −1 km −1 moment rate). To the north, partial locking increases from Monarch Peak to Buck Ridge (3.4 × 10 15 Nm yr −1 km −1 moment rate). Further north, the fault branches into two fault segments: the southern Calaveras fault and the San Andreas Fault. The amount of locking is relatively low north of the branching point and it increases sharply to a fully locked fault (6 × 10 15 Nm yr −1 km −1 moment rate) near the Santa Cruz Mountain segment. The length of the detected aseismic zone is about 150 km here. In sum, the aseismic creep behaviors along the Creeping section are complex and likely reflect the complexity in the frictional properties (Moore & Rymer 2007) and the local stress field (Provost & Houston 2001). The bounds give more statistical meaning to the moment accumulation rate estimation. For example, we found an upper bound of the moment accumulation rate of 6 × 10 15 Nm yr −1 km −1 near Monarch Peak. Given the relatively large range of the lower bound obtained in the inverse problem, however, we cannot rule out the possibility of no moment accumulation at the central Creeping Section with >90 per cent probability. Future work should include reducing the data uncertainties with more accurate observations. The moment accumulation rates illustrated in Fig. 10(c) place our analysis into perspective with fully locked faults in California. To the south, the Cholame-Carrizo segment, with a long-term slip rate of 34 mm yr −1 and a locking depth of 12 km, is estimated to have a moment accumulation rate of 12.2 × 10 15 Nm yr −1 km −1 . To the north, the branching Santa Cruz Mountain segment, with a long-term slip rate of 17 mm yr −1 and a locking depth of 12 km, has moment accumulation is 6.1 × 10 15 Nm yr −1 km −1 . In comparison with nearby faults, the central portion of the Creeping Section is accumulating seismic moment at rates that are 5-23 per cent of the fully locked Carrizo segment. The moment accumulation rate analysed for the profile of the central Creeping Section is marked as a black line with an open circle in Fig. 10(b). This value is obtained from the analysis in Section 4.1. It is remarkable that the moment accumulation rate obtained by the simple integral method matches the elastic model inversion so accurately.
We compared our moment accumulation rate results with published rates for the Creeping Section and adjacent fault segments to infer the relative earthquake hazard in this region (Field et al. 2014). Smith & Sandwell (2003) estimated that the Creeping Section has low moment accumulation rate per unit length of fault 1.6 × 10 15 Nm yr −1 km −1 . Maurer & Johnson (2014) found that the 100 km long Creeping Section has a median potency rate of ∼2.5 × 10 7 m 3 yr −1 . If we average their potency rate by the fault length and assume a 30 GPa shear modulus, we obtained a moment accumulation rate per unit length of fault 0.75 × 10 15 Nm yr −1 km −1 . We note that dislocation models at the north end of the Creeping Section differ slightly depending on how the geometry of the faults is modelled at the fault junction (Rolandone et al. 2008;Ryder & Bürgmann 2008;Maurer & Johnson 2014). Previous work, as well as this study, all pointed out that in order to individually estimate the moment accumulation rate for these two closely-spaced fault strands, it is necessary to use the fault deformation models that took the long-term slip rates of both faults into consideration.

D I S C U S S I O N A N D C O N C L U S I O N S
In a previous study to recover the high resolution interseismic velocity field, we developed a model-based GPS/InSAR integration approach to combine long wavelength signal from GPS with the short wavelength signal from InSAR (Tong et al. 2013). This study explores new approaches for using the geodetic data near the faults to evaluate the moment accumulation rate. The moment accumulation rate across the central Creeping Section was derived without inversion by simply integrating the surface velocity measurements. This novel method was validated by direct comparisons to the momentbounding inversion results. The integral method developed here can be applied to dense GPS arrays or transects or to InSAR measurements in remote or inaccessible areas.
Surface creep rates have been used in previous studies to infer the degree to which the central part of the Creeping Section is partially locked (Titus et al. 2006;Johnson 2013;Tong et al. 2013;  Lienkaemper et al. 2014). As noted in previous studies, there are significant variations in the surface creep rate along the Creeping Section with a shape that looks like an inverted U-high creep rate in the centre tapering to zero creep rate at the ends (Rymer et al. 1984;Titus et al. 2006;Tong et al. 2013). The creep rate data are from a compilation of all the available measurements since the 1970s with various resolutions due to different techniques used (Tong et al. 2013). We incorporated this additional data set by extrapolating surface creep rates to derive an independent estimate of the moment accumulation rate (Fig. 10). The moment accumulation rate from surface creeping rate data was estimated by assuming the surface creep rate is uniform from the surface down to the seismogenic depth. As shown in Fig. 10, there are significant differences between our model and the interpolation result, especially near the Monarch Peak, consistent with Johnson (2013). It should be noted that although the interpolation accounts for the noise in the creep data, the apparent partial locking at Monarch Peak could also result from artefacts introduced by the interpolation schemes. This discrepancy suggests that the creep rate is changing with depth, in good agreement with previous studies on the Hayward fault Shirzaei & Bürgmann 2013) and the Parkfield segment (Murray et al. 2001). The moment reduction performed in the UCERF3 exercise attempted to account for this depth variation although there is no generally accepted theory for how aseismic slip rate should vary with depth. Our study emphasizes that it is critical to study the off-fault deformation signal when evaluating the moment accumulation rate.
Our best-fitting model indicates that the central portion of the Creeping Section is accumulating seismic moment at rates that are about 5-23 per cent of the fully locked Carrizo segment. The moment accumulation rate reaches a minimum near Monarch Peak and increases both southward and northward. If we suppose that the entire section of the creeping fault from north of Parkfield to south of San Juan Bautista (∼150 km in length) fails and the accumulated seismic moment is released in the form of one major earthquake, the moment magnitude of that earthquake could reach M w 7.5 assuming a 500-yr recurrence interval, compared to a M w 7.9 for a completely locked fault assuming the same length and the same recurrence interval.
Streaks of seismicity along the Creeping Section have been interpreted as alignments on the boundaries between locked and creeping patches, where repeating earthquakes represent weak asperities at a border between much larger locked and creeping patches on the fault plane (Sammis & Rice 2001). We calculated the total moment released by the earthquakes on the creeping fault using the earthquake catalogue from 1984 to 2011 from Waldhauser & Schaff (2008). Seismicity within 1 km from the fault trace were used. Fig. 10(a) shows the distribution of the moment release rate averaged over 27 yrs. It can be seen that the accumulated seismic moment has not been completely released by the recorded earthquakes in the last three decades. We can also quantify the seismic moment budget over a longer time period. Toppozada & Branum (2004) presented a history of M > 5.5 earthquakes in the Parkfield and the Creeping section from 1800 to 2001. They calculated the cumulative seismic moment released from 1857 to 2001 from the Parkfield to the Bitterwater region to be approximately 3.25 × 10 19 Nm. Assuming the moment accumulation rate we derived in this study represents a long-term average, we estimate that the accumulated seismic moment over the same area (from Parkfield to Bitterwater) and the same time period to be 3.9 × 10 19 Nm, which is 120 per cent of the moment release. After removing the contribution of seismic moment from historical earthquakes from 1857 to 2001 (Toppozada et al. 2002;Toppozada & Branum 2004), the moment deficit is reduced to 6.5 × 10 18 Nm (∼M w 6.5). When considering the 2004 Parkfield earthquake and its post-seismic afterslip (Johanson et al. 2006), the moment gap reduces to even a smaller value of 3.44 × 10 18 Nm (∼M w 6.3).
This moment deficit could have two possible explanations: (1) there could be significant earthquakes that were undocumented throughout the historical record, or (2) the aseismic slip rate along the Creeping Section has changed over time. Taking the post-seismic deformation after the 2004 Parkfield earthquake as an example, the moment release from afterslip reached as large as M w 6.1 (Johanson et al. 2006). Episodic creep events can also modify the moment budget depending on how deep and frequent the creep events originate (Shirzaei & Bürgmann 2013). It remains to be seen whether the afterslip following great earthquakes or transient slip during the interseismic period can fill in the seismic moment gap.
Given this small, but finite, moment gap (M w 6.3), would it be possible for a rupture to propagate through the Creeping Section? Dynamic rupture simulations that incorporate thermal pressurization weakening (Noda & Lapusta 2013) predict that a rupture can even propagate across 40 km wide sections that are fully creeping during the interseismic period. Kaneko et al. (2010) have shown that the dimensions of a low-coupling fault patch play a critical role in defining permanent barriers to large earthquakes. The likelihood of a through going rupture decreases as the width of the velocity weakening zone increases (Kaneko et al. 2010). Future research on rupture dynamic simulations could consider three scenarios for rupture across the Creeping Section using the lower bound, the best fit, and the upper bound models, respectively. Because the seismic moment of the next big earthquake is related to both the accumulation rate and the history of past ruptures (recurrence intervals), historical records and paleoseismological evidence of past earthquakes on the Creeping Section would be critical for understanding the likelihood of a through going rupture.

A C K N O W L E D G E M E N T S
We thank Takeshi Sagiya and Gareth Funning for their careful and constructive feedbacks, which helped impove our manuscript. We thank Duncan Agnew and Bob Parker for discussions on the inversion problem. This research was supported by the National Science Foundation (EAR-0838252, EAR-1147435, EAR-1424374, and EAR-1439697) and the Southern California Earthquake Center. This material is based on data services provided by the UNAVCO Facility with support from the National Science Foundation (NSF) and National Aeronautics and Space Administration (NASA) under NSF Cooperative Agreement No. EAR-0735156. We thank the European Space Agency (ESA) and the Japan Aerospace and Exploration Agency (JAXA) for providing the SAR data used in our study. The ERS data is obtained through the WInSAR data archive and the ALOS-1 data is obtained through the Alaska Satellite Facility.

A P P E N D I X : A N A LY T I C D E R I VAT I O N S O F T H E I N T E G R A L M E T H O D
The seismic moment accumulation rate M per unit length of fault L is given by the well-known formula where D is the thickness of the locked zone, μ is the shear modulus, and S is the slip deficit rate on the fault. This is the standard 'textbook' formula, although usually authors consider the coseismic moment due to coseismic slip. In the general case where the slip deficit rate s varies with depth z, the moment rate M is given by where D m is the maximum slip deficit depth. (Note we have assumed a uniform shear modulus with depth.) The objective of the following analysis is to show that the total moment rate per unit length of fault can be measured directly from geodetic data; no slip vs. depth model is needed. The only assumptions are that the strike-slip fault is 2-D and the Earth behaves as an elastic half-space. From general dislocation theory, the along-strike velocity (i.e. y-direction of velocity) as a function of distance from the fault is given by (Fig. 2) Next we estimate that the integral of this displacement, times distance from the x-origin, is proportional to the total moment rate. We call this proxy integral Q and later show how it is related to the moment rate M. We integrate to an upper limit W and then take the limit as W → ∞.
After re-arranging the order of integration, one finds The integral over x can be done analytically.
In the limit as W → ∞ the second term on the right side is zero because z has an upper bound of D m so the total integral is simply 1. Overall, we find this proxy is Comparing eq. (A2) with eq. (A7) it is clear that the geodetic moment rate can be directly related to the integral of the velocity times the distance from the origin.
In other words, the geodetic moment rate can be estimated directly from measurements of velocity as a function of distance from the fault. As a check of this result, we can use the solution for the y-velocity for uniform slip rate of magnitude S extending from the surface of the Earth to a depth D as given in eq. (A1). We know that the moment rate per length for this model is simply μSD so we just need to check that the integration of (A8) provides the same answer. The y-velocity for this model is given by For this case the moment per length is given by This integration can be performed analytically using the formula In the limit as W → ∞ the final result is This agrees with our original estimate of moment except for a factor of 2. Indeed to calculate the total geodetic moment we should integrate over both sides of the fault. The final result is The main utility of this formula is to demonstrate that geodetic measurements of y-velocity across an infinitely-long strike slip fault provide a direct estimate of the geodetic moment. It is unnecessary to attempt the unstable inverse problem to calculate slip versus depth and then integrate this function. We used eq. (5) to integrate filtered velocities from InSAR to estimate the short-wavelength component of the moment accumulation rate. To understand the effect of the filtering, we derived an analytical expression for how the Gaussian high-pass filter modulates the moment accumulation rate estimate.
From eqs (4) and (A13) we have The second part of the equation can be written as Following l(x, σ ) = 1