\documentclass[10pt]{article}

\usepackage{fullpage}
\usepackage{setspace}
\usepackage{parskip}
\usepackage{titlesec}
\usepackage[section]{placeins}
\usepackage{xcolor}
\usepackage{breakcites}
\usepackage{lineno}
\usepackage{hyphenat}





\PassOptionsToPackage{hyphens}{url}
\usepackage[colorlinks = true,
            linkcolor = blue,
            urlcolor  = blue,
            citecolor = blue,
            anchorcolor = blue]{hyperref}
\usepackage{etoolbox}
\makeatletter
\patchcmd\@combinedblfloats{\box\@outputbox}{\unvbox\@outputbox}{}{%
  \errmessage{\noexpand\@combinedblfloats could not be patched}%
}%
\makeatother


\usepackage{natbib}




\renewenvironment{abstract}
  {{\bfseries\noindent{\abstractname}\par\nobreak}\footnotesize}
  {\bigskip}

\titlespacing{\section}{0pt}{*3}{*1}
\titlespacing{\subsection}{0pt}{*2}{*0.5}
\titlespacing{\subsubsection}{0pt}{*1.5}{0pt}


\usepackage{authblk}


\usepackage{graphicx}
\usepackage[space]{grffile}
\usepackage{latexsym}
\usepackage{textcomp}
\usepackage{longtable}
\usepackage{tabulary}
\usepackage{booktabs,array,multirow}
\usepackage{amsfonts,amsmath,amssymb}
\providecommand\citet{\cite}
\providecommand\citep{\cite}
\providecommand\citealt{\cite}
% You can conditionalize code for latexml or normal latex using this.
\newif\iflatexml\latexmlfalse
\providecommand{\tightlist}{\setlength{\itemsep}{0pt}\setlength{\parskip}{0pt}}%

\AtBeginDocument{\DeclareGraphicsExtensions{.pdf,.PDF,.eps,.EPS,.png,.PNG,.tif,.TIF,.jpg,.JPG,.jpeg,.JPEG}}

\usepackage[utf8]{inputenc}
\usepackage[ngerman,english]{babel}



\usepackage{float}








\begin{document}

\title{Integrated Sentinel-1 InSAR and GNSS time-series along the San Andreas
fault system}



\author[1]{Xiaohua Xu}%
\author[2]{David T. Sandwell}%
\author[3]{Emilie Klein}%
\author[4]{Yehuda Bock}%
\affil[1]{University of California, San Diego}%
\affil[2]{UCSD}%
\affil[3]{2Laboratoire de géologie - CNRS UMR 8538, Ecole normale supérieure - PSL University, Paris, France  }%
\affil[4]{Scripps Institution of Oceanography}%


\vspace{-1em}



  
  \date{May 22, 2024}


\begingroup
\let\center\flushleft
\let\endcenter\endflushleft
\maketitle
\endgroup





\selectlanguage{english}
\begin{abstract}
Measuring crustal strain and seismic moment accumulation, is crucial for
understanding the growth and distribution of seismic hazards along major
fault systems. Here we develop a methodology to integrate 4.5 years
(2015 - 2019.5) of Sentinel-1 Interferometric Synthetic Aperture Radar
(InSAR) and continuous Global Navigation Satellite System (GNSS) time
series to achieve 6 to 12-day sampling of surface displacements at
\textasciitilde{}500 m spatial resolution over the entire San Andreas
fault system (SAFS). We decompose the line-of-sight InSAR displacements
into three dimensions by combining the deformation azimuth from a
GNSS-derived interseismic fault model. We then construct strain rate
maps using a smoothing interpolator with constraints from elasticity.
The resulting deformation field reveals a wide array of crustal
deformation processes including: on- and off-fault secular and transient
tectonic deformation; creep rates on all the major faults; and vertical
signals associated with hydrological processes. The strain rate maps
show significant off-fault components that were not captured by
GNSS-only models. These results are important in assessing the seismic
hazard in the region.%
\end{abstract}%



\sloppy


\textbf{Integrated Sentinel-1 InSAR and GNSS time-series along the San
Andreas fault system}

\emph{Xiaohua Xu\textsuperscript{1}, David T.
Sandwell\textsuperscript{1}, Emilie Klein\textsuperscript{2} and Yehuda
Bock\textsuperscript{1}.}

\emph{\textsuperscript{1}Institute of Geophysics and Planetary Physics,
Scripps Institution of Oceanography, University of California San Diego,
La Jolla, CA, USA}

\selectlanguage{ngerman}\emph{\textsuperscript{2}Laboratoire de géologie - CNRS UMR 8538, Ecole
normale supérieure - PSL University, Paris, France}

\emph{Abstract}

Measuring crustal strain and seismic moment accumulation, is crucial for
understanding the growth and distribution of seismic hazards along major
fault systems. Here we develop a methodology to integrate 4.5 years
(2015 - 2019.5) of Sentinel-1 Interferometric Synthetic Aperture Radar
(InSAR) and continuous Global Navigation Satellite System (GNSS) time
series to achieve 6 to 12-day sampling of surface displacements at
\textasciitilde{}500 m spatial resolution over the entire San Andreas
fault system (SAFS). We decompose the line-of-sight InSAR displacements
into three dimensions by combining the deformation azimuth from a
GNSS-derived interseismic fault model. We then construct strain rate
maps using a smoothing interpolator with constraints from elasticity.
The resulting deformation field reveals a wide array of crustal
deformation processes including: on- and off-fault secular and transient
tectonic deformation; creep rates on all the major faults; and vertical
signals associated with hydrological processes. The strain rate maps
show significant off-fault components that were not captured by
GNSS-only models. These results are important in assessing the seismic
hazard in the region.

\emph{Plain Language Summary}

Seismic hazard models rely on accurate measurements of small
interseismic motion over large space on the Earth's crust. Traditional
geodetic models based on Global Navigation Satellite System (GNSS) data
cannot resolve small scale deforming patterns, mainly due to expensive
and limited station deployment. Interferometric Synthetic Aperture Radar
(InSAR) has become the emerging tool for mapping the surface
deformation, with its advantages of low-cost and full-coverage. Yet
InSAR measurement, compared to GNSS, comes with larger bias from the
atmospheric noise, especially over length scales greater than 80 km.
Here we combined the two methods to resolve fine spatial scales and
achieve high accuracy. ~Our results are presented as deformation
time-series over the entire San Andreas fault system. From these
deformation time series we have estimated fault creep rates and strain
accumulation. One important finding is that there is significant
off-fault strain, though we suspect they are mainly due to hydrological
processes. These results will advance our knowledge of the earthquake
cycle, strain/moment accumulation, and the associated seismic hazards.

\emph{Key points}

\begin{itemize}
\tightlist
\item
  A practical approach is developed to integrate Sentinel-1 InSAR and
  GNSS time-series over the entire San Andreas fault system
\item
  The product is used to estimate fault creep and three components of
  horizontal crustal strain rate which shows notable off-fault portion
\item
  Challenges remain in splitting tectonic and hydrologic sources and
  whether hydrologic strain will increase seismic hazards
\end{itemize}

1. Introduction

An improved understanding of how energy and seismic moment accumulate in
the crust and upper mantle along major plate boundaries is essential for
forecasting the size and timing of major earthquakes
{[}\emph{Smith-Konter \& Sandwell} , 2009; \emph{Field et al.,}
2015;\emph{Weiss et al.,} 2020; \emph{Rollins et al.,} 2020{]}. Recent
studies have shown that most damaging earthquakes occur in areas where
the crustal strain rate exceeds 100 nanostrain/yr {[}e.g., \emph{Elliot
et al} ., 2016.; \emph{Zeng et al.} , 2018; \emph{Bayona et al.,}
2021{]}. Many of these areas are heavily populated and have had major
destructive earthquakes. Moreover, one of the largest uncertainties in
California earthquake hazard models (i.e., UCERF-3 {[}\emph{Field et al}
., 2014;\emph{Field et al.} , 2015{]}) is the amount of plate boundary
deformation that is accommodated by off-fault strain and whether this
strain is accumulating as elastic or as plastic deformation. Therefore,
accurate strain rate measurements are needed to improve earthquake
forecasts. Achieving an ideal 100-nanostrain/yr accuracy at a 10-km
resolution (i.e., a typical fault locking depth in California) requires
a horizontal velocity model that has an accuracy of 1 mm/yr. Moreover,
moderate earthquakes, fault creep, and other transient processes produce
temporal variations in strain rate that commonly exceed
100-nanostrain/yr {[}\emph{Holt and Shcherbenko} , 2013; \emph{Klein et
al.,} 2019{]}. Currently continuous GNSS measurements can provide vector
deformation better than the required 1 mm/yr accuracy but not with 10 km
spatial resolution. InSAR provides very high spatial resolution but
cannot achieve the 1 mm/yr accuracy, mainly due to the perturbations
from atmospheric noise {[}\emph{Emardson et al.} , 2003{]}. In addition,
current InSAR systems provide only two components of surface deformation
and thus cannot uniquely distinguish between horizontal and vertical
strain {[}\emph{Shen and Liu} , 2020{]}. Here we develop a path to
achieving the time-dependent strain rate mapping objective by combining
4.5 years of measurements from InSAR and GNSS along the San Andreas
fault system (SAFS).

The accuracy and spatial resolution of the current strain rate models
derived from GNSS velocities can be assessed by comparing results from
various groups. An accuracy analysis was performed as part of the
developing SCEC Community Geodetic Model (CGM-V1) {[}\emph{Sandwell et
al.,} 2016a{]}. The 17 models were taken from previous publications
{[}\emph{Zeng \& Shen,} 2016; \emph{Shen et al.,} 2015;
Smith-\emph{Konter \& Sandwell,} 2009; \emph{Tong et al.,}
2014;\emph{Tape et al.,} 2009; \emph{Petersen et al.,} 2008;
\emph{Petersen et al.,} 2014; \emph{Platt \& Becker,} 2010;
\emph{McCaffrey,} 2005;\emph{Loveless \& Meade,} 2011; \emph{Hackl et
al.,} 2009;\emph{Parsons,} 2006; \emph{Parsons et al.,} 2013;
\emph{Kreemer et al.,} 2014; \emph{Flesch et al.,} 2000; \emph{Field et
al.} , 2014;\emph{Sandwell \& Wessel} , 2016{]}. The mean and standard
deviation of the 10 ``best'' well-correlated models is shown in Fig. 1a.
Note the standard deviation (Fig. 1b) commonly exceeds 50 nanostrain/yr
especially above the major faults where the uncertainties can exceed 100
nanostrain/yr. These significant deviations among the models are not due
to inaccuracies in the GNSS data but to the incomplete spatial sampling
of the GNSS stations, which is typically 10-20 km in California
{[}\emph{Wei et al.,} 2010{]}.

To further characterize this lack of spatial resolution of the strain
rate field, we analyzed the strain rate models by computing cross
spectra in their overlapping region in Southern California.
Radially-averaged cross spectra were computed between every pair of
models using Generic Mapping Tools (GMT) {[}\emph{Wessel et
al.,}2019{]}. As illustrated in Fig. 1c, there is a large variation in
the coherence between these models. Most models agree well at long
wavelengths, but generally disagree at short wavelengths, except those
having very similar, or identical, fault models like \emph{Shen et al} .
{[}2015{]} and \emph{Zeng \& Shen} {[}2016{]} or \emph{Tong et
al.}{[}2014{]} and \emph{Smith-Konter \& Sandwell} {[}2009{]}. The
disagreements are due to different physical modeling approaches, assumed
fault geometries, and slightly different GNSS velocity data sets. For
most pairs, the coherence is high at very long wavelength and decreases
to zero coherence at \textasciitilde{}10 km. The 0.2 coherence threshold
of the median of all the cross spectra is located at 30-40 km
wavelengths.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image1/image1}

\end{center}
\end{figure}

Fig. 1. Second invariant of the horizontal strain rates from GNSS models
{[}\emph{Sandwell et al.} , 2016a{]}. a) Average of 10 ``best'' models
based on a spatial correlation analysis. b) Standard deviation of the 10
models. c) Radially-averaged cross spectra between all pairs of SCEC
CGM-V1 strain rate models, performed in the area denoted by the white
box in a). The mean and median of all the spectra have a 0.2 coherence
at wavelengths of 30-40 km. Our objective for accurate strain-rate
measurement is to maintain at least 0.2 coherence among the models, at a
10 km spacing, when integrated with InSAR.

2. Methods

Our InSAR/GNSS integration approach is an extension of standard
published methods (e.g., \emph{Tong et al.} , 2013; \emph{Weiss et
al.,}2020) although in addition to secular velocity, we also calculate
line-of-sight (LOS) displacement time series {[}e.g., \emph{Neely et
al.,} 2019{]}. The GNSS weekly displacements were derived by means of a
median filter {[}\emph{Klein et al.} , 2019{]} of daily time series
estimated as part of a NASA MEaSUREs project {[}\emph{Bock et al.} ,
2016{]}. Moreover, we use the secular velocity from a GNSS-only
interseismic model {[}\emph{Zeng \& Shen} , 2017{]} to create
semi-vertical vector InSAR time series from the LOS displacements. A
brief description of the method follows:

\begin{enumerate}
\tightlist
\item
  Gather Sentinel-1 Terrain Observations with Progressive Scans (TOPS)
  data from multiple tracks and re-assemble into common re-defined
  frames, typically 250 km by 500 km.
\item
  Geometrically co-register all SAR acquisitions and construct all
  interferograms with perpendicular baseline \textless{} 150 m and
  temporal separation \textless{} 90 days {[}\emph{Xu et al.} ,
  2017;\emph{Sandwell et al.} , 2016b{]}.
\item
  Mask bodies of water and areas of persistent low coherence regions and
  replace them with nearest-neighbors {[}\emph{Shanker \& Zebker} ,
  2009{]}. This step improves the phase unwrapping accuracy which is
  done with Statistical-Cost, Network-Flow Algorithm for Phase
  Unwrapping (SNAPHU) {[}\emph{Chen \& Zebker} , 2002{]}.
\item
  Perform elevation dependent atmospheric phase correction
  {[}\emph{Elliott et al.} , 2008{]}. Compute the difference between the
  remaining InSAR phase and projected GNSS weekly solutions
  {[}\emph{Klein et al.} , 2019{]}, interpolate this difference, filter
  at 80-km wavelength and remove this difference from each
  interferogram.
\item
  Construct time-series using a coherence-based SBAS approach integrated
  with atmospheric phase correction using common-scene stacking
  {[}\emph{Tymofyeyeva \& Fialko} , 2015; \emph{Tong \& Schimdt} ,
  2016;\emph{Xu et al.} , 2017{]}.
\item
  Subtract a horizontal GNSS velocity model {[}e.g., \emph{Zeng \& Shen}
  , 2017{]} from the time-series to create semi-vertical InSAR
  time-series.
\end{enumerate}

Since Sentinel-1 TOPS data is acquired under burst acquisition mode and
there is occasional inconsistency in data coverage, especially in the
early days of the mission, the frame boundaries in step 1) are a
compromise between spatial coverage and acquisition numbers. The total
number of interferograms generated here is 5230, connecting acquisitions
from 910 dates over 9 tracks. Enhanced spectral diversity
{[}\emph{Prats-Iraola et al.} , 2012{]} is not performed in step 2,
since it will remove an expected tectonic signal that will eventually
supply a third InSAR component {[}\emph{Li et al.,} 2021{]}. The
estimated mis-registration could be up to 2/1000 pixel/yr along SAFS and
spread across the scenes, where a constant shift from ESD is inadequate,
while the performance of bivariate approach {[}\emph{Wang et al.,}
2017{]} is yet to be evaluated. Moreover, the common scene stacking time
series approach (step 5) is capable of mitigating along-track orbital
errors by absorbing burst discontinuities, that are random in time, into
atmospheric phase screens {[}\emph{Xu et al.} , 2017{]}. The
nearest-neighbor interpolation in step 3) is implemented so phases are
allowed to vary properly along very long coastlines, and stay connected
through snowy Sierras and heavy vegetations in northern California.
Elevation dependency in step 4) is assumed as a bivariate quadratic
polynomial thus spatial variations in atmospheric contribution are
accounted for. The relatively large, 80-km wavelength filter, that is
applied to the GNSS correction for each interferogram, is sufficient to
absorb the large-scale atmospheric and orbital errors affecting the
InSAR displacements and also accommodate areas such as the Central
Valley having sparse GNSS coverage. A remove-restore approach
{[}\emph{Tong et al.} , 2013{]} using a purely horizontal secular
velocity model is not used because it is incompatible with the
significant vertical signal in Central Valley. Not only is this vertical
deformation distributed over hundreds of kilometers, but it also has a
sharp transition around the edges of the sedimentary basins. The
interpolation of discrepancies between GNSS and InSAR is adopted here,
taking advantage that though the vertical deformation changes
dramatically over a large area, the differences from the two type of
observations maybe systematic and vary slow enough in space to be well
evaluated. When the final velocity is computed, the first and last four
records are not used, mainly because the atmospheric correction approach
gains less constraints when acquisitions are non-evenly distributed.

3. Analyses, Results and Discussions

\emph{Average Velocity}

The integrated Sentinel-1 InSAR and GNSS time-series reveal a complex
time-averaged velocity field (Fig. 2, 3) over the entire San Andreas
fault system. The most prominent feature is the large-scale tectonic
motion between the two plates with sharp transitions at creeping
strike-slip faults and broader transitions at more deeply-locked faults.
Superimposed on this tectonic pattern are numerous nontectonic
deformations mostly related to hydrological and hydrothermal processes.
These non-tectonic features are highlighted in the semi-vertical
component (Fig. 2d). Our velocity maps are available as kmz files
(http://topex.ucsd.edu/gmtsar/tar/San\_Andreas\_Xu\_et\_al.kmz ) so that
the user can examine in detail the correlations between vertical
deformation and natural and man-made surface features (See Supplementary
Text S1 and Movie S1 for much more details).\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image2/image2}

\end{center}
\end{figure}

\textbf{Fig. 2.} LOS deformation velocity (descending) along the San
Andreas fault system. a) Total LOS deformation velocity with red color
denoting motion toward the satellite and blue away from the satellite. A
second color scale is added to highlight the major subsidence in the
Central Valley. b) LOS velocity from \emph{gpsgridder} {[}\emph{Sandwell
\& Wessel} , 2016{]} projected into the LOS (for illustration purposes
only). b) Velocity variation/uncertainty is the deviation of the time
series at each pixel from a linear regression fit. The velocity
discontinuities between swaths are due to the change in look angle. d)
Semi-vertical velocity is the LOS with horizontal velocity model
{[}\emph{Zeng \& Shen} , 2017{]} removed. The triangles in d) mark the
location of GNSS stations plotted in Fig 4b (black) and Fig 4c (gray).\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image3/image3}

\end{center}
\end{figure}

\textbf{Fig. 3} . LOS deformation velocity (ascending) along the San
Andreas fault system. a) Total LOS deformation velocity with red color
denoting motion toward the satellite and blue away from the satellite. A
second color scale is added to highlight the major subsidence in the
Central Valley. b) Velocity variation/uncertainty derived from linear
regression on the resulting time-series. c) LOS velocity from gpsgridder
{[}\emph{Sandwell \& Wessel} , 2016{]} projected into the LOS (for
illustration purpose only).

\emph{LOS Time Series}

We have compared the integrated InSAR/GNSS LOS time series with the
GNSS-only time series and find both excellent and poor agreements (Table
S2). Note the two time-series will not match exactly as our approach is
essentially trusting GNSS at longer wavelengths while InSAR at short
wavelengths. Fig. 4b shows several examples where there is excellent
agreement between the InSAR and GNSS time-series. The typical deviation
of GNSS daily solutions from its moving average is around 3-4 mm, while
the integrated InSAR deviates only slightly larger, around 3-6 mm (Fig.
4a). The vertical component could largely reproduce this deviation and
can be considered as the major source of discrepancy. However, there are
a number of cases where the integrated InSAR GNSS time-series fail to
match the GNSS solutions (Fig. 4c). In most cases this is a direct
result of anomalous displacements or artifacts in the GNSS displacement
time series. An extreme example is station p800 in Los Angeles showing
systematic effects due to the growth of thick vegetation nearly
concealing the GNSS antenna and causing a significant drift in the north
direction over a period of about 2 years (overlapping with the 4.5 years
of our time series) until the vegetation is cleared in May 2018 (add
figure in supplement). Station CRHS also in Los Angeles and affected by
vegetation and whose data after 2006 are suspect and should not be used
(personal communication, Chris Walls, UNAVCO NOTA engineer) -- the
station was subsequently abandoned. Station CUHS in the Cuyama Valley is
subject to significant subsidence of about 30 mm/yr with a seasonal
signature resulting in non-tectonic horizontal displacements with
amplitudes up to 30 mm. Station CJMG in the San Gabriel Mountains has
annual horizontal artifacts starting in 2017 with an amplitude of 30 mm
(peak to peak). Another example is station P789 on the San Andreas fault
in the transition zone between locked and creeping section is
experiencing uplift since 2013 of about 2.5 mm/yr but not fully
sufficient to explain its misfit with InSAR. For the purpose of this
study these GNSS data were not excluded. Rather, we allowed for
discrepancies to exist through the integration, especially when a single
station shows large mismatch, because the GNSS correction model, applied
to each interferogram, is smoothed with a robust filter that down weighs
the anomalous misfits. Thus, in cases a GNSS station starts to behave
anomalously, the InSAR time-series remain largely unaffected (Fig. 4c).
In further studies, we recommend a more conservative approach in
choosing stations for the GNSS/InSAR integration. More assessment and
results are available in Table S2 and Fig S1 in the supplements.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image4/image4}

\end{center}
\end{figure}

Fig. 4. Comparisons of InSAR LOS time series with GNSS time series
projected into the LOS. a) Histograms of deviation from GNSS moving
average (window 4 weeks). Green colors are for GNSS daily solution, blue
for integrated InSAR \& GNSS time-series, and red for the semi-vertical
component in the integrated InSAR \& GNSS time-series (horizontal GNSS
model removed, comparison made with GNSS vertical component only). b)
Subplots of InSAR and GNSS time series at 4 GNSS stations showing good
agreement after integration. Uncertainties of InSAR time-series are
standard deviations taken in 500-m boxes. c) Subplots of InSAR and GNSS
time series at 5 GNSS stations showing poor agreement after integration.

\emph{Fault Creep Rate}

Fault creep rate can be estimated from just two components of LOS
deformation if one assumes there is no fault-normal displacement
{[}\emph{Xu et al.} , 2018{]}. By making that assumption we measure
fault-parallel and vertical creep rate along each individual strand
within the San Andreas Fault system (Fig. 5, Table S1). The method,
which is similar to the \emph{Bufford \& Harsh} {[}1980{]} method, is to
fit a linear model to the LOS data over a flat area on each side of the
fault (up to 3 km) and record the offset of the linear models at the
fault location. The along-fault offset locations have a typical spacing
of 10 km and roughly follow the locations from the study of \emph{Tong
et al.} {[}2013{]}. Some areas of fault creep have a velocity step up to
1 km wide; In these cases, the linear fit is confined to the flatter
part of the LOS data on each side of the fault (e.g., Fig. 5a, at
Hayward fault). The fault creep rates estimated in this study are
presented together with a compilation by \emph{Field et al.} {[}2009{]}
(Fig. 5a). Fault vertical differential (Fig. 5b) shows no systematic
pattern along these faults. However, not accounting for the vertical
offset will bias the horizontal offset estimates. A few larger vertical
offsets (Fig 5b) are associated with subsidence at fault step-overs.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image5/image5}

\end{center}
\end{figure}

Fig. 5. Comparison of creep rate estimates along major segments of the
San Andreas fault system. a) Red circles are right-lateral creep rates
from decomposing two-look InSAR LOS estimates into fault parallel
horizontal and vertical. Detailed estimates for each segment are
available in Table S1. A compilation of creep estimates from \emph{Field
et al.} {[}2009{]} are also plotted with markers for each individual
study listed in the upper right box. Survey types are listed as:
AA-alignment array, CM-creep meter, Cult-cultural offset features,
Geod-small geodetic array, Mod-inferred from model, Tri-trilateration.
The upper left box shows an example from a descending InSAR track
covering Hayward fault, showing how the estimates are performed. b)
Fault vertical differential in the creep rate estimates. Positive is
defined as eastern side uplift (northern side up for Garlock fault).

\emph{Velocity Decomposition}

In order to compute the horizontal components of strain rate we must
first construct 3 components of surface velocity (east, north and up)
from only two components of LOS velocity. To accomplish this, we adopt
the assumption that the direction of surface displacement matches the
direction from the GNSS velocity model {[}\emph{Wright et al.} ,
2004;\emph{Tymofyeyeva \& Fialko} , 2018; \emph{Shen \& Liu} , 2020{]}.
This utilizes the local azimuth of a horizontal GNSS velocity predicted
from the interseismic slip model of \emph{Zeng \& Shen} {[}2017{]} as a
constraint to reduce the required degrees of freedom in the measurements
from three to two. The decomposed velocity field (Fig. 6) reveals
detailed spatial variations in the east-west component together with
dramatic vertical motions. Some of the variations in horizontal
component are associated with hydrological signals (e.g., inside Central
Valley). Subsidence inside the Central Valley is greater than 30 cm/yr
with other areas like Geysers, Heber and Cerro Prieto Geothermal Fields
standing out in the map. Uplifts that occur at ground water recharging
sites (Long Beach, Santa Clara Valley, etc.) are usually a direct
reflection of human activities, some of which cannot be well represented
by a single velocity (e.g., Fig 2c., GNSS station SACY). Since the
velocity is derived as a linear fit, a positive rate does not
necessarily result from an overall uplift. Due to the near-polar orbits,
InSAR satellites are generally not very sensitive to north-south motion.
Thus, the decomposed north-south component absorbs most of its
information from the GNSS model, resulting in an oversmoothed field.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image6/image6}

\end{center}
\end{figure}

Fig 6. Surface vector velocity maps for the San Andreas fault system in
the ITRF2014 coordinate system. a), b) and c) are the decomposed east,
north and vertical component respectively. d) and e) are GNSS horizontal
velocity model from \emph{Zeng \& Shen} {[}2017{]}. f) are stations
distributions, with white triangle denoting campaign sites used in the
GNSS model and red triangles denoting continuous sites used both in GNSS
model and InSAR/GNSS integration. Decomposition of InSAR GNSS integrated
velocity maps uses information of horizontal deformation direction
provided by the model shown in d) and e).

\emph{Strain Rate}

The 3-D velocity maps derived from the InSAR/GNSS combination have small
spatial scale variations that are associated with hydrological processes
and residual atmospheric error. These models must be smoothed to extract
the underlying tectonic strain rate {[}e.g., \emph{Weiss et
al.,}2020{]}. This is accomplished by using a smoothing interpolator
that has constraints from elasticity to couple the two horizontal
velocity components {[}\emph{Haines and Holt,} 1993{]} as implemented in
the program \emph{gpsgridder} in GMT {[}\emph{Sandwell \& Wessel} ,
2016{]}. We first subsample both the GNSS and InSAR/GNSS horizontal
velocity fields shown in Fig. 6 at a spacing of 2.5 km resulting in
40,016 2-D velocity estimates. The \emph{gpsgridder} algorithm uses
singular value decomposition of the elasticity Greens functions and we
select only those associated with the largest 800 eigenvalues. We then
compute the three horizontal strain rate components from the smoothed
horizontal velocities. The comparison of the GNSS only and InSAR/GNSS
for several components of strain rate (second invariant, maximum shear
and dilatation) are shown in Fig 7.

The San Andreas fault is well delineated by the concentrated high-strain
(\textgreater{} 200 nanostrain/yr). Compared to the prediction of the
GNSS-only strain rate model, the integrated InSAR/GNSS model reveals
larger strain concentration, especially at locations where there is
known surface creep (e.g., creeping section of SAFS and Hayward fault).
Maacama and Bartlett Springs both have much larger strain concentrations
in the InSAR/GNSS model although this may be a consequence of the
assumption that the InSAR deformation is parallel to the GNSS
deformation. Dilatational strain predicted by the integrated InSAR/GNSS
model shows extension in the Central Valley associated with the
subsidence from groundwater extraction. The shear component is also
larger in areas of known hydrological signals. There is a prominent
shear strain anomaly in the InSAR/GNSS strain rate maps (Fig. 7e, f)
that is to the east of and parallel to the creeping section. However,
there is not enough GNSS data in the region to determine if this is a
real strain-rate anomaly or artifact from the interpolation approach.
Overall, the strain is more widely distributed in the integrated
InSAR/GNSS model, with significant components being off-fault. The
remaining questions are how much of these are from tectonic motion and
how much from hydrological activities, and whether these off-fault
components are steady over time or just transients. Answering these will
be the key leading to accurate assessment of the strain/moment
accumulation and the associated seismic hazards.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image7/image7}

\end{center}
\end{figure}

Fig. 7. Strain rate maps for San Andreas fault system. a) is the second
invariant of strain rate derived with gpsgridder using GNSS data only.
Black and blue bars represent compressional and extensional principal
strain rates respectively, and they are clipped at 100 nanostrain/yr. b)
is same as a) but using decomposed InSAR+GNSS horizontal velocity
components. c) Differences between a) and b). d) and e) Maximum shear
strain rate (defined as\(\max\left({\dot{\varepsilon}}_{1}-{\dot{\varepsilon}}_{2}\right)/2\),\emph{Savage et al.,} 2001)
maps with f) being their difference. g) and h) Dilatation rate (defined
as\(\left({\dot{\varepsilon}}_{1}+{\dot{\varepsilon}}_{2}\right)\),\emph{Savage et al.,} 2001) with i) being their
difference. In all plots, the dashed bounding box denotes the area with
Sentinel-1 coverage from two look directions.

\par\null

5. Conclusion

The main conclusions are:

\begin{itemize}
\tightlist
\item
  Current strain rate models based only on point GNSS measurements are
  reliable at longer wavelengths (\textgreater{} 30-40 km) but have
  large differences at shorter wavelengths so they cannot achieve the
  100 nanostrain/yr accuracy threshold for assessing seismic hazard. In
  addition, there are significant transient deformation processes (e.g.,
  hydrologic, magmatic, and tectonic) having strain rate signals that
  exceed the 100 nanostrain/yr threshold needed for hazard assessment.
\item
  We developed an approach to integrate Sentinel-1 InSAR and GNSS
  time-series (4.5 years) over the entire San Andreas fault system from
  the ascending and descending look directions. This analysis was
  enabled by the frequent, high-quality, observations from the
  Sentinel-1 satellites.
\item
  A video tour of the high-resolution LOS velocity maps reveals a wide
  array of deformation processes including: active faults and stepovers;
  extraction and recharge of groundwater, petroleum and geothermal
  fluids; and continuous expansion of the surface of dry lake beds
  (Movie S1, and available at
  https://www.youtube.com/watch?v=SxNLQKmHWpY ).
\item
  The two components of average LOS velocity are used to refine and
  update estimates of creep rate along the major strands of the San
  Andreas fault system.
\item
  The two LOS components are decomposed into 3-components of velocity by
  assuming the direction of deformation matches those predicted by a
  GNSS-only velocity model.
\item
  The higher spatial resolution vector velocity maps are used to
  estimate the three components of horizontal crustal strain rate. The
  results show significant off-fault strain, yet challenges remain in
  separating the contribution from tectonic and hydrologic sources and
  whether hydrologic strain will increase seismic hazards.
\item
  Given the 20-year plus observation plan of the twin Sentinel-1
  satellites, as well as continued GNSS operations, these high spatial
  resolution, time-dependent products will continue to improve.
\end{itemize}

Acknowledgements

The authors want to thank European Space Agency for the extraordinary
data open policy on the Sentinel-1 mission and want to thank Alaska
Satellite Facility and UNAVCO for archiving the data and the precise
orbital products. This study was funded by the NASA Earth Surface and
Interior Program (80NSSC19K1043 and NNX17AD99G), the NASA MEaSUREs
Program (NNH17ZDA001N), the National Science Foundation, Office of
Advanced Cyberinfrastructure (OAC-1834807), and the Southern California
Earthquake Center (SCEC) (20074). SCEC is funded by the NSF cooperative
Agreement EAR-1033462 and USGS Cooperative Agreement G12AC20038. All
data are originally available through Sentinel Scientific Data Hub
(https://scihub.copernicus.eu) and Alaska Data Search Vertex
(https://search.asf.alaska.edu). Processed data are currently available
at https://topex.ucsd.edu/gmtsar/insargen and on the move to SCEC as
part of the Community Geodetic Model
(https://www.scec.org/research/cgm).

References

Shanker, A. P., \& Zebker, H. A. (2009). Sparse two-dimensional phase
unwrapping using regular grid methods.~\emph{IEEE Geoscience and Remote
Sensing Letters} ,~\emph{6} (3), 519-522.

Bayona, J.A., Savran, W., Strader, A., Hainzl, S., Cotton, F. and
Schorlemmer, D., 2021. Two global ensemble seismicity models obtained
from the combination of interseismic strain measurements and
earthquake-catalogue information. \emph{Geophysical Journal
International} , 224(3), pp.1945-1955.

Bock, Y., S. Kedar, A. W. Moore, P. Fang, J. Geng, Z. Liu, D. Melgar, S.
E. Owen, M. B. Squibb, F. Webb (2016), Twenty-Two Years of Combined GPS
Products for Geophysical Applications and a Decade of Seismogeodesy,
International Association of Geodesy Symposia, Springer International
Publishing, doi:10.1007/1345\_2016\_220.

Burford, R. O., \& Harsh, P. W. (1980). Slip on the San Andreas fault in
central California from alinement array surveys.~\emph{Bulletin of the
Seismological Society of America} ,~\emph{70} (4), 1233-1261.

Chen, C. W., \& Zebker, H. A. (2002). Phase unwrapping for large SAR
interferograms: Statistical segmentation and generalized network
models.\emph{IEEE Transactions on Geoscience and Remote Sensing} ,
40(8), 1709-1719.

Elliott, J. R., Biggs, J., Parsons, B., \& Wright, T. J. (2008). InSAR
slip rate determination on the Altyn Tagh Fault, northern Tibet, in the
presence of topographically correlated atmospheric
delays.~\emph{Geophysical Research Letters} ,~\emph{35} (12).

Elliott, J.R., Walters, R.J. and Wright, T.J., 2016. The role of
space-based observation in understanding and responding to active
tectonics and earthquakes. \emph{Nature Communications} , 7, p.13844.

Emardson, T. R., et al. (2003), Neutral atmospheric delay in
interferometric synthetic aperture radar applications: Statistical
description and mitigation, \emph{Journal of Geophysical Research-Solid
Earth} , 108(B5).

Field, E. H., Dawson, T. E., Felzer, K. R., Frankel, A. D., Gupta, V.,
Jordan, T. H., \ldots{} \& Wills, C. J. (2009). Uniform California
earthquake rupture forecast, version 2 (UCERF 2).~\emph{Bulletin of the
Seismological Society of America} ,~\emph{99} (4), 2053-2107.

Field, E. H., Arrowsmith, R. J., Biasi, G. P., Bird, P., Dawson, T. E.,
Felzer, K. R., \ldots{} \& Zeng, Y. (2014). Uniform California
Earthquake Rupture Forecast, Version 3 (UCERF3)---The Time-Independent
Model.\emph{Bulletin of the Seismological Society of America} , 104(3),
1122-1180.

Field, E.H., Biasi, G.P., Bird, P., Dawson, T.E., Felzer, K.R., Jackson,
D.D., Johnson, K.M., Jordan, T.H., Madden, C., Michael, A.J. and Milner,
K.R., 2015. Long-term time-dependent probabilities for the third Uniform
California Earthquake Rupture Forecast (UCERF3). \emph{Bulletin of the
Seismological Society of America} , 105(2A), pp.511-543.

Flesch, L. M., Holt, W. E., Haines, A. J., \& Shen-Tu, B. (2000).
Dynamics of the Pacific-North American plate boundary in the western
United States. \emph{Science} , 287(5454), 834-836.

Hackl, M., R. Malservaisi, and S. Wdowinski (2009), Strain pattern from
dense GPS networks, \emph{Nat. Hazards Earth Syst.} , 9., 1177-1187.

Haines, A.J. and Holt, W.E., 1993. A procedure for obtaining the
complete horizontal motions within zones of distributed deformation from
the inversion of strain rate data. \emph{Journal of Geophysical
Research: Solid Earth} , 98(B7), pp.12057-12082.

Holt, W. E., \& Shcherbenko, G. (2013). Toward a continuous monitoring
of the horizontal displacement gradient tensor field in Southern
California using cGPS observations from Plate Boundary Observatory
(PBO).~\emph{Seismological Research Letters} ,~\emph{84} (3), 455-467.

Klein, E., Bock, Y., Xu, X., Sandwell, D. T., Golriz, D., Fang, P., \&
Su, L. (2019). Transient deformation in California from two decades of
GPS displacements: Implications for a three-dimensional kinematic
reference frame.~\emph{Journal of Geophysical Research: Solid Earth}
,~\emph{124} (11), 12189-12223.

Kreemer, C., Blewitt, G., \& Klein, E. C. (2014). A geodetic plate
motion and global strain rate model. \emph{Geochemistry, Geophysics,
Geosystems} , 15(10), 3849-3889.

Li, X., J\selectlanguage{ngerman}ónsson, S., \& Cao, Y. (2021). Interseismic deformation from
Sentinel-1 burst-overlap interferometry: Application to the southern
Dead Sea fault. \emph{In review for GRL}

Loveless, J. P., and B. J. Meade (2011), Stress modulation on the San
Andreas fault by interseismic fault system interactions, \emph{Geology}
, 39(11), 1035-1038.

McCaffrey, R. (2005), Block kinematics of the Pacific--North America
plate boundary in the southwestern United States from inversion of GPS,
seismological, and geologic data, \emph{J. Geophys. Res} ., 110, B07401,
doi:10.1029/2004JB003307.

Neely, W. R., Borsa, A. A., \& Silverii, F. (2019). GInSAR: A cGPS
correction for enhanced InSAR time series.~\emph{IEEE Transactions on
Geoscience and Remote Sensing} ,~\emph{58} (1), 136-146.

Parsons, T. (2006), Tectonic stressing in California modeled from GPS
observations, \emph{J. Geophys. Res} . 111, doi:10.1029/2005JB003946.

Parsons, T., Johnson, K. M., Bird, P., Bormann, J., Dawson, T. E.,
Field, E. H., \ldots{} \& Zeng, Y. (2013). Appendix C---Deformation
models for UCERF3.

Petersen, M. D., Y. Zeng, K. M. Haller, R. McCaffrey, W. C. Hammond, P.
Bird, M. Moschetti, Z. Shen, J. Bormann, and W. Thatcher {[}2014{]}
Geodesy- and geology-based slip-rate models for the Western United
States (excluding California) national seismic hazard maps, U.S. Geol.
Surv. Open-File Rep., 2013-1293, 38 pages doi.org/10.3133/ofr20131293.

Platt J. P. and T. W. Becker (2010), Where is the real transform
boundary in California?, \emph{Geochem. Geophys. Geosyst} ., 11, Q06012,
doi:10.1029/2010GC003060.

Prats-Iraola, P., Scheiber, R., Marotti, L., Wollstadt, S., \& Reigber,
A. (2012). TOPS interferometry with TerraSAR-X.~\emph{IEEE Transactions
on geoscience and remote sensing} ,~\emph{50} (8), 3179-3188.

Rollins, C., Wright, T.J., Weiss, J.R., Hooper, A.J., Walters, R.J.,
Lazecky, M., Morishita, Y. and Maghsoudi, Y., 2020, December. Tectonic
strain rates across the central Alpine-Himalayan Belt from Sentinel-1
InSAR and GNSS observations, and implications for seismic hazard. In AGU
Fall Meeting 2020. AGU.

Savage, J. C., Gan, W., \& Svarc, J. L. (2001). Strain accumulation and
rotation in the Eastern California Shear Zone.~\emph{Journal of
Geophysical Research: Solid Earth} ,~\emph{106} (B10), 21995-22007.

Sandwell, D.T. and Wessel, P., (2016). Interpolation of 2\selectlanguage{english}-D vector data
using constraints from elasticity. \emph{Geophysical Research Letters} ,
43(20).

Sandwell, D.T., Zeng, Y., Shen, Z.K., Crowell, B., Murray, J.,
McCaffrey, R. and Xu, X., (2016a). The SCEC community geodetic model V1:
Horizontal velocity grid. In SCEC Annu. Meeting.

Sandwell D. T. Xu X. Mellors R. Wei M. Tong X., and Wessel P. (2016b).
GMTSAR: An InSAR Processing System Based on Generic Mapping Tools,
Second Ed., available at
http://topex.ucsd.edu/gmtsar/tar/GMTSAR\_2ND\_TEX.pdf (last accessed
December 2019).

Shen, Z.K. and Liu, Z., 2020. Integration of GPS and InSAR Data for
Resolving 3\selectlanguage{english}-Dimensional Crustal Deformation. \emph{Earth and Space
Science} , 7(4), e2019EA001036.

Shen, Z. K., Wang, M., Zeng, Y., \& Wang, F. (2015). Optimal
interpolation of spatially discretized geodetic data. \emph{Bulletin of
the Seismological Society of America} , 105, 2117-2127.

Smith\selectlanguage{english}-Konter, B., \& Sandwell, D. (2009). Stress evolution of the San
Andreas fault system: Recurrence interval versus locking
depth.~\emph{Geophysical Research Letters} ,~\emph{36} (13).

Tape, C., P. Muse, M. Simons (2009), D. Dong, and F. Webb, Multiscale
estimation of GPS velocity fields, \emph{Geophys. J. Int} ., 179,
945-971.

Tong, X., D. Sandwell, and B. Smith-Konter (2013), High-resolution
interseismic velocity data along the San Andreas Fault from GPS and
InSAR, \emph{J. Geophys. Res. Solid Earth} , 118,
doi:10.1029/2012JB009442.

Tong, X., B. Smith-Konter, and D. T. Sandwell (2014), Is there a
discrepancy between geological and geodetic slip rates along the San
Andreas Fault System?, \emph{J. Geophys. Res. Solid Earth} , 119,
doi:10.1002/2013JB010765.

Tong, X., \& Schmidt, D. (2016). Active movement of the Cascade
landslide complex in Washington from a coherence-based InSAR time series
method.~\emph{Remote Sensing of Environment} ,~\emph{186} , 405-415.

Tymofyeyeva, E., \& Fialko, Y. (2015). Mitigation of atmospheric phase
delays in InSAR data, with application to the Eastern California Shear
Zone. \emph{Journal of Geophysical Research: Solid Earth} , 120(8),
5952-5963.

Tymofyeyeva, E., \& Fialko, Y. (2018). Geodetic evidence for a blind
fault segment at the southern end of the San Jacinto Fault
Zone.\emph{Journal of Geophysical Research: Solid Earth} , 123,
878--891. https://doi.org/10.1002/2017JB014477.

Uieda, L. (2018). Verde: Processing and gridding spatial data using
Green's functions.~\emph{Journal of Open Source Software} ,~\emph{3}
(30), 957.

Wang, K., Xu, X., \& Fialko, Y. (2017). Improving burst alignment in
TOPS interferometry with bivariate enhanced spectral
diversity.~\emph{IEEE Geoscience and Remote Sensing Letters} ,~\emph{14}
(12), 2423-2427.

Wei, M., D. T. Sandwell, and B. Smith-Konter (2010), Optimal combination
of InSAR and GPS for measuring interseismic crustal
deformation,\emph{J.} \emph{Adv. in Space Res} ., doi:
10.1016/j.asr.2010.03.013.

Weiss, J. R., Walters, R. J., Morishita, Y., Wright, T. J., Lazecky, M.,
Wang, H., \ldots{} \& Parsons, B. (2020). High\selectlanguage{english}-resolution surface
velocities and strain for Anatolia from Sentinel\selectlanguage{english}-1 InSAR and GNSS
data.~\emph{Geophysical Research Letters} ,~\emph{47} (17),
e2020GL087376.

Wessel, P., Luis, J. F., Uieda, L., Scharroo, R., Wobbe, F., Smith, W.
H. F., \& Tian, D. (2019). The generic mapping tools version
6.~\emph{Geochemistry, Geophysics, Geosystems} ,~\emph{20} (11),
5556-5564.

Wright, T. J., Parsons, B. E., \& Lu, Z. (2004). Toward mapping surface
deformation in three dimensions using InSAR. \emph{Geophysical Research
Letters} , 31(1).

Xu, X., Sandwell, D. T., Tymofyeyeva, E., Gonzalez-Ortega, A., \& Tong,
X. (2017). Tectonic and anthropogenic deformation at the Cerro Prieto
geothermal step-over revealed by Sentinel-1A InSAR. \emph{IEEE
Transactions on Geoscience and Remote Sensing} , 55(9), 5284-5292.

Xu, X., Ward, L. A., Jiang, J., Smith\selectlanguage{english}-Konter, B., Tymofyeyeva, E.,
Lindsey, E. O., \ldots{} \& Sandwell, D. T. (2018). Surface creep rate
of the southern San Andreas fault modulated by stress perturbations from
nearby large events.~\emph{Geophysical Research Letters} ,~\emph{45}
(19), 10-259.

Zeng Y. and Z.-K. Shen (2016). A Fault-Based Model for Crustal
Deformation, Fault Slip Rates and Off-Fault Strain Rate in
California,\emph{Bulletin of the Seismological Society of America} ,,
106(2), doi:10.1785/0120140250.

Zeng, Y., \& Shen, Z. K. (2017). A Fault\selectlanguage{english}-Based Model for Crustal
Deformation in the Western United States Based on a Combined Inversion
of GPS and Geologic Inputs.~\emph{Bulletin of the Seismological Society
of America} ,~\emph{107} (6), 2597-2612.

Zeng, Y., Petersen, M. D., \& Shen, Z.-K. (2018). Earthquake potential
in California-Nevada implied by correlation of strain rate and
seismicity. \emph{Geophysical Research Letters} , 45, 1778--1785.
https://doi.org/10.1002/ 2017GL075967

\textbf{Hosted file}

\verb`xu_et_al_safs_supps_v7.docx`  available at \url{https://authorea.com/users/529071/articles/597002-integrated-sentinel-1-insar-and-gnss-time-series-along-the-san-andreas-fault-system}

\selectlanguage{english}
\FloatBarrier
\end{document}

