Shallow deformation of the San Andreas fault 5 years following the 2004 Parkfield earthquake (Mw6) combining ERS2 and Envisat InSAR

This study focuses on the shallow deformation that occurred during the 5 years following the Parkfield earthquake (28/09/2004, Mw 6, San Andreas Fault, California). We use Synthetic Aperture Radar interferometry (InSAR) to provide precise measurements of transient deformations after the Parkfield earthquake between 2005 and 2010. We propose a method to combine both ERS2 and ENVISAT interferograms to increase the temporal data sampling. Firstly, we combine 5 years of available Synthetic Aperture Radar (SAR) acquisitions including both ERS-2 and Envisat. Secondly, we stack selected interferograms (both from ERS2 and Envisat) for measuring the temporal evolution of the ground velocities at given time intervals. Thanks to its high spatial resolution, InSAR could provide new insights on the surface fault motion behavior over the 5 years following the Parkfield earthquake. As a complement to previous studies in this area, our results suggest that shallow transient deformations affected the Creeping-Parkfield-Cholame sections of the San Andreas Fault after the 2004 Mw6 Parkfield earthquake.


Transient Deformations at Parkfield Section
During the time period from 1992-2004, the creep rate spatially decreased along the SAF: from about 2 cm/yr at the Creeping section, to 1.4 ± 0.3 cm/yr at the Northwestern Parkfield section, to 0.6 ± 0.3 cm/yr at the southeastern Parkfield section, and finally to ~0 cm/yr at the Cholame-Carrizo section 2,3,11 . Additionally, two-colored electronic distance meters and InSAR time series analysis have revealed a temporal change of the creep rate along the Parkfield section during the same time period (notably an episodic creep acceleration between 1999 and 2000 at the location of the PKEQ epicenter 11 ). These transients creep indicate that the Parkfield section experienced an irregular spatio-temporal strain release during the last decade preceding the PKEQ 11,13 .
ScientiFic RepoRtS | (2018) 8:6032 | DOI: 10.1038/s41598-018-24447-3 The post-seismic relaxation following the PKEQ has been documented through the use of global positioning systems (GPS), InSAR, and seismological analyses. The seismicity distribution indicates that the post-seismic relaxation affects a larger area than the co-seismic rupture trace, including the Creeping and Cholame sections. Turner et al. 14 noticed a relative increasing rate of micro-seismicity (about 20%) between 2004 and 2007 with respect to the mean rate calculated over the period from 1986-2011 along the first 20 km of the Creeping section starting from the 1966 Parkfield Mw 6 epicenter (Fig. 1). This is often associated with an increase of the creep rate at depth 15,16 . In addition, the PKEQ has triggered non-volcanic tremors (NVTs, Fig. 1) at the transition between the Parkfield and Cholame sections [17][18][19] , which are associated with slow slip events occurring at 20 km depth 17 , or viscoelastic relaxation of the lower crust in response to the PKEQ 18 . These studies show that the PKEQ induced long-term perturbations of crustal properties in the SAF zone at depth [17][18][19][20] .
Geodetic studies also suggest that, following the PKEQ, the Parkfield section exhibits a complex spatio-temporal behavior, including both after-slip mechanisms at shallow depths (between 0 and 20 km) and viscoelastic relaxation at larger depths (more than 20 km). Moreover, the equivalent seismic moment released by these mechanisms over the 5 years following the PKEQ is more than twice the one released during the PKEQ [21][22][23][24][25] . From the comparison between the average fault slip rate before and after the PKEQ derived from InSAR and GPS (between the 1986-2004 and 2006-2011 periods) and seconded with numerical models, Barbot et al. 12 detected the presence of creep at depth south of the PKEQ epicenter at the transition between the Parkfield and the Cholame sections, which might indicate the presence of a soft barrier at depth, acting as an efficient arrest to earthquake ruptures for the Mw6 earthquakes but allowing propagation of larger ruptures.
Significant insights into the fault properties have been provided from the analysis of the spatio-temporal evolution of the transient deformations along the Creeping-Parkfield-Cholame sections. In complement to previous studies in this area, here we intend to better characterize the spatio-temporal distribution of the shallow fault displacement field after the PKEQ between 2005 and 2010.

Methods and Data Processing
We constructed both ERS2 and Envisat interferograms from raw Synthetic Aperture Radar (SAR) data, aiming at producing mean deformation rates on annual time-spans. Processing is done via the GAMMA Software 26 . We selected the descending track 256 for both ERS2 and Envisat (C-band, i2 mode, polarization VV, frame n°2889, Fig. 1). Barbot et al. 12 processed Envisat data (track 256, frame 2889) to produce a mean deformation velocity map around the Parkfield section between 2006 and 2010. In this study we focus on a longer period (2005-2010) with a larger dataset, including additional Envisat ASAR and ERS2 acquisitions. Our dataset is composed of 37 raw images from ERS2 and 25 raw images from Envisat ASAR (Fig. 2). To increase the temporal data sampling, we stack unwrapped interferograms jointly from ERS2 and Envisat data. Then, we derive the evolution of the along-fault creep distribution per sub-period instead of a full-period times-series analysis, such as the one proposed by de Michele et al. 11 . We co-register both ERS2 and Envisat Single Look Complex images upon one unique geometric model. Then we compute ERS2 and Envisat interferograms separately. Once the interferograms are computed and multilooked (2 in range and 10 in azimuth), orbital and topographic fringes are removed using precise orbital information and Shuttle Radar Topography Mission (SRTM) at 30 m resolution per pixel.
Two major biases can affect SAR interferograms when dealing with tectonic analyses. Firstly, non-stationary tropospheric SAR signal delay can add an undesired InSAR signal component. This delay can be correlated to topography, in which case it can be modeled and removed by linear regression with a digital elevation model (SRTM). Data stacking is also an effective way to reduce the non-topography-related atmospheric noise. Secondly, random temporal changes on the surface of the Earth can reduce the signal to noise ratio (snr), which we refer to as SAR signal coherence. Therefore, a spatial filter is applied to the wrapped interferograms to increase the snr 27 and facilitate the interferogram unwrapping procedure performed using the minimum cost flow algorithm 28,29 . We mitigate the atmospheric bias for each interferogram by estimating a linear phase trend with respect to the elevation model that we remove from each interferogram 30 . We initially calculate about 400 interferograms, most of which present low signal coherence, either due to temporal changes on the surface or to geometric signal decorrelation. As a further step, to avoid the use of low signal coherence interferograms, we perform a visual inspection. Among the total number of computed interferograms, we select 43 coherent interferograms that we use for further analyses. This significant loss of data is mainly due to poor signal coherence, mainly generated by random temporal changes in the vegetation cover, and/or persistent widespread atmospheric bias, notably at the location of the Creeping section area.
We derive the spatio-temporal evolution of the ground velocity distribution along the fault trace by stacking unwrapped interferograms from both ERS2 and Envisat. We define 5 periods of interest depending on the data availabilities (Fig. 2). We present 5 different data stacks. We calculate: (1)  We use the COSI-Corr 31 «Stacking Profile tool» to extract the velocity offsets across the fault from the stacked interferograms. The COSI-Corr Stacking Profile tool performs a linear regression on the InSAR stack profiles at each side of the fault trace, yielding one offset value at each measured location. We took measurements at less than 1 km apart from the fault trace all along the fault line as suggested in de Michele et al. 11 for this specific area (~2 km across fault baseline equivalent). The ground velocity offsets across the fault are then projected along strike and compared with creep-meter values over the same time period (available from USGS website). Creep-meter devices consist of two piers spaced 30-60 m apart and located on opposite sides of the fault connected to each other by a wire 1,32 . The cumulative horizontal along strike displacement is derived from the time evolution of the angle of the wire from the strike of the fault. From this cumulative displacement, we extracted creep rate per sub-period (the same sub periods as the interferogram stacks). As a simplification, we assume that the fault displacement at less than 1 km from the fault line is mainly due to pure strike-slip displacement. This assumption is roughly corroborated from a preliminary analysis of the available GPS time series. The velocity distribution along the fault derived from InSAR is converted from rad/yr (LOS) to cm/yr (strike-slip dextral motion). This assumption allows us to make the comparison between InSAR and creep-meters easier.
As a final step, all data are referenced to station P538 (set to 0 cm/yr), as P538 station is both close to the fault line and steady through time (see top picture, Fig. 3). . Indeed, we highlight that some sub-period stacks are more affected by atmospheric bias and noise because they are constructed from a lower number of interferograms. This can be seen around the westward limit of the Creeping section (see Fig. 3). Nonetheless, the local offset measured in the near field of the fault should still be considered as an interpretable measure. This is because, while InSAR local signal coherence might affect the precision on the offset measurements, we can exclude atmospheric biases at this very local scale. We report that, the 2009-2010 stack contains an unwrapping error at the location of the 2004 earthquake epicenter. Therefore, interpretations at this specific location shall be taken carefully.

SNR analysis (Stacks overview
To further investigate the InSAR data reliability and evaluate our pure horizontal fault motion hypothesis mentioned above, we compare InSAR data and GPS for each period using a profile across the Parkfield section, more densely covered by GPS stations (Figs 3 and 4). Firstly, the GPS velocities are projected on the LOS and then converted in strike-slip geometry following the same procedure applied to our InSAR data (case-1, red circles, Fig. 4). Secondly, as a proxy to evaluate the pure strike fault motion hypothesis, we compare the InSAR in strike-slip geometry with the fault strike component of the GPS horizontal velocities (case-2, blue circles, Fig. 4). Qualitatively, from the visual comparison between InSAR and GPS (Fig. 4) one can observe the good fit between the two, especially for the 2005-2010 period. Nonetheless, some differences between GPS case-1 and case-2 can be noticed in 2005 at the location of MIDA station (Fig. 4). We speculate that this difference might be due to the vertical or horizontal fault-normal components of the ground deformation that affects our fault motion hypothesis at this location in 2005. However, we do not observe other significant differences between the GPS (red dots) and the pure along strike fault motion derived from the GPS (blue dots) for other sub-periods. Additionally, from the analysis of the root mean square error (rmse) between the InSAR profiles in Fig. 4 (Fig. 4), might be considered a good score.
The high spatial resolution (about 50 m per pixel) allows us to derive a detailed along fault evolution of the surface creep rate (Fig. 5). The across-fault offsets (within 1 km of the fault) extracted from each sub-period stack are reported along the fault, according to the distance from A to B (Fig. 5). For each sub-period, we compare each InSAR profile to the creep-meter rate and to the 1992-2004 inter-seismic along-fault velocity profile derived by de Michele et al. 11 (blue line). From the comparison between the spline interpolation (black line) and the measured offset (red cross) we can assess the snr of the systematic velocity offset extraction procedure that we used. While we observe a good snr between the offset measurements with respect to the spline interpolation for However, we observe that all the velocity profiles spatially mimic the broad SAF segmentation at the Parkfield (i.e. transitional creeping behavior) and Cholame (i.e. locked behavior) sections.
Finally, to better visualize the creep evolution with time, we split the fault trace into 10 subsections, each 10 km long (Fig. 6). We referenced them from "a" to "j" along the Creeping, Parkfield, and Cholame sections of the SAF. For each subsection and sub-period, we focus our attention on the comparison between the mean creep rate evolution with respect to the inter-seismic references from Ryder et al. 3 for the Creeping section and from de Michele et al. 11 for the Parkfield and Cholame sections. As shown in Figs 5 and 6, the surface creep rate along the three sections of the SAF exhibits a complex spatio-temporal evolution following the PKEQ.

Observations for each sub-period and sub-sections
The Creeping section. Northwest of the Parkfield section, the 20 km of the Creeping section included in our analysis experienced an irregular evolution of the surface creep rate following the PKEQ (Figs 5 and 6ab). In  Creeping section. Moreover, during 2008, we observe important lateral variations in the creep rate distribution between Middle-Mountain and Gold-Hill (Fig. 5).
The Cholame section. In 2005 at the Parkfield-Cholame transition section, surface creep sharply decreases from 5 ± 1.2 cm/yr at the Gold-Hill location to ± 0.5 cm/yr at km 60 (Fig. 5). We observe a similar behavior in 2006-2007 (5 ± 0.5 cm/yr at Gold Hill, 1 ± 0.5 cm/yr at kilometer 60), but the transition zone seems to extend a few km southeastward compared to 2005 (Fig. 5). We infer that the surface creep has decayed toward the Cholame

Discussion
From the combined stacks of ERS and Envisat SAR interferograms, we measure the transient surface creep at the Creeping, Parkfield, and Cholame sections of the SAF during the 5 years following the 2004 earthquake. Based on previous studies, such transients in the surface creep signal were not detected during the 10 years prior to the PKEQ. These observed perturbations are characterized by both high spatial variability (with large creep values up to ~7 ± 1.5 cm/yr in 2005) and a complex temporal evolution with oscillatory like behavior. They extend from Slack-Canyon (Creeping section) to the first 10 km of the Cholame section ( Fig. 5a-g), making the total portion of the SAF perturbed after the PKEQ about twice the fault rupture length reported by Rymer et al. 33 for the coseismic. This observation is in agreement with the work of Turner et al. 14  We speculate that this oscillatory spatio-temporal behavior of the SAF might be due to a complex combination of PKEQ postseismic lower crust relaxation mechanisms in addition to the rate-and-state frictional fault behavior in response to the local or distant tectonic and/or non-tectonic sources of stress-change that radiate along the SAF. However, the relation between the shallow fault motion and deep fault motion is not yet clear and it still represents a scientific debate. The geological and structural complexity of the shallow depth crust shall be the set of additional phenomena that should be investigated further before any conclusion can be made on the fault properties from these data. Notably, we do not observe, for the central Parkfield segment, a convergence between the creep measurements of this study and the mean creep rate distribution documented by de Michele et al. 11 , even in 2009-2010, almost 6 years after the PKEQ. It is difficult to conclude if that is due to a long term change of the fault frictional properties or if this lasting transient creep is just part of more complex seismic cycle strain loading/release processes not yet fully documented. There is the possibility that these shallow transients along the analyzed segments might be somehow related to the fault at depth. These relations should be investigated in further studies. Nonetheless, the mapping of the transient deformations on the SAF along these three sections might bring constraints on the charge-discharge cycle of the fault system and perhaps give some hints to the understanding of the Parkfield earthquake sequence variability.
Creep meter arrays initially appear to be good candidates for such monitoring, as they present both a long time recording (from the 60 s to now) and high sampling rate of the ground deformation. However, while we observe that the midterm spatio-temporal trend between InSAR and the creep meters is similar (see Fig. 5, 2005-2010 period), when it comes to short term measurements, creep meter estimates might become less reliable in terms of fault motion. As shown in the San Francisco Bay area faults, creep meters might undergo short term rate changes due to various sources of environmental changes 34 . We also speculate that this difference might be exacerbated by the difference between the creep-meter across fault baseline scale (~60 m) and the ~2 km across fault baseline in InSAR (see Fig. 4). In other words, despite the ability of creep-meter devices to sample the fault creep rate with a high frequency, they are limited to fully constraining the amplitude of the surface creep rate for short periods and for small across fault deformation gradients compared to the InSAR or GPS techniques. The use of the combined present day SAR missions, such as the Sentinels 1A/1B (and future missions such as the Sentinels 1D and 1E), thanks to their high revisit time and large spatial coverage, might fulfill the monitoring needs as a complement to ground measurements from field investigations and permanent GPS stations.
In this study we made the hypothesis of purely strike slip horizontal fault motion to favor the comparison between InSAR and the creep-meter array. This hypothesis is corroborated by the good agreement between fault strike GPS velocities components and InSAR data. However, previous GPS based studies over the central portion of the SAF highlighted complex 3 dimensional deformations patterns around the fault 2 . Consequently, this hypothesis needs to be verified further, maybe using recent InSAR based technics, such as Multiple Aperture Interferometry, and combining ascendant descendent InSAR acquisitions.

Conclusions
By combining ERS2 and Envisat interferograms, we could improve the spatio-temporal sampling of the surface creep evolution for the Creeping-Parkfield-Cholame sections of the SAF, between 2005 and 2010. This methodological improvement, along with the high spatial resolution of the SAR interferograms, allowed us to estimate a map of the near field SAF surface displacement after the PKEQ. These measurements suggest the presence of transient surface deformations on this complex section of the SAF that contrast both in amplitude and distribution with what was documented during the 10 years prior to the PKEQ with InSAR. We observe that, after the PKEQ, the shallow creep rate on the Creeping, Parkfield, and Cholame sections has been highly affected from 2005 to 2010, with surface creep values reaching up to 7 ± 1.5 cm/yr with an oscillatory like spatio-temporal evolution that reveal a complex and not steady behavior of the SAF on the Creeping, Parkfield, and Cholame sections. In these terms, the measurement of the near field surface fault motion with a high spatio-temporal resolution appears as the first step in the improvement of our understanding of the geophysical process driving the surface faults kinematics behavior. This leads to the fact that monitoring both the long-term shallow fault motion and transient surface creep rate all along the Creeping, Parkfield, and Cholame sections still presents a primary scientific interest. In these terms, future geodetic based studies in this area shall benefit from the use of combined present day SAR missions, such as the European Sentinels 1A/1B missions (and future missions such as the Sentinels 1D and 1E).