WATER RESOURCES RESEARCH
Supporting Information for Quantifying Ground
Deformation in the Los Angeles and Santa Ana Coastal
Basins due to Groundwater Pumping
Bryan Riel
1
, Mark Simons
2
, Daniel Ponti
3
, Piyush Agram
1
, Romain Jolivet
4
Contents of this file
1. Text S1
2. Figure S2.1
Corresponding author: B. Riel, Jet Propulsion Laboratory, California Institute of Technology,
Pasadena, California, USA. (briel@jpl.nasa.gov)
1
Jet Propulsion Laboratory, Pasadena,
California, USA.
2
Seismological Laboratory, California
Institute of Technology, Pasadena, California,
USA.
3
United States Geological Survey, Menlo
Park, California, USA.
4
École Normale Supérieure de Paris, Paris,
France.
D R A F T March 16, 2018, 10:28am D R A F T
X - 2
RIEL ET AL.: SUPP: LA BASIN GROUND DEFORMATION
3. Figure S2.2
4. Figure S2.3
5. Figure S2.4
6. Table S2.1
S1. Spatiotemporal InSAR Time Series Analysis
In
Riel et al.
[2014], a new method for detecting transient signals in geodetic time series was
proposed that modeled time series as a linear combination of displacement functions chosen
from a
dictionary
(matrix) of functions that resemble secular, seasonal, and transient signals.
The displacement functions for transient signals were third-order time-integrated B-splines
(B
i
-splines) which exhibit one-sided behavior of a particular timescale. For reconstruction of
transient signals of unknown onset times and durations, a highly overcomplete, non-orthogonal
collection of B
i
-splines of various timescales are used to populate the dictionary. Then, transient
detection becomes a least squares minimization of a joint cost function that consists of a data
misfit term and a sparsity-inducing regularization term to limit the total number of dictionary
elements needed to reconstruct the original time series [
Tibshirani
, 1996;
Chen et al.
, 1998;
Riel
et al.
, 2014]. For a geodetic time series
d
∈
R
M
×
1
and a dictionary
G
∈
R
M
×
P
consisting of
P
elements, the resulting cost function would be:
φ
(
m
)
=
‖
Gm
−
d
‖
2
2
+
λ
‖
m
‖
1
,
(1)
where
m
is the vector of coefficients for each element in
G
and the
`
1
-norm regularization term
(scaled by the parameter
λ
) minimizes the number of non-zero values in
m
. Enforcing a sparse
set of dictionary elements is critical for proper decomposition of any time series into its relevant
timescales and enhancing interpretability of the non-zero elements in
m
.
D R A F T March 16, 2018, 10:28am D R A F T
RIEL ET AL.: SUPP: LA BASIN GROUND DEFORMATION
X - 3
Riel et al.
[2014] also introduced a spatial weighting scheme for analyzing a network of GPS
stations through an iterative re-weighting approach that ensures that the estimated
m
i
for a
station
i
is consistent with
m
j
from its neighboring stations. The “consistency” is controlled by
an exponential distance weighting function between station
i
and stations
j
=
1
, . . .,
N
,
j
,
i
in
a network of
N
GPS stations. For InSAR time series analysis, an identical procedure can be
used where each pixel in an interferogram is a “station”, and spatial weighting functions can be
computed for each pixel and the other pixels in the interferogram. However, for an interferogram
with
N
pixels, this approach would require the evaluation of
N
2
weighting functions, which
would be very computationally intensive. An alternative approach would be to formulate a
simultaneous
estimation problem where we estimate
m
i
,
i
=
1
, . . .,
N
all at once in a large least
squares problem. In this approach, the global linear model for
K
interferograms,
{
G
}{
m
}
=
{
d
}
,
would be constructed as:
G
1
0
. . .
0
0 G
1
. . .
0
.
.
.
.
.
.
.
.
.
0
0 0 0 G
1
.
.
.
.
.
.
.
.
.
.
.
.
G
K
0
. . .
0
0 G
K
. . .
0
.
.
.
.
.
.
.
.
.
0
0 0 0 G
K
m
1
m
2
.
.
.
m
N
=
d
1
1
d
2
1
.
.
.
d
N
1
.
.
.
d
1
K
d
2
K
.
.
.
d
N
K
,
(2)
where
G
k
is the
k
-th row of the temporal dictionary corresponding to interferogram
k
,
m
i
is the
coefficient vector corresponding to the
i
-th pixel, and
d
i
k
is the
i
-th pixel of interferogram
k
. The
problem size becomes
{
G
} ∈
R
(
N
·
K
)×(
N
·
P
)
,
{
m
} ∈
R
(
N
·
P
×
1
)
, and
{
d
} ∈
R
(
N
·
K
×
1
)
, where
{
d
}
is a
“flattened” 1-D vector containing the entire InSAR stack.
D R A F T March 16, 2018, 10:28am D R A F T
X - 4
RIEL ET AL.: SUPP: LA BASIN GROUND DEFORMATION
S1.1. Large scale time series analysis using the alternating direction method of multipliers
Applying the cost function in Equation 1 to the linear model in Equation 2 would require
solving a very large convex optimization problem for a length
(
N
·
P
)
solution vector
{
m
}
.
For most modern InSAR time series, this problem cannot successfully be solved on a single
computer. However, very powerful distributed optimization algorithms have been developed in
recent years that take advantage of rapidly improving cluster computing frameworks to minimize
the computational requirement of any one compute node. One such algorithm is the alternating
direction method of multipliers (ADMM) described by
Boyd et al.
[2011]. The most general
form of the problem solved by ADMM is
minimize
f
(
m
)
+
g
(
z
)
subject to
Am
+
Bz
=
c
(3)
with optimization variables
m
and
z
and parameter arrays
A
,
B
, and
c
. Here, the primary variable
we want to solve for,
m
, has been split into two parts,
m
and
z
, such that the convex functions
f
and
g
are
separable
across the splitting [
Boyd et al.
, 2011]. This separability is crucial for
exploiting large scale parallel computing frameworks. The constraint function is then used to
enforce consistency between
m
and
z
. We form the augmented Lagrangian
L
ρ
(
m
,
z
,
y
)
=
f
(
m
)
+
g
(
z
)
+
y
T
(
Am
+
Bz
−
c
)
+
(
ρ
/
2
)‖
Am
+
Bz
−
c
‖
2
2
,
(4)
where
y
is a
dual
variable (from dual ascent optimization methods) and
ρ >
0
is a penalty
parameter. ADMM performs the following iterations:
m
k
+
1
=
argmin
m
L
ρ
(
m
,
z
k
,
y
k
)
(5)
z
k
+
1
=
argmin
z
L
ρ
(
m
k
+
1
,
z
,
y
k
)
(6)
y
k
+
1
=
y
k
+
ρ
(
Am
k
+
1
+
Bz
k
+
1
−
c
)
.
(7)
D R A F T March 16, 2018, 10:28am D R A F T
RIEL ET AL.: SUPP: LA BASIN GROUND DEFORMATION
X - 5
For the sparsity-inducing regularization problem in Equation 1,
f
is the misfit cost function and
g
is the
`
1
-norm regularization term. In this case, we would solve Equation 3 with
A
and
B
set
to the identity and negative identity matrices, respectively, and
c
set to 0.
The original ADMM formulation takes advantage of the separability of the convex functions
f
and
g
by splitting the optimization problem in one of two ways: 1) splitting across the data,
or 2) splitting across the parameters [
Boyd et al.
, 2011]. The former is useful for a modest
number of parameters and a very large number of training data while the latter is useful for the
reverse situation. In Section 1.1.1, we will show how we can combine both splitting approaches
in a consistent way to efficiently analyze InSAR time series. For the current discussion, we
will proceed with the second approach by splitting across the parameters for our linear model
{
G
}{
m
}
=
{
d
}
. For this approach, we partition the variable vector
{
m
}
as
{
m
}
=
(
m
1
, . . .,
m
N
)
where the individual
m
i
’s do not need to be the same sizes. Correspondingly, we partition
{
G
}
=
(
G
1
, . . .,
G
N
)
such that
{
G
}{
m
}
=
Õ
N
i
=
1
G
i
m
i
. This partitioning allows us to formulate
the separable cost function
minimize
‖
N
’
i
=
1
G
i
m
i
−{
d
}‖
2
2
+
λ
N
’
i
=
1
‖
m
i
‖
1
.
The ADMM problem in Equation 3 can then be re-formulated as:
minimize
‖
z
i
−{
d
}‖
2
2
+
λ
N
’
i
=
1
‖
m
i
‖
1
subject to
G
i
m
i
−
z
i
=
0
,
i
=
1
, . . .,
N
,
(8)
where the variables
z
i
are introduced with the same length as the global data vector
{
d
}
. Each
update of Equations 5 - 7 can be done in parallel where a global reduction step (i.e., an
Allreduce
summing operation) is performed prior to the
z
k
+
1
update to combine the individual
G
i
m
i
D R A F T March 16, 2018, 10:28am D R A F T
X - 6
RIEL ET AL.: SUPP: LA BASIN GROUND DEFORMATION
predictions. Each individual worker would also need a copy of the InSAR stack contained in the
global
{
d
}
vector.
Up to this point, the formulation of the ADMM problem has only considered the
`
1
-norm
sparsity-inducing regularization function. However, we still need to enforce spatial coherency
of the
m
vectors for each pixel as is done for GPS networks. One straightforward way to achieve
this is to augment the regularization function
g
with an additional function that penalizes the
spatial structure of the estimated
{
m
}
vector through a covariance matrix,
S
, such that:
g
i
(
m
i
)
=
λ
‖
m
i
‖
1
+
β
‖
S
i
m
i
‖
2
2
,
where
β
is an additional penalty parameter controlling the strength of the spatial penalty. This
joint
`
1
- and
`
2
-norm regularization function is called the
elastic net
problem and has proven to
be very useful for a large number of machine learning and statistics problems [
Zou and Hastie
,
2005;
Jenatton et al.
, 2011]. The original elastic net formulation where
S
=
I
is used for
selecting sparse
groups
of predictors whereas the
`
1
-norm-only regularization tends to select
one predictor per group. This group selection is very useful for our case since we want to choose
a sparse set of B
i
-splines that are spatially grouped together. Appropriate construction of
S
can
then enforce spatial coherency beyond the original elastic net.
In practice, performing ADMM for InSAR time series analysis methods is tractable due to
the sparsity structure of the
{
G
}
array where the vast majority of its elements are zero. We can
use powerful sparse linear algebra libraries to efficiently perform the parallel ADMM iterations.
To maintain tractability, we also wish to construct the spatial covariance matrix
S
with a sparse
number of non-zero elements. If we were to apply the spatial weighting function discussed
for GPS networks, we would have a dense
(
N
·
P
)×(
N
·
P
)
covariance matrix. Instead, since
interferograms are continuous images of ground deformation, we can use a discrete Laplacian
D R A F T March 16, 2018, 10:28am D R A F T
RIEL ET AL.: SUPP: LA BASIN GROUND DEFORMATION
X - 7
smoothing operator to construct
S
such that the spatial distribution of a B
i
-spline of a particular
onset time and duration is smooth. Thus,
S
is very sparse, and we can enhance the inherent
grouping effect of the elastic net with a simple smoothing operation.
S1.1.1. Recursive time series analysis with ADMM
In the previous section, we discussed how we can split the global optimization problem across
parameters in order to use parallel processing computing frameworks. However, each worker
still needs to have a copy of the entire InSAR stack in
{
d
}
. When the number of interferograms
becomes large, or when each interferogram has a large number of pixels, this requirement can
be computationally taxing. We can mitigate this effect by processing the InSAR stack in a
sequential
or
recursive
manner, i.e. one interferogram at a time. This sequential approach
has been discussed in applications that aim to estimate sparse signals in a real-time manner
as data are acquired [e.g.
Angelosante et al.
, 2010;
Kopsinis et al.
, 2011]. In addition to real-
time estimation capabilities, these applications also require low memory usage in deployment
situations where onboard memory is limited. For our purposes, performing time series analysis
on an interferogram-by-interferogram basis using the ADMM method while splitting across
parameters greatly reduces the computational burden on any one compute node.
Following the approach outlined by
Angelosante et al.
[2010], we can re-write the misfit cost
function,
(‖{
G
}{
m
}−{
d
}‖
2
2
)
, as (
{
m
}
T
R
K
{
m
}−
2
{
m
}
T
r
K
)
where
R
K
=
K
’
k
{
G
}
T
k
{
G
}
k
,
r
k
=
K
’
k
{
G
}
T
k
{
d
}
k
,
(9)
and the
k
subscript indicates the subset of the
{
G
}
and
{
d
}
arrays corresponding to interferogram
k
, and
R
K
∈
R
(
N
·
P
)×(
N
·
P
)
and
r
K
∈
R
(
N
·
P
)×
1
. We can then sequentially update
R
K
and
r
K
for a
D R A F T March 16, 2018, 10:28am D R A F T
X - 8
RIEL ET AL.: SUPP: LA BASIN GROUND DEFORMATION
new interferogram by using
R
K
+
1
=
R
K
+
{
G
}
T
K
+
1
{
G
}
K
+
1
,
r
K
+
1
=
r
K
+
{
G
}
T
K
+
1
{
d
}
K
+
1
.
(10)
For a given InSAR time series of
K
interferograms, we can recursively build
R
K
and
r
K
without
needing to load the entire time series into memory and distribute copies of the time series to
each worker. Now, the ADMM minimization problem for the elastic net regularization with
partitioning across parameters becomes:
minimize
[
{
z
}
T
R
K
{
z
}−
2
{
z
}
T
r
K
]
+
N
’
i
=
1
[
λ
‖
m
i
‖
1
+
β
‖
S
i
m
i
‖
2
2
]
subject to
{
m
}−{
z
}
=
0
.
(11)
In the above equation, note that
{
z
}
and
{
m
}
are both vectors of length
(
N
·
P
)
and are variables
in the global optimization problem. By maintaining linear separability in the regularization
function, the
{
m
}
-update using Equation 5 can still be performed in parallel. While the update
of
{
z
}
using Equation 6 is in the global domain,
R
K
is still very sparse and we can make
use of parallel sparse linear algebra packages to efficiently perform the update. Here, we use
the Portable, Extensible Toolkit for Scientific Computation (PETSc) [
Balay et al.
, 1997] in
conjunction with the Multifrontal Massively Parallel Sparse direct Solver (MUMPS) [
Amestoy
et al.
, 2001].
S1.1.2. Selection of penalty parameters
Equation 11 requires selection of two penalty parameters for the
`
1
- and
`
2
-norms,
λ
and
β
,
respectively. One possible approach would be
K
-fold cross validation, which would involve
randomly partitioning the data vector into
K
subsets. Then, for a given
λ
and
β
, we would
solve for
{
m
}
using
K
−
1
subsets for training and use the last subset to test the out of sample
performance of the trained model. We would then iterate over the
K
subsets such that each
D R A F T March 16, 2018, 10:28am D R A F T
RIEL ET AL.: SUPP: LA BASIN GROUND DEFORMATION
X - 9
subset is used as the testing set. The combination of
λ
,
β
that gave the lowest testing area
over the
K
folds would then be the optimal penalty parameters. However, this two-dimensional
exploration of penalty parameters could be quite costly for a large number of data and parameters.
Alternatively, we can reduce the computational cost by using an
independent
source of data as a
testing set and not have to perform a two-dimensional grid search
K
times. In this work, we can
use data from GPS stations within the coastal basins and outlying areas as an independent test of
predictive power of our InSAR time series model. Thus, we perform a single two-dimensional
grid search for
λ
and
β
and choose the combination that results in a time series model that best
agrees with the GPS data.
S1.2. Discussion on Recursive ADMM for InSAR Time Series Analysis
A high fraction of InSAR time series methods and their applications have relied on a pixel-by-
pixel approach where the time history of each pixel is solved independently of its neighboring
pixels. This approach is simple, easily parallelizable, and has proven to be useful for a large
number of ground deformation studies. However, these pixel-by-pixel approaches implicitly
ignore spatially correlated data errors, such as those from atmospheric effects which can be
coherent over length scales of tens to hundreds of kilometers [
Emardson et al.
, 2003;
Lohman
and Simons
, 2005]. Omission of these errors from the time series analysis can lead to significant
changes in the time series model parameters and, perhaps more importantly, their uncertainties.
The formulation of the ADMM minimization problem allows us to directly incorporate data
covariance matrices for these spatially correlated data errors. The recursive update in Equation 10
would be modified by inserting a covariance matrix between the
{
G
}
k
and
{
d
}
k
multiplications.
Accounting for these errors, in conjunction with enforcing spatial coherency in the elastic net
D R A F T March 16, 2018, 10:28am D R A F T
X - 10
RIEL ET AL.: SUPP: LA BASIN GROUND DEFORMATION
regularization scheme, greatly improves the robustness of the time series model to data artifacts
such as unwrapping errors and partial scene coverage due to incomplete SAR frames.
Recursive processing of the interferograms is highly suitable for real-time analysis of InSAR
data. The shorter repeat times of modern SAR platforms has lead to an unprecedented volume
of data being collected to monitor ground motion with high temporal resolution. However,
most modern software packages designed for InSAR time series analysis reprocess the entire
time series each time a new interferogram is acquired. This strategy is not sustainable in
terms of computational time and resources. By maintaining a “state” representative of the
current time series model (encapsulated within the
{
m
}
vector), and by keeping the state
invariant to any specific observation epoch, we can efficiently update our time series model
by recursively updating the ADMM arrays (
R
K
and
r
K
) and rapidly estimating a new state,
{
m
}
. For non-real-time applications, the recursiveness allows for simple parallelization of
interferogram assimilation since each worker can compute their own
R
K
and
r
K
, and a global
reduction operation can combine the arrays before estimation of the time series model.
While the combination of partitioning across the global parameter vector and recursive updat-
ing of the ADMM arrays (
R
k
and
r
K
) greatly increases the computational efficiency of ADMM
for InSAR time series analysis, we are still limited by the requirement that every pixel in an
interferogram requires its own set of temporal dictionary elements. Analysis of a high spatial
resolution InSAR stack with a hundred or more dictionary functions leads to a very large number
of parameters and memory requirements of several tens to over a hundred gigabytes. However,
optimized computing frameworks that specialize in distributed storage and distributed process-
ing, such as Apache Hadoop or Apache Spark, can easily handle such memory requirements
D R A F T March 16, 2018, 10:28am D R A F T
RIEL ET AL.: SUPP: LA BASIN GROUND DEFORMATION
X - 11
by intelligently loading and caching problem data. Thus, future work involves migrating our
ADMM implementation to an appropriate distributed processing framework.
S2. Comparison between Recursive ADMM Time Series Analysis and Other
Decomposition Methods
In many recent InSAR time series studies, time series decomposition techniques based on fea-
ture extraction and dimensionality reduction have been used to isolate geophysical signals from
large geodetic datasets. Two popular techniques are principal component analysis (PCA) and the
closely related independent component analysis (ICA). Both techniques are characterized as
non-
parameteric
where the decomposed components are derived entirely from information obtained
from the data. Both techniques aim to learn inherent correlations between features/variables in
general multivariate data (e.g., singular value decomposition for PCA). PCA has been used on
InSAR time series for hydrology studies ([
Chaussard et al.
, 2014]), modeling volcanic inflation
events ([
Lin et al.
, 2010]), and GPS time series for filtering common mode errors ([
Dong et al.
,
2006]). ICA use for geodesy has been a recent development and has been applied to InSAR time
series for drought analysis ([
Chaussard et al.
, 2017]), GPS time series for detecting mixtures
of transient and steady signals ([
Gualandi et al.
, 2016]), etc. A full review of studies utilizing
PCA/ICA is well outside the scope of this paper, but both techniques have proven to be useful
and efficient for extracting geophysical signals from large time series. In this section, we briefly
assess the performance of ICA for the InSAR time series acquired over the Central and Santa Ana
Coastal Basins and discuss the advantages and disadvantages of ICA compared to our proposed
method.
Prior to the ICA decomposition, we first construct monotonic-in-time time series for each
pixel using a standard Short Baseline Subset (SBAS) approach [
Casu et al.
, 2006]. The time
D R A F T March 16, 2018, 10:28am D R A F T
X - 12
RIEL ET AL.: SUPP: LA BASIN GROUND DEFORMATION
epochs for the time series correspond to the dates of the SAR acquisitions. Then, ICA is
performed on the InSAR time series using the FastICA algorithm ([
Hyvärinen and Oja
, 2000])
to obtain twenty independent components. Figure S2.2 shows the spatial and temporal signature
of three of the “largest” components (i.e., components with time series with the largest dynamic
range). While the components have clear associations with distinct geophysical processes (i.e.,
subsidence in Figure S2.2A and seasonal oscillations in Figure S2.2B), we can also observe
substantial mixing between components. All three components show a mix of seasonal and
long-term signals in both the spatial and temporal data. The third component, which appears
to be associated with the uplift signal around 2005, shows uplift earlier (2003) than when uplift
actually occurred (2005-2006), which suggests some type of trade-off between independent
components. These results indicate that the InSAR data was not properly decomposed into truly
independent processes. While it is not shown here, PCA decompositions gives qualitatively
similar results. In general, geophysical signals that propagate in space (such as the seasonal
oscillations with spatially varying time delays) will tend to be split into multiple components,
which complicates interpretation of single components. Since our primary goal is to cleanly
separate the short- and long-term signals for all data sets (hydraulic head, GPS, and InSAR) in
a consistent manner, we believe that the proposed method is more suitable for our purposes.
As mentioned in the previous paragraph, methods such as ICA or PCA usually perform the de-
composition on pixel time series rather than interferometric observations (line-of-sight displace-
ment between two SAR acquisitions). Thus, a time series estimation approach (such as SBAS)
must be performed prior to decomposition, which may introduce additional noise/uncertainty.
This two-step approach is likely not much of an issue when the signals of interest are large (like
the ones presented here), but it may add difficulty when attempting to detect more subtle signals.
D R A F T March 16, 2018, 10:28am D R A F T
RIEL ET AL.: SUPP: LA BASIN GROUND DEFORMATION
X - 13
Our proposed method works on the interferometric observations directly (by multiplying the
temporal dictionary, G, by an appropriate connectivity matrix as discussed in
Hetland et al.
[2012]). Additionally, by using parameterized time basis functions (e.g., integrated B-splines),
we can easily deal with missing InSAR data by subsetting the matrix G appropriately. Thus, the
B-splines act as both smoothing and interpolating splines where the smoothing is controlled by
the sparsity regularization and the interpolation happens automatically since the B-splines can
be evaluated at any arbitrary time epoch. When using PCA/ICA using traditional approaches,
one often needs to perform some type of inpainting or spatial interpolation of interferograms
prior to decomposition.
Using parameterized time functions in our proposed method, in conjunction with the linear
model formulation, allows for straightforward quantification of uncertainties for reconstructed
time series since we can estimate the covariance of the B-spline amplitudes via the standard
model resolution matrix. Uncertainty quantification for PCA/ICA methods would likely be more
empirical in nature and would involve some form of bootstrapping or randomizing of the data,
which would only provide some measure of uncertainty on component amplitudes and not on
errors in the decomposition (i.e., variable mixing of components for different realizations of the
data).
In summary, since much of the interpretation of the multi-disciplinary data in this manuscript
relies on consistent decomposition into short- and long-term signals, the time series analysis
method proposed in the previous sections is more appropriate for interpreting the hydrologi-
cal processes at different time scales than PCA/ICA methods. Recent improvements in ICA
decomposition have resulted in more “independence” between components ([
Gualandi et al.
,
2016]), which would be promising for similar datasets with spatially-propagating processes and
D R A F T March 16, 2018, 10:28am D R A F T
X - 14
RIEL ET AL.: SUPP: LA BASIN GROUND DEFORMATION
unsteady temporal signatures. We leave further comparison between decomposition methods
for future work.
S3. Supplemental tables and figures
Table S2.1:
Water Replenishment District (WRD) and Orange County Water District (OCWD)
well/port locations, depths, and time delays of peak seasonal hydraulic head.
Well ID
Latitude
Longitude
Port depth (ft)
Delay (days)
AMD-1_1_WB1_MP9
33.855213
-117.856559
-1037.00
120
AMD-2_1_WB1_MP9
33.843828
-117.880102
-1294.00
118
AMD-2_1_WB1_MP8
33.843828
-117.880102
-1154.00
118
AMD-2_1_WB1_MP7
33.843828
-117.880102
-1014.00
115
AMD-2_1_WB1_MP6
33.843828
-117.880102
-824.00
110
AMD-2_1_WB1_MP5
33.843828
-117.880102
-659.00
110
AMD-2_1_WB1_MP4
33.843828
-117.880102
-512.00
114
AMD-2_1_WB1_MP3
33.843828
-117.880102
-387.00
115
AMD-2_1_WB1_MP2
33.843828
-117.880102
-262.00
120
AMD-4_1_WB1_MP9
33.844661
-117.932781
-1124.00
76
AMD-4_1_WB1_MP8
33.844661
-117.932781
-1059.00
76
AMD-4_1_WB1_MP6
33.844661
-117.932781
-794.00
69
AMD-4_1_WB1_MP5
33.844661
-117.932781
-702.00
67
AMD-4_1_WB1_MP4
33.844661
-117.932781
-561.00
68
AMD-4_1_WB1_MP3
33.844661
-117.932781
-381.00
68
AMD-4_1_WB1_MP2
33.844661
-117.932781
-296.00
69
AMD-5_1_WB1_MP5
33.828215
-117.901659
-497.00
99
AMD-5_1_WB1_MP4
33.828215
-117.901659
-415.00
98
AMD-5_1_WB1_MP3
33.828215
-117.901659
-301.00
98
AMD-7_1_WB1_MP12
33.822943
-117.937393
-1169.00
72
AMD-7_1_WB1_MP11
33.822943
-117.937393
-1074.00
71
AMD-7_1_WB1_MP10
33.822943
-117.937393
-934.00
70
AMD-7_1_WB1_MP8
33.822943
-117.937393
-694.00
60
AMD-7_1_WB1_MP7
33.822943
-117.937393
-580.00
59
AMD-7_1_WB1_MP5
33.822943
-117.937393
-371.00
59
AMD-7_1_WB1_MP4
33.822943
-117.937393
-311.00
58
AMD-7_1_WB1_MP3
33.822943
-117.937393
-271.00
57
AMD-7_1_WB1_MP2
33.822943
-117.937393
-221.00
61
AMD-8_1_WB1_MP9
33.828845
-117.989207
-1164.00
59
AMD-8_1_WB1_MP7
33.828845
-117.989207
-859.00
53
AMD-8_1_WB1_MP6
33.828845
-117.989207
-764.00
55
AMD-8_1_WB1_MP4
33.828845
-117.989207
-525.00
50
AMD-8_1_WB1_MP3
33.828845
-117.989207
-315.00
68
BPM-1_1_WB1_MP8
33.855578
-118.026439
-1267.00
67
BPM-1_1_WB1_MP6
33.855578
-118.026439
-890.00
56
D R A F T March 16, 2018, 10:28am D R A F T
RIEL ET AL.: SUPP: LA BASIN GROUND DEFORMATION
X - 15
BPM-1_1_WB1_MP4
33.855578
-118.026439
-613.00
56
BPM-1_1_WB1_MP3
33.855578
-118.026439
-458.00
56
BPM-2_1_WB1_MP13
33.820536
-118.016350
-1931.00
72
BPM-2_1_WB1_MP12
33.820536
-118.016350
-1764.00
72
BPM-2_1_WB1_MP11
33.820536
-118.016350
-1614.00
82
BPM-2_1_WB1_MP9
33.820536
-118.016350
-1367.00
60
BPM-2_1_WB1_MP8
33.820536
-118.016350
-1243.00
59
BPM-2_1_WB1_MP7
33.820536
-118.016350
-1028.00
53
BPM-2_1_WB1_MP6
33.820536
-118.016350
-903.00
53
BPM-2_1_WB1_MP5
33.820536
-118.016350
-778.00
51
BPM-2_1_WB1_MP4
33.820536
-118.016350
-581.00
47
CB-1_1_WB2_MP8
33.840929
-117.952539
-1333.00
75
CB-1_1_WB2_MP6
33.840929
-117.952539
-1053.00
73
CB-1_1_WB2_MP5
33.840929
-117.952539
-873.00
63
CB-1_1_WB2_MP4
33.840929
-117.952539
-663.00
61
CB-1_1_WB2_MP3
33.840929
-117.952539
-443.00
59
CER1_1
33.847258
-118.058747
-1122.62
55
CER1_2
33.847258
-118.058747
-967.62
49
CER1_3
33.847258
-118.058747
-577.62
51
CER2_2
33.877723
-118.104265
-861.62
50
CER2_3
33.877723
-118.104265
-686.62
48
CMP1_1
33.893792
-118.192925
-1311.59
70
CMP1_2
33.893792
-118.192925
-1091.59
71
COSM-1_1_WB1_MP15
33.687495
-117.914829
-1764.00
61
COSM-1_1_WB1_MP14
33.687495
-117.914829
-1599.00
60
COSM-1_1_WB1_MP13
33.687495
-117.914829
-1435.00
50
COSM-1_1_WB1_MP12
33.687495
-117.914829
-1215.00
52
COSM-1_1_WB1_MP11
33.687495
-117.914829
-1103.00
47
COSM-1_1_WB1_MP10
33.687495
-117.914829
-983.00
47
COSM-1_1_WB1_MP9
33.687495
-117.914829
-853.00
48
COSM-1_1_WB1_MP8
33.687495
-117.914829
-723.00
47
COSM-1_1_WB1_MP7
33.687495
-117.914829
-621.00
48
COSM-1_1_WB1_MP6
33.687495
-117.914829
-541.00
47
COSM-1_1_WB1_MP5
33.687495
-117.914829
-451.00
43
COSM-1_1_WB1_MP4
33.687495
-117.914829
-351.00
37
DWN1_1
33.921290
-118.137080
-1079.59
95
DWN1_2
33.921290
-118.137080
-849.59
87
DWN1_3
33.921290
-118.137080
-489.59
58
DWN1_4
33.921290
-118.137080
-279.59
51
FFS-1_1_WB2_MP7
33.864158
-117.897841
-1300.00
94
FFS-1_1_WB1_MP7
33.864158
-117.897841
-1299.00
123
FFS-1_1_WB2_MP6
33.864158
-117.897841
-1160.00
93
FFS-1_1_WB1_MP6
33.864158
-117.897841
-1159.00
126
FFS-1_1_WB2_MP5
33.864158
-117.897841
-1060.00
91
FFS-1_1_WB1_MP5
33.864158
-117.897841
-1059.00
125
FFS-1_1_WB2_MP4
33.864158
-117.897841
-820.00
93
D R A F T March 16, 2018, 10:28am D R A F T
X - 16
RIEL ET AL.: SUPP: LA BASIN GROUND DEFORMATION
FFS-1_1_WB1_MP4
33.864158
-117.897841
-819.00
124
FFS-1_1_WB2_MP3
33.864158
-117.897841
-530.00
98
FFS-1_1_WB1_MP3
33.864158
-117.897841
-529.00
128
FFS-1_1_WB2_MP2
33.864158
-117.897841
-361.00
102
FVM-1_1_WB1_MP15
33.723313
-117.942072
-1323.00
45
FVM-1_1_WB1_MP14
33.723313
-117.942072
-1233.00
46
FVM-1_1_WB2_MP14
33.723313
-117.942072
-1233.00
47
FVM-1_1_WB1_MP11
33.723313
-117.942072
-1003.00
45
FVM-1_1_WB2_MP11
33.723313
-117.942072
-1003.00
41
FVM-1_1_WB1_MP10
33.723313
-117.942072
-896.00
45
FVM-1_1_WB2_MP9
33.723313
-117.942072
-814.00
38
FVM-1_1_WB2_MP8
33.723313
-117.942072
-632.00
33
FVM-1_1_WB2_MP7
33.723313
-117.942072
-560.00
34
FVM-1_1_WB1_MP6
33.723313
-117.942072
-500.00
41
FVM-1_1_WB2_MP6
33.723313
-117.942072
-500.00
34
FVM-1_1_WB1_MP5
33.723313
-117.942072
-450.00
42
FVM-1_1_WB2_MP5
33.723313
-117.942072
-450.00
36
FVM-1_1_WB1_MP4
33.723313
-117.942072
-360.00
43
FVM-1_1_WB2_MP4
33.723313
-117.942072
-360.00
34
GGM-1_1_WB1_MP10
33.777852
-117.948037
-1519.00
62
GGM-1_1_WB1_MP9
33.777852
-117.948037
-1264.00
60
GGM-1_1_WB1_MP8
33.777852
-117.948037
-1074.00
50
GGM-1_1_WB1_MP7
33.777852
-117.948037
-954.00
47
GGM-1_1_WB1_MP5
33.777852
-117.948037
-744.00
46
GGM-1_1_WB1_MP4
33.777852
-117.948037
-552.00
48
GGM-1_1_WB1_MP3
33.777852
-117.948037
-465.00
48
GGM-2_1_WB1_MP10
33.775245
-118.011053
-1629.00
65
GGM-2_1_WB1_MP9
33.775245
-118.011053
-1489.00
64
GGM-2_1_WB1_MP8
33.775245
-118.011053
-1254.00
63
GGM-2_1_WB1_MP7
33.775245
-118.011053
-1149.00
56
GGM-2_1_WB1_MP6
33.775245
-118.011053
-1049.00
52
GGM-2_1_WB1_MP5
33.775245
-118.011053
-954.00
52
GGM-2_1_WB1_MP4
33.775245
-118.011053
-719.00
48
GGM-2_1_WB1_MP3
33.775245
-118.011053
-462.00
47
HBM-1_1_WB1_MP9
33.737435
-117.993264
-1130.00
57
HBM-1_1_WB1_MP8
33.737435
-117.993264
-1038.00
56
HBM-1_1_WB1_MP7
33.737435
-117.993264
-924.00
58
HBM-1_1_WB1_MP6
33.737435
-117.993264
-702.00
47
HBM-1_1_WB1_MP5
33.737435
-117.993264
-562.00
47
HBM-1_1_WB1_MP4
33.737435
-117.993264
-483.00
47
HBM-1_1_WB1_MP2
33.737435
-117.993264
-191.00
49
HBM-2_1_WB1_MP6
33.683092
-117.977290
-447.00
61
HBM-2_1_WB1_MP4
33.683092
-117.977290
-307.00
63
HBM-2_1_WB1_MP3
33.683092
-117.977290
-247.00
65
HBM-6_1_WB1_MP8
33.704608
-118.014292
-578.00
68
HBM-6_1_WB1_MP7
33.704608
-118.014292
-508.00
68
D R A F T March 16, 2018, 10:28am D R A F T
RIEL ET AL.: SUPP: LA BASIN GROUND DEFORMATION
X - 17
HBM-6_1_WB1_MP6
33.704608
-118.014292
-296.00
79
HBM-6_1_WB1_MP5
33.704608
-118.014292
-264.00
79
HBM-6_1_WB1_MP4
33.704608
-118.014292
-215.00
80
IDM-1_1_WB1_MP7
33.714789
-117.769596
-763.00
75
IDM-1_1_WB2_MP7
33.714789
-117.769596
-763.00
58
IDM-1_1_WB1_MP6
33.714789
-117.769596
-703.00
74
IDM-1_1_WB2_MP6
33.714789
-117.769596
-703.00
59
IDM-1_1_WB1_MP5
33.714789
-117.769596
-631.00
73
IDM-2_1_WB1_MP12
33.721587
-117.831767
-1404.00
71
IDM-2_1_WB1_MP11
33.721587
-117.831767
-1259.00
71
IDM-2_1_WB1_MP10
33.721587
-117.831767
-1182.00
71
IDM-2_1_WB1_MP9
33.721587
-117.831767
-1055.00
69
IDM-2_1_WB1_MP8
33.721587
-117.831767
-890.00
71
IDM-2_1_WB1_MP7
33.721587
-117.831767
-713.00
59
IDM-2_1_WB1_MP6
33.721587
-117.831767
-613.00
64
IDM-2_1_WB1_MP5
33.721587
-117.831767
-493.00
57
IDM-2_1_WB1_MP4
33.721587
-117.831767
-353.00
57
LAM-1_1_WB1_MP10
33.795906
-118.067250
-1253.00
59
LAM-1_1_WB1_MP9
33.795906
-118.067250
-1153.00
60
LAM-1_1_WB1_MP8
33.795906
-118.067250
-1073.00
64
LAM-1_1_WB1_MP6
33.795906
-118.067250
-834.00
63
LAM-1_1_WB1_MP5
33.795906
-118.067250
-572.00
53
LAM-1_1_WB1_MP4
33.795906
-118.067250
-472.00
53
LAM-1_1_WB1_MP3
33.795906
-118.067250
-272.00
46
LAM-1_1_WB1_MP2
33.795906
-118.067250
-222.00
46
LB1_3
33.798002
-118.088391
-968.27
54
LB2_1
33.850093
-118.201093
-939.58
53
LB6_4
33.810438
-118.149150
-453.47
45
LB6_5
33.810438
-118.149150
-353.47
43
LMR1_3
33.882950
-118.007100
-616.15
65
LW1_1
33.853585
-118.152247
-948.13
60
NOR1_2
33.911958
-118.055411
-902.62
116
NOR1_3
33.911958
-118.055411
-632.62
113
RH1_4
33.974822
-118.115437
-277.39
93
SAR-1_1_WB2_MP10
33.832691
-117.866790
-1014.00
105
SAR-1_1_WB2_MP9
33.832691
-117.866790
-914.00
100
SAR-1_1_WB2_MP8
33.832691
-117.866790
-894.00
99
SAR-1_1_WB2_MP7
33.832691
-117.866790
-829.00
103
SAR-1_1_WB2_MP6
33.832691
-117.866790
-584.00
111
SAR-1_1_WB2_MP5
33.832691
-117.866790
-519.00
115
SAR-1_1_WB2_MP4
33.832691
-117.866790
-367.00
122
SAR-1_1_WB2_MP3
33.832691
-117.866790
-327.00
123
SAR-1_1_WB2_MP2
33.832691
-117.866790
-297.00
123
SAR-2_1_WB2_MP12
33.818845
-117.871762
-1351.00
98
SAR-2_1_WB2_MP11
33.818845
-117.871762
-1231.00
89
SAR-2_1_WB2_MP10
33.818845
-117.871762
-1101.00
86
D R A F T March 16, 2018, 10:28am D R A F T
X - 18
RIEL ET AL.: SUPP: LA BASIN GROUND DEFORMATION
SAR-2_1_WB1_MP10
33.818845
-117.871762
-1099.00
109
SAR-2_1_WB2_MP9
33.818845
-117.871762
-1021.00
84
SAR-2_1_WB1_MP9
33.818845
-117.871762
-1019.00
106
SAR-2_1_WB2_MP8
33.818845
-117.871762
-981.00
82
SAR-2_1_WB1_MP8
33.818845
-117.871762
-979.00
107
SAR-2_1_WB2_MP7
33.818845
-117.871762
-881.00
75
SAR-2_1_WB1_MP7
33.818845
-117.871762
-879.00
99
SAR-2_1_WB2_MP5
33.818845
-117.871762
-611.00
89
SAR-3_1_WB2_MP11
33.799513
-117.878657
-1393.00
94
SAR-3_1_WB2_MP10
33.799513
-117.878657
-1269.00
82
SAR-3_1_WB2_MP9
33.799513
-117.878657
-1199.00
83
SAR-3_1_WB2_MP8
33.799513
-117.878657
-1074.00
75
SAR-3_1_WB2_MP7
33.799513
-117.878657
-954.00
68
SAR-3_1_WB2_MP5
33.799513
-117.878657
-644.00
87
SAR-3_1_WB2_MP4
33.799513
-117.878657
-514.00
85
SAR-4_1_WB2_MP9
33.785616
-117.883404
-1168.00
58
SAR-4_1_WB2_MP8
33.785616
-117.883404
-1068.00
52
SAR-4_1_WB2_MP6
33.785616
-117.883404
-868.00
65
SAR-4_1_WB2_MP5
33.785616
-117.883404
-738.00
69
SAR-4_1_WB1_MP4
33.785616
-117.883404
-598.00
61
SAR-4_1_WB2_MP4
33.785616
-117.883404
-598.00
57
SAR-4_1_WB1_MP3
33.785616
-117.883404
-478.00
57
SAR-4_1_WB2_MP3
33.785616
-117.883404
-478.00
47
SAR-5_1_WB2_MP10
33.770381
-117.892245
-1543.00
59
SAR-5_1_WB2_MP9
33.770381
-117.892245
-1293.00
47
SAR-9_1_WB1_MP9
33.752549
-117.905783
-1262.00
49
SAR-9_1_WB1_MP7
33.752549
-117.905783
-877.00
41
SAR-9_1_WB1_MP5
33.752549
-117.905783
-606.00
41
SAR-9_1_WB1_MP4
33.752549
-117.905783
-491.00
43
SBM-1_1_WB1_MP7
33.759091
-118.043498
-916.00
61
SBM-1_1_WB1_MP6
33.759091
-118.043498
-706.00
57
SC-4_1_WB1_MP7
33.791031
-117.826774
-660.00
75
SC-4_1_WB1_MP6
33.791031
-117.826774
-573.00
81
SC-5_1_WB1_MP10
33.764610
-117.843246
-1430.00
82
SC-5_1_WB1_MP9
33.764610
-117.843246
-1238.00
75
SC-5_1_WB1_MP8
33.764610
-117.843246
-1024.00
58
SC-5_1_WB1_MP7
33.764610
-117.843246
-937.00
61
SC-5_1_WB1_MP6
33.764610
-117.843246
-807.00
60
SC-5_1_WB1_MP5
33.764610
-117.843246
-670.00
62
SC-6_1_WB1_MP7
33.738024
-117.877766
-1124.00
51
SC-6_1_WB1_MP6
33.738024
-117.877766
-964.00
42
SC-6_1_WB1_MP4
33.738024
-117.877766
-542.00
39
SGA1_1
33.944895
-118.176882
-1346.59
84
SGA1_2
33.944895
-118.176882
-1226.59
77
SGA1_3
33.944895
-118.176882
-816.59
70
SGA1_4
33.944895
-118.176882
-471.59
45
D R A F T March 16, 2018, 10:28am D R A F T
RIEL ET AL.: SUPP: LA BASIN GROUND DEFORMATION
X - 19
WMM-1_1_WB1_MP15
33.752478
-117.972011
-1624.00
53
WMM-1_1_WB1_MP10
33.752478
-117.972011
-1214.00
52
WMM-1_1_WB2_MP9
33.752478
-117.972011
-1065.00
50
WMM-1_1_WB1_MP5
33.752478
-117.972011
-744.00
41
WMM-1_1_WB2_MP3
33.752478
-117.972011
-483.00
44
WMM-1_1_WB2_MP2
33.752478
-117.972011
-361.00
48
D R A F T March 16, 2018, 10:28am D R A F T
X - 20
RIEL ET AL.: SUPP: LA BASIN GROUND DEFORMATION
References
Amestoy, P. R., I. S. Duff, J.-Y. L’Excellent, and J. Koster (2001), A fully asynchronous mul-
tifrontal solver using distributed dynamic scheduling,
SIAM Journal on Matrix Analysis and
Applications
,
23
(1), 15–41.
Angelosante, D., J. A. Bazerque, and G. B. Giannakis (2010), Online adaptive estimation of
sparse signals: Where rls meets the
`
1
-norm,
Signal Processing, IEEE Transactions on
,
58
(7),
3436–3447.
Balay, S., W. D. Gropp, L. C. McInnes, and B. F. Smith (1997), Efficient management of
parallelism in object oriented numerical software libraries, in
Modern Software Tools in
Scientific Computing
, edited by E. Arge, A. M. Bruaset, and H. P. Langtangen, pp. 163–202,
Birkhäuser Press.
Boyd, S., N. Parikh, E. Chu, B. Peleato, and J. Eckstein (2011), Distributed optimization
and statistical learning via the alternating direction method of multipliers,
Foundations and
Trends
®
in Machine Learning
,
3
(1), 1–122.
Casu, F., M. Manzo, and R. Lanari (2006), A quantitative assessment of the SBAS algorithm
performance for surface deformation retrieval from DInSAR data,
Remote Sensing of Envi-
ronment
,
102
(3-4), 195–210.
Chaussard, E., R. Bürgmann, M. Shirzaei, E. Fielding, and B. Baker (2014), Predictability
of hydraulic head changes and characterization of aquifer-system and fault properties from
InSAR-derived ground deformation,
Journal of Geophysical Research: Solid Earth
,
119
(8),
6572–6590.
Chaussard, E., P. Milillo, R. Bürgmann, D. Perissin, E. J. Fielding, and B. Baker (2017),
Remote sensing of ground deformation for monitoring groundwater management practices:
D R A F T March 16, 2018, 10:28am D R A F T
RIEL ET AL.: SUPP: LA BASIN GROUND DEFORMATION
X - 21
Application to the Santa Clara Valley during the 2012–2015 California drought,
Journal of
Geophysical Research: Solid Earth
,
122
(10), 8566–8582.
Chen, S. S., D. L. Donoho, and M. A. Saunders (1998), Atomic decomposition by basis pursuit,
SIAM Journal on Scientific Computing
,
20
(1), 33–61.
Dong, D., P. Fang, Y. Bock, F. Webb, L. Prawirodirdjo, S. Kedar, and P. Jamason (2006),
Spatiotemporal filtering using principal component analysis and Karhunen-Loeve expansion
approaches for regional GPS network analysis,
Journal of Geophysical Research: Solid Earth
,
111
(B3).
Emardson, T., M. Simons, and F. Webb (2003), Neutral atmospheric delay in interferomet-
ric synthetic aperture radar applications: Statistical description and mitigation,
Journal of
Geophysical Research: Solid Earth
,
108
(B5).
Gualandi, A., E. Serpelloni, and M. Belardinelli (2016), Blind source separation problem in
GPS time series,
Journal of Geodesy
,
90
(4), 323–341.
Hetland, E., P. Musé, M. Simons, Y. Lin, P. Agram, and C. DiCaprio (2012), Multiscale InSAR
time series (MInTS) analysis of surface deformation,
Journal of Geophysical Research
,
117
.
Hyvärinen, A., and E. Oja (2000), Independent component analysis: algorithms and applications,
Neural networks
,
13
(4-5), 411–430.
Jenatton, R., J.-Y. Audibert, and F. Bach (2011), Structured variable selection with sparsity-
inducing norms,
The Journal of Machine Learning Research
,
12
, 2777–2824.
Kopsinis, Y., K. Slavakis, and S. Theodoridis (2011), Online sparse system identification and
signal reconstruction using projections onto weighted balls,
Signal Processing, IEEE Trans-
actions on
,
59
(3), 936–952.
D R A F T March 16, 2018, 10:28am D R A F T
X - 22
RIEL ET AL.: SUPP: LA BASIN GROUND DEFORMATION
Lin, Y.-n. N., A. P. Kositsky, and J.-P. Avouac (2010), PCAIM joint inversion of InSAR and
ground-based geodetic time series: Application to monitoring magmatic inflation beneath the
Long Valley Caldera,
Geophysical Research Letters
,
37
(23).
Lohman, R. B., and M. Simons (2005), Some thoughts on the use of insar data to constrain models
of surface deformation: Noise structure and data downsampling,
Geochemistry, Geophysics,
Geosystems
,
6
(1).
Riel, B., M. Simons, P. Agram, and Z. Zhan (2014), Detecting transient signals in geodetic time
series using sparse estimation techniques,
Journal of Geophysical Research: Solid Earth
,
119
(6), 5140–5160.
Tibshirani, R. (1996), Regression shrinkage and selection via the lasso,
J. R. Stat. Soc. B
,
58
(1),
267–288.
Zou, H., and T. Hastie (2005), Regularization and variable selection via the elastic net,
Journal
of the Royal Statistical Society: Series B (Statistical Methodology)
,
67
(2), 301–320.
D R A F T March 16, 2018, 10:28am D R A F T
RIEL ET AL.: SUPP: LA BASIN GROUND DEFORMATION
X - 23
Figure S2.1:
Spatial distribution of largest B
i
-spline coefficients determined from the recursive
ADMM InSAR time series analysis. For each image, the color at each pixel corresponds to the
amplitude of the specified B
i
-spline. The largest estimated coefficient is associated with a rapid
subsidence signal with a time scale of 2.4 years centered in the beginning of 2007 (A). We can
also observe high values for the uplift event due to heavy rainfall in (D) and (E). In (H), we
see that the uplift in the Santa Fe Spring oil field and subsidence in the Pomona basin are well
described by a long-term B
i
-spline centered on 1999.
D R A F T March 16, 2018, 10:28am D R A F T