Quantifying Ground Deformation in the Los Angeles and Santa Ana Coastal Basins Due to Groundwater Withdrawal - Riel_et_al-2018-Water_Resources

RESEARCH ARTICLE

10.1029/2017WR021978

Quantifying Ground Deformation in the Los Angeles and

Santa Ana Coastal Basins Due to Groundwater Withdrawal

Bryan Riel

1,2

, Mark Simons

, Daniel Ponti

, Piyush Agram

, and Romain Jolivet

Jet Propulsion Laboratory, Pasadena, CA, USA,

Seismological Laboratory, California Institute of Technology, Pasadena,

CA, USA,

United States Geological Survey, Menlo Park, CA, USA,

Ecole Normale Sup

erieure de Paris, Paris, France

Abstract

We investigate complex surface deformation within the Los Angeles and Santa Ana Coastal

Basins due to groundwater withdrawal and subsequent aquifer compaction/expansion. We analyze an 18

year interferometric synthetic aperture radar (InSAR) time series of 881 interferograms in conjunction with

global positioning system (GPS) data within the groundwater basins. The large data set required the devel-

opment of a distributed time series analysis framework able to automatically decompose both the InSAR

and GPS time series into short-term and long-term signals. We find that short-term, seasonal oscillations of

ground elevations due to annual groundwater withdrawal and recharge are unsteady due to changes in

seasonal withdrawal by major water districts. The spatial pattern of seasonal ground deformation near the

center of the basin corresponds to a diffusion process with peak deformation occurring at locations with

highest groundwater production. Long-term signals occur over broader areas and are ultimately caused by

long-term changes in groundwater production. Comparison of the geodetic data with hydraulic head data

from major water districts suggests that different regions of the groundwater system are responsible for

different temporal components in the observed ground deformation. Short-term, seasonal ground deforma-

tion is caused by compaction of shallower aquifers used for the majority of groundwater production

whereas long-term ground deformation is correlated with delayed compaction of deeper aquifers and

potential compressible clay layers. These results demonstrate the potential for geodetic analysis to be an

important tool for groundwater management to maintain sustainable pumping practices.

1. Introduction

In regions over and adjacent to active aquifer systems, ground deformation can occur as a result of ground-

water withdrawal, long-term drought effects, heavy rainfall, and artificial recharge of the aquifers (e.g., Gallo-

way & Burbey, 2011; Todd & Mays, 1980). Ground deformation is a response to changes in pore pressure in

the aquifers, which changes the effective stress on the aquifer system’s granular matrix and causes contrac-

tion or expansion of the pore spaces. From a natural hazard perspective, land subsidence following ground-

water withdrawal and compaction of aquifer systems is of particular interest when the subsidence is long

term, leading to increased strain on infrastructure, potential formation of earth fissures and surface faults,

and changes to surface water drainage (Galloway & Burbey, 2011). Long-term subsidence can occur as a

result of a slow decline in groundwater levels or permanent compaction due to effective stress levels

exceeding preconsolidation stress levels (Wilson & Gorelick, 1996). More generally, ground deformation due

to extraction of fluids in subsurface reservoirs can be described by short-term, elastic responses and long-

term, inelastic compression or poroelastic rebound processes (Galloway & Burbey, 2011; Todd & Mays,

1980). Many municipal water districts closely monitor hydraulic head data to track pore pressures in order

to maintain sustainable pumping practices and prevent effective stress levels within aquifers from exceed-

ing their preconsolidation levels.

Measurements of ground deformation complement hydraulic head data for groundwater monitoring. Con-

tinuous monitoring of ground elevations at discrete points can be achieved through leveling and global

positioning system (GPS) data. Additionally, GPS data provides measurements of horizontal motions, which

can be useful for quantifying and modeling aquifer properties (Bawden et al., 2001; Galloway & Burbey,

2011). Recently, interferometric synthetic aperture radar (InSAR) has proven to be a very useful remote sens-

ing technique for acquiring spatially dense ground deformation measurements for deformation driven by

hydrologic and geothermal fluid processes. Interferograms have been used to observe subsidence in urban

Key Points:

Groundwater pumping induces

ground deformation on multiple

spatial and temporal scales

Different aquifer regions display

distinct modes of deformation

InSAR time series analysis reveals fine

scale geologic structures and

localized subsidence features

Supporting Information:

Supporting Information S1

Correspondence to:

B. Riel,

briel@jpl.nasa.gov

Citation:

Riel, B., Simons, M., Ponti, D., Agram, P.,

& Jolivet, R. (2018). Quantifying ground

deformation in the Los Angeles and

Santa Ana Coastal Basins due to

groundwater withdrawal.

Water

Resources Research

, 3557–3582.

https://doi.org/10.1029/

2017WR021978

Received 13 OCT 2017

Accepted 17 APR 2018

Accepted article online 20 APR 2018

Published online 18 MAY 2018

V

2018. American Geophysical Union.

RIEL ET AL.

3557

Water Resources Research

areas due to groundwater withdrawal (Bawden et al., 2001; Buckley et al., 2003), quantify spatial variations

of geologic structures in irrigated areas (Valentine et al., 2001), and assess sustainability of groundwater

withdrawal and reinjection practices in geothermal plants (Massonnet et al., 1997). From an inverse model-

ing perspective, the spatial density of ground observations from InSAR data has also allowed for estimation

of aquifer system storage parameters (Hoffmann et al., 2003). Studies have shown that InSAR data acquired

over sufficient timespans can be used to quantify seasonal deformation caused by the annual cycle of

groundwater pumping and recharge, as well as longer-term subsidence from accumulated overdraft of

aquifers (Amelung et al., 1999; Bawden et al., 2001; Galloway et al., 1998; Galloway & Hoffmann, 2007).

Another class of InSAR techniques uses ‘‘stacks’’ of many coregistered interferograms collected over a finite

time period to construct a full time series of ground deformation. These techniques have been used to mea-

sure the evolution of land subsidence in the Santa Clara Valley in California (Chaussard et al., 2014; Schmidt

€

urgmann, 2003), subsidence and uplift in Phoenix, Arizona (Miller & Shirzaei, 2015), seasonal uplift and

subsidence in the Los Angeles area (Lanari et al., 2004; Watson et al., 2002), etc.

Ground deformation within and around the Los Angeles area has been measured for several decades using

GPS data from the Southern California Integrated GPS Network (SCIGN) and InSAR data. Historically, the pri-

mary driver for acquiring geodetic data was to quantify the rate of tectonic contraction across the region

and the rate of elastic loading on potentially seismogenic faults such as those involved in the 1987 Whittier

Narrows and 1994 Northridge earthquakes (e.g., Watson et al., 2002). However, many of the geodetic signals

used to study these fault systems are contaminated or completely obscured by nontectonic processes such

as groundwater pumping and oil extraction (Bawden et al., 2001). Several geodetic studies have thus aimed

to quantify the total contribution of nontectonic sources of deformation for the Los Angeles area.

Bawden et al. (2001) used a series of interferograms from 1997 to 1999 to observe several anthropogenic

deformation processes in the basins surrounding the Los Angeles area, including seasonal uplift and subsi-

dence due to groundwater withdrawal in the Santa Ana Coastal Basin which is the primary source of

groundwater for Orange County. The larger magnitude of the subsidence signal as compared to the uplift

(60 mm for the former, 50 mm for the latter) implied a net subsidence signal in the basin thought to be due

to inelastic compaction of lower permeability aquitards within the aquifer system. Watson et al. (2002),

Lanari et al. (2004), and Zhang et al. (2012) extended the analysis by Bawden et al. (2001) by including more

interferograms over a longer time span. In particular, Lanari et al. (2004) applied the small baseline subset

(SBAS) algorithm (Berardino et al., 2002) to produce a time series model from 1995 to 2002. Cross correla-

tion of the spatially varying time series with a reference sinusoid was performed to compute a time shift for

each ground point. This time shift map revealed sharp boundaries for the region of the basin responding to

the annual groundwater withdrawal and recharging, and heterogeneous time shifts within the basin also

suggested lateral variability in hydraulic conductivity. The spatially dense measurements of ground defor-

mation provided by InSAR time series can therefore elucidate subtle characteristics of time-dependent

ground deformation within the basin due to changes in aquifer pressure caused by anthropogenic and nat-

ural processes. In these studies, seasonal oscillations were generally quantified in an average sense under

the assumption that the oscillations were purely sinusoidal with a period of 1 year.

In this study, we first explore the geologic and hydrologic setting of the coastal basins in the Los Angeles

area. Following a brief discussion on how groundwater-level changes drive ground deformation in section

2.2, we examine the time history of groundwater levels for the main aquifer systems in section 3.1 using

hydraulic head data from the Water Replenishment District (WRD) in Los Angeles and the Orange County

Water District (OCWD). We show time series that exhibit highly complex time histories with time-varying

seasonal variations and various longer-term trends related to background groundwater levels. We then per-

form an initial comparison between groundwater levels and surface deformation using GPS data to show

that deformation signals with short-term variations are driven by short-term variations in the principal aqui-

fer system while long-term deformation signals appear to be driven by long-term variations in the deeper

portions of the aquifer system.

In section 4, we perform an InSAR time series analysis on an expanded data set that includes interferograms

from 1992 to 2011 to investigate the spatial variations in short-term, seasonal ground deformation. By using

the assumption that seasonal ground motions can be represented by a linear combination of sinusoids, we

can compare the timing and amplitude of peak seasonal deformation of ground points within the basin to

study the effects of aquifer geometry, groundwater pumping practices, and hydraulic conductivity on the

Water Resources Research

10.1029/2017WR021978

RIEL ET AL.

3558

ground response. This method of time series analysis is similar to the ones used in Bawden et al. (2001) and

Watson et al. (2002). However, since both hydraulic head and ground deformation time series show nonsi-

nusoidal seasonal variations, in section 5, we develop a new method for InSAR time series analysis that

automatically decomposes the time series into generic long-term and short-term signals. We compare the

decomposed InSAR time series to all available monitoring wells to better understand the relationship

between deformation signals of different durations and different regions of the groundwater system.

Finally, in section 6, we discuss how these results are related to groundwater management with regards to

sustainable pumping practices and we relate several of our InSAR observations to physical processes based

on groundwater dynamics.

2. Background

2.1. Hydrogeology and Structure of the Los Angeles Basins

The Los Angeles area consists of several basins containing groundwater systems, including the Los Angeles

Central and West Coast Basins and the Santa Ana Coastal Basin in Orange County (Figure 1). The Central

and West Coast Basins lie within Los Angeles County and are monitored by WRD whereas the Santa Ana

Coastal Basin lies within Orange County and is monitored by OCWD. The Los Angeles Central Basin (hereaf-

ter referred to as the Central Basin) is separated from the West Coast Basin by the northern portion of the

Newport-Inglewood Fault Zone (NIF), which acts locally as a barrier for fluid flow between the two basins

(Thiros et al., 2010). In this study, we focus on ground deformation within the Central and Santa Ana Coastal

Basins. While this region is tectonically active (about 4 mm/yr of uniaxial contraction along thrust faults;

Bawden et al., 2001), motion along the major faults bounding the basins is relatively minimal compared to

the seasonal motion caused by groundwater fluctuations (less than 1 mm/yr for NIF; Southern California

Earthquake Data Center, 2012).

Both the Central and Santa Ana Coastal Basins have groundwater systems with similar structural characteris-

tics. The forebay area occupies about 38 percent of the Santa Ana Coastal Basin and occupies the eastern

region that meets the Santa Ana Mountains (Figure 1). Here unconfined vertical movement of water is not

restricted by laterally extensive clay layers (Woodside & Westropp, 2015). Groundwater recharge occurs

mainly in the forebay area. The confined area (also known as the

‘‘pressure area’’) is considerably larger and extends from the western

edge of the forebay area to the Pacific Ocean (Thiros et al., 2010).

Here there are laterally continuous thick layers of silt and clay that

restrict vertical flow of groundwater, causing aquifers to be under con-

fining pressure. While the actual groundwater system consists of sev-

eral contiguous aquifer units, clay interbeds, and confining aquitards,

the OCWD has developed a simplified model consisting of shallow,

middle (principal), and deep aquifer systems to approximate ground-

water flow (Figure 2; Ehman et al., 2014; Woodside & Westropp, 2015).

The shallow aquifer system generally spans depths of up to 60 m for

most of the basin and is unconfined in the forebay region. The princi-

pal aquifer system, which supports over 90% of the basin pumping in

Orange Country, is generally greater than 300 m thick for much of the

basin. The deep aquifer system, which is limited in production capabil-

ity due to the presence of amber colored groundwater (which

requires extra treatment to remove colors and odors from the water),

defines water-storing units up to 600 m in depth in the center of the

Santa Ana Coastal Basin (Woodside & Westropp, 2015).

2.2. Aquifer Compaction Theory

The relationship between ground deformation and aquifer pressure

for

confined

aquifers can be explained using groundwater flow theory

based on the Principle of Effective Stress (Terzaghi, 1923). Effective

stress can be expressed as

Figure 1.

Location and tectonic setting of the Los Angeles Central and West

Coast Basins and the Santa Ana Coastal Basin. The thick gold lines represent

major faults in the area, including the Newport-Inglewood Fault (NIF), the Whit-

tier Fault (WF), the Palos Verdes Fault (PVF), and the Hollywood Fault (HF). The

thin black lines indicate county boundaries, and the thick blue line corresponds

to the Santa Ana River. The Prado Dam on the east side of the Santa Ana Moun-

tains is the primary flood control facility for the downstream Santa Ana River.

The dashed red line indicates the approximate boundary between the forebay

and confined areas of the groundwater system modified from estimates by the

WRD and OCWD. The inset shows the location of the study area along the

California coast.

Water Resources Research

10.1029/2017WR021978

RIEL ET AL.

3559

5

2

;

(1)

where

and

are components of the effective and total stress tensors, respectively,

is the pore fluid

pressure, and

is the Kronecker delta where

5

;

5

;

(

. Assuming a Newtonian fluid and strains pri-

marily for the

component (due to the natural horizontal orientation of aquifer units), we can simplify

equation (1) to

5

2

(2)

By also assuming that changes in the total/overburden stress are negligible for our study period (since

hydraulic head changes do not change the state of the confining units), changes in effective stress can then

be simply expressed as

D

52

D

(3)

Note that the assumption of constant overburden stress does not apply to unconfined aquifers since by def-

inition they are not under confining stress. It can be shown that changes in stress can be related to changes

in hydraulic head as (Heath, 1982)

D

52

D

;

(4)

where

is the hydraulic head,

is the density of water, and

is the gravitational acceleration constant. In

order to relate pore fluid pressure changes to ground deformation, we use the definition for one-

dimensional skeletal compressibility,

, as the ratio of vertical strain to vertical effective stress:

5

2

D

D

;

(5)

where the

D

is the change in thickness of a control volume with initial thickness

(Galloway & Burbey, 2011).

This equation thus relates the compaction and expansion of sediments to changes in effective stress; integrated

compaction over the entire depth of the aquifer system is equivalent to what we measure as land subsidence

(Chen et al., 2016). The skeletal specific storage,

5

, can be used to combine equations (4) and (5) as

5

5

D

D

;

(6)

Figure 2.

Stratigraphic framework of the Los Angeles Basin (modified from Sheet 19 in Ehman et al., 2014). Aquifer structure is interpreted from connected seism

lines within the basin (red line in map inset). Boundaries of stratigraphic sequences are indicated by solid colored lines, and the locations and dept

hs of several

United States Geological Survey (USGS) monitoring wells projected onto the seismic line are also shown. The purple shaded area indicates the approxi

mate depth

extent of the ‘‘principal’’ aquifer system, i.e., the depth range where most groundwater withdrawal occurs. The orange shaded area indicates the dee

p aquifer

system while the unshaded sequences are associated with the shallow aquifer system. One second in two-way travel time is approximately equal to 1.1 km

Water Resources Research

10.1029/2017WR021978

RIEL ET AL.

3560

where

is the skeletal storage coefficient. The storage coefficient (storativity),

, is the sum of

and a com-

ponent of the aquifer system storage attributed to the compressibility of water; therefore, when water is

assumed to be incompressible,

5

. This same relationship holds for the specific storage,

, and the skele-

tal specific storage,

. In unconfined aquifers, the storativity is also controlled by the specific yield,

, such

that

5

1

. However, in unconfined aquifers, the contribution from the specific yield is usually much

larger than that due to compression of the aquifer.

Consolidation experiments on typical sediments have shown that specific storage for a given material can

behave very differently depending on whether the effective stress is above or below the preconsolidation

stress (i.e., the previous maximum effective stress). This boundary essentially separates the regimes of elastic

and inelastic deformation depending on the head level. To account for these two regimes, two separate

skeletal specific storages are used:

5

for

pre

for

pre

;

(

(7)

where

and

are the elastic and inelastic skeletal specific storage coefficients, respectively. In the

inelastic regime when head drops below the preconsolidation level (the previous level of minimum hydrau-

lic head), irreversible loss of water storage occurs, and in the case where this state persists, the skeletal spe-

cific storage is expected to vary proportionally to the logarithm of the effective stress (Galloway & Burbey,

2011). From a water management perspective, the preconsolidation level is used to define sustainable

pumping rates in order for aquifer deformation to remain elastic.

3. Groundwater Level and Ground Deformation Time History

3.1. Monitoring Groundwater Levels With Hydraulic Head Time Series

Since the primary driver of ground deformation within the coastal basins is changes in aquifer pressure (sec-

tion 2.2), we gather data from several multiport monitoring wells located within the coastal basins that mea-

sure hydraulic head levels at various aquifer depths. Here we use twelve monitoring wells operated by WRD

and thirty wells operated by OCWD (Figure 3). For both sets of wells, the depths of the sample ports allow

us to determine hydraulic head levels for distinct aquifer units. For unconfined aquifers, the hydraulic head

is the elevation of the water table where the hydrostatic pressure of the water is equal to the atmospheric

pressure. For confined aquifers, the hydraulic head level corresponds to the piezometric (or potentiometric)

surface, which is the elevation water would rise to in a piezometric well (Todd & Mays, 1980). Groundwater

levels as measured by hydraulic head time series for OCWD well sta-

tion SAR-9 (located in the middle of the Santa Ana Coastal Basin

encompassing the confined aquifers) show large annual oscillations,

particularly at the depth of the principal aquifer system (Figure 4).

These oscillations are caused by the annual cycle of groundwater

recharge (artificial and natural) and withdrawal. From the early 1990s

to the mid-2000s, groundwater levels typically peaked around March

after the rainy season and are at their lowest toward the end of the

summer during the heaviest periods of groundwater withdrawal. This

time period also corresponds to the largest annual fluctuations in

groundwater levels due to water storage programs that encouraged

increased groundwater withdrawal during the summer months when

regional demand for imported water is higher (Woodside & Westropp,

2015). Superimposed over the annual water-level fluctuations are

longer-term variations. Generally, we can observe a decrease in overall

water levels, which started in 1970 after the basin was essentially

refilled after replenishment from Colorado River water (Woodside &

Westropp, 2015). We can also observe two transient increases in

groundwater levels in 2005 and in 2012, both of which followed peri-

ods of decreased pumping activity due to increased water flow (base

flow from the upstream watershed plus storm runoff) into the Santa

Figure 3.

Distribution of continuous GPS and well data used in this study. The

white circles correspond to WRD wells used in this study while the white dia-

monds correspond to OCWD wells. The black triangles show the GPS coverage

in this area. Hydraulic head time series for OCWD well station SAR-9 is shown in

Figure 4. GPS data were acquired from the Scripps Orbit and Permanent Array

Center (SOPAC; http://sopac.ucsd.edu/). All GPS data used in this study have

been postprocessed for outlier removal and common mode filtering.

Water Resources Research

10.1029/2017WR021978

RIEL ET AL.

3561

Ana River. Note that during other periods of increased water flow into

the Santa Ana River (e.g., around 1995 and 1998), groundwater levels

did not see a corresponding increase because groundwater produc-

tion was not adjusted in response to the increased inflow (see section

6.1 for further discussion).

3.2. Groundwater Diffusion and Varying Response Times

Since groundwater is pumped at discrete points within the Central

and Santa Ana Coastal Basins, we expect spatial variations in hydraulic

head changes in response to pressure gradients. Low permeabilities

in clay interbeds and confining units can cause time delays in hydrau-

lic head levels from one aquifer unit to another, as well as intraaquifer

delays. Additionally, variations in aquifer thickness can affect the prop-

agation of hydraulic head changes.

Three-dimensional flow of ground water in porous media can be

described by the differential equation for a control volume (e.g.,

Jacob, 1950; Keranen et al., 2014):

1

1

52

1

;

(8)

where

, and

are the rectangular components of the instantaneous bulk fluid velocity in the control

volume, and

;

is a source term that can vary in time and in space. For the purpose of this discus-

sion, we assume one-dimensional flow in the

-direction. From equation (8), it can be shown that propaga-

tion of hydraulic head in an aquifer can be written as (Hantush, 1962)

1

1

5

2

;

(9)

where

5

is the hydraulic conductivity and

5

is the aquifer thickness. In the special case of flow

through an aquifer with uniform thickness and conductivity, equation (9) can be written as (Fetter, 2000)

5

2

;

(10)

For the following discussion, let us assume a point source,

;

5

2

, and assume

(

) is periodic

to represent groundwater recharge and withdrawal. Therefore, for a time-varying source with a temporal

frequency

5

cos

, the solution to equation (10) will have the approximate form (Guenther & Lee,

1996):

;

2

ffiffiffiffiffi

2

cos

ffiffiffiffiffiffi

2

6

(11)

In the above equation, there is a frequency dependent attenuation term that would damp out the diffusion

of the hydraulic head at higher source frequencies. In other words, head variations due to a point source

perturbation in water pressure (such as a production well) would decay more rapidly away from the source

for a faster withdrawal/pumping cycle. The

factor in the periodic term in the above equation controls

the rate of diffusion, or hydraulic response time, of the pressure perturbation. For larger

, we would expect

a shorter response time. By solving for the homogeneous solution to equation (10), we can also estimate

the material-dependent response time for confined aquifers as (Alley et al., 2002)

5

;

(12)

where

is the response time and

is a characteristic thickness for a specific aquifer unit. By defining the

hydraulic diffusivity as the ratio

, it can be seen that aquifers with a higher diffusivity would experience

shorter response times.

To assess timing differences in the vertical dimension, we use the SAR-9 hydraulic head data to estimate

the amplitude and time to peak signal for the annual oscillations in head. Since the head data contain both

seasonal variations and long-term trends, we model the time series as a linear combination of sinusoids

Figure 4.

Hydraulic head time series for OCWD well SAR-9 for selected ports in

the shallow, principal, and deep aquifer systems. The dots represent the raw

head data while the solid lines are interpolated data using smoothing splines.

The green bars represent the total water flow into the Santa Ana River (base

plus storm flow) measured at Prado Dam (Woodside & Westropp, 2015). The

largest head variations occur in the principal aquifer system, where the

majority of groundwater withdrawal occurs.

Water Resources Research

10.1029/2017WR021978

RIEL ET AL.

3562

with annual and semiannual periods for the seasonal signal and third-order integrated B-splines (hereafter

referred to as B

-splines; see Hetland et al., 2012) for the transient, long-term trends:

5

5

cos

1

sin

1

5

2

;

(13)

where

5

5 years and

5

1 year, and

2

represent the

-splines centered at time

(Hetland et al.,

2012). Here we partition the data time span into 32 evenly spaced knot times

so that the

-splines each

have an effective duration of 2

2

2

. We estimate the coefficients

, and

simultaneously using

regularized least squares for the cost function:

5

argmin

2

1

2

;

(14)

where

is the time series data,

is the temporal design matrix containing the sinusoids and B

-splines

along the columns,

is the vector of coefficients of the elements in

, and

is a prior covariance matrix

for regularization. For this analysis, we set

to be the identity matrix with an optimal scaling coefficient

selected using

-fold cross validation. After estimating

for each well, we can compute the amplitude and

phase delay for each seasonal component as

5

ffiffiffiffiffiffiffiffiffiffiffiffiffi

1

;

(15)

5

tan

2

;

(16)

where

is the amplitude and

is the phase delay, or time to peak signal. We repeat this procedure for the

time series at each port depth to estimate the time to reach peak head levels for the seasonal signal between

1996 and 2006 (supporting information Table S2.1) and for the transient increase starting in 2004–2005

Figure 5.

Depth profile of (left) amplitude and (right) phase delay of SAR-9 hydraulic head data. The blue lines correspond to

the average seasonal oscillations while the red lines correspond to the water-level recovery signal initiating in 2004/2005 due

to decreased pumping activity. The horizontal shaded areas centered on the blue and red lines indicate the uncertainties asso-

ciated with those values. The vertical shaded areas represent the depths corresponding to the three aquifer system layers as

estimated by the OCWD three-layer model, and the brown hatched regions indicate the approximate range of depths of the

aquitards separating the aquifer system layers. The amplitudes for both the seasonal oscillation and transient increase are

largest in (left) the principal zone; correspondingly, the phase delays are shorter in (right) the principal zone.

Water Resources Research

10.1029/2017WR021978

RIEL ET AL.

3563

caused by a decrease in overall groundwater withdrawal (Figure 5). For both the annual and transient peak

signals, we observe that water levels in the principal aquifer system reach peak levels earlier than both the

shallow and deep aquifer systems. The peak-to-peak seasonal increases in water levels are much higher in the

principal aquifer system, which is expected since 90% of groundwater withdrawal is from the confined princi-

pal aquifer system.

The depth-dependent timing for the annual signal and 2004–2005 transient signal are noticeably different.

It appears to take longer for hydraulic head levels to equilibrate for the 2005 increase than for the annual

cycle. These observations are consistent with the solution to 1-D diffusion in equation (11). The transient

increase in groundwater levels following 2005 has an effective period of

4 years while the seasonal cycle

of recharge and withdrawal has a period of

1 year. The former has a slower rate of diffusion but lower

attenuation of head amplitudes within the aquifer system whereas the latter has a higher diffusion speed

but much more rapid attenuation of amplitude with depth. We conclude that the timing differences

between these two processes are entirely due to variations in the temporal character of the source.

Propagation of hydraulic head is described by a three-dimensional diffusion process. We thus investigate

timing differences in the horizontal direction by estimating the average time of peak hydraulic head for all

wells in the principal aquifer system using the same modeling procedure in equation (13). Because the

OCWD model for the depths of the shallow, principal, and deep aquifer systems is only an approximation of

the actual basin structure, we classify monitoring wells as being located within the ‘‘principal’’ aquifer sys-

tem primarily by the characteristics of their hydraulic head time series. Specifically, for each well, we choose

ports with the largest seasonal variations in head levels (typically more than 2 times the variations of the

shallowest ports) and with similar times of peak hydraulic head (1–2 months differences in peak head times)

since diffusion should be relatively fast in the principal aquifer system. We then compute the median peak

time of those ports to represent the peak time for that well.

Monitoring wells in the center of the basin have earlier peak times than wells in the margins of the basin

(Figure 6). The principal aquifer system generally reaches peak groundwater levels between February and

March in the central region of the basin while the margins peak between March and April. Since the princi-

pal aquifer system is thickest and most productive in the center of the basin, the majority of groundwater

production is obtained from units in the basin center. Thus, the earlier peak times for the central monitoring

wells are mostly controlled by the timing of nearby pumping activity. Within about 10–15 km from the

basin center, and particularly in the Santa Ana Coastal Basin, we can observe a gradient in peak times from

Figure 6.

Time to peak seasonal signal for hydraulic head data from well ports located in the principal aquifer system.

Large-scale production wells (wells with seasonal withdrawal volume amplitudes greater than 3

3

[300 ac ft]) are

indicated by black diamonds. Wells in the center of the basin closer to areas of concentrated pumping reach their peak

signals earlier in the year compared to wells in the margins of the basin. Inset shows timing of wells projected onto the

A–A

line. The dashed white line indicates the approximate boundary between the forebay and confined areas.

Water Resources Research

10.1029/2017WR021978

RIEL ET AL.

3564

the center toward the margins, which suggests horizontal propagation of pressure away from the main

pumping sites. Comparison of the propagation speeds of groundwater in the vertical and horizontal direc-

tions in the center of the basin for the annual cycle reveals that propagation speeds are approximately 1

order of magnitude slower in the vertical direction than the horizontal direction (

30 m/d versus

300 m/

d, respectively). This discrepancy is consistent with the general observation that vertical conductivity is

much lower than horizontal conductivity due to the horizontal orientation of bedding and the likely pres-

ence of laterally continuous clay beds within the main aquifer units that impede vertical groundwater flow.

This propagation signal becomes a little less clear in the Central Basin in Los Angeles County, which is likely

affected by the existence of a shelf that extends east from the NIF limiting the interconnectivity of aquifer

units.

The later peak times in the eastern margin of the basin (closer to the recharge zone) are later than expected

for normal horizontal diffusion. Groundwater flow in this area is substantially different than in the center of

the basin due to the unconfined and semiconfined nature of aquifer units that are pumped in this area.

Due to the bowl-shaped structure of the basin and the overall thinning of aquifers from the center of the

basin toward the margins, groundwater withdrawal in the center of the basin occurs primarily from

younger

shallower units while withdrawal in the margins occurs primarily from

older

, deeper units that are likely

folded and deformed (see Figure 2). These stratigraphic differences will also be associated with large varia-

tions in storage coefficient between confined and unconfined aquifer units. Therefore, the later delay times

for the wells in the eastern margin may be due to the effect of greater vertical propagation of groundwater

from older units to younger units and large lateral variations in storage coefficient. Other factors could

include lateral variations in hydraulic conductivity or the presence of faults acting as barriers to fluid flow.

3.3. Ground Deformation History

From equation (6), we know that ground deformation over a confined aquifer is approximately proportional

to changes in hydraulic head under certain simplifying assumptions, mainly that aquifer or aquitard com-

pression and expansion is elastic for effective stress levels less than the previous maximum effective stress.

The proportionality constant in the elastic relationship is the skeletal storage coefficient,

. However,

ground deformation can also be affected by inelastic compaction of clays in aquitards and clay-rich inter-

beds. In this case, ground deformation will not be correlated with hydraulic head and is expected to vary

exponentially with time (Chaussard et al., 2014). The exponential relationship is derived from the theory of

hydrodynamic consolidation and is used to describe the delayed response of fine-grained materials after

effective stress levels have surpassed previous maximum levels.

We start with the hypothesis that all deformation within the basin is elastic. For short-term signals driven by

the annual cycle of groundwater withdrawal and recharge, this assumption is most likely to be true since

groundwater levels over the past two decades have been consistently higher than levels in the first half of

the twentieth century (Woodside & Westropp, 2015). Additionally, the common driving mechanism for

inelastic deformation associated with compressible materials typically occurs over a long time span,

although rapid, substantial stress increases can also lead to rapid inelastic deformation. In order to test the

hypothesis of purely elastic deformation, we must compare short-term deformation signals with short-term

variations in hydraulic head and long-term deformation signals with long-term variations in hydraulic head.

As discussed in section 1, the observed deformation on the ground surface is a result of the integrated com-

paction of aquifers and aquitards along the entire depth of the aquifer system. However, we also know that

there is a time delay for a pressure perturbation to diffuse throughout an aquifer system (e.g., Figure 5). This

time delay is dependent on hydraulic diffusivity along the diffusion path and the rate of pressure change at

the pressure source (e.g., groundwater withdrawal rate). Therefore, different regions of the aquifer systems

will be compacting/expanding at different times. By comparing ground deformation time series with

hydraulic head time series at various depths, we can estimate the effective depth at which aquifer deforma-

tion is most correlated with ground deformation. In the ideal case, this depth would approximate the depth

of a hypothetical localized pressure source.

We examine vertical displacement data from the GPS station SACY, which is part of SCIGN and is located

about 1.4 km away from the SAR-9 well (Figure 3). We decompose both the GPS vertical displacements and

the SAR-9 head time series into long-term and short-term (seasonal) signals using a modified form of equa-

tion (13). Instead of using sinusoids to model the seasonal signals, we use a linear combination of third-

Water Resources Research

10.1029/2017WR021978

RIEL ET AL.

3565

order B-splines (different than the integrated B-splines used for transient signals) to allow for seasonal sig-

nals with time-varying amplitudes and phase delays. We assign the temporal support of the B-splines such

that the seasonal signal each year is described by a linear combination of B-splines spaced 0.2 years apart.

With this approach, we could reconstruct seasonal signals with wide variations from year to year. In this

study, we construct

in equation (14) such that B

-splines are independent while B-splines are correlated

with other B-splines that share the same centroid time within any given year (e.g., B-splines centered in

March are correlated with other B-splines centered in March). We assign the correlation strength for the B-

splines to be exponentially decaying in time with a decay time of two years. The decay time was experimen-

tally chosen in order to maintain the flexibility of the B-splines to model time-varying seasonal signals while

still enforcing a level of coherency from year to year.

After decomposing the SACY vertical displacements and the SAR-9 hydraulic head time series at multiple

port depths into long-term and short-term components, we compute the Pearson correlation coefficients

between each head time series and the GPS time series. We fit the correlation coefficient depth profiles

with third-order polynomials in order to reduce the noise of the coefficient estimation and estimate the

depth of maximum correlation between hydraulic head and ground deformation. We then compute an

optimal scaling factor between the GPS and head time series at the depth of maximum correlation (analo-

gous to

in equation (6)).

For both the long-term and short-term components, the vertical ground displacement data is well matched

by hydraulic head variations at a given depth (Figure 7). The transition from regular groundwater fluctua-

tions (caused by the annual cycle of recharge and withdrawal) to unsteady oscillations between 2008 and

Figure 7.

Cross correlation analysis between SAR-9 head data and SACY GPS data. The upper plots correspond to short-

term, seasonal signals for both data sets while the lower plots correspond to long-term signals. We iterate over the time

series at each port depth, compute the Pearson correlation coefficient between well and GPS data, and estimate the

storage coefficient (scaling parameter) to best match the well and GPS data. We perform 100 bootstrap trials with random

subsets of the data in order to estimate uncertainties in the correlation coefficients. The right plots show the mean

Pearson correlation coefficient (and 95% confidence interval) between the time series at each port depth and the SACY

GPS time series. The dashed black line corresponds to a third-order polynomial fit to the mean correlation coefficients,

and the red star indicates the location of maximum correlation. The depths of best correlation are different for the

short-term and long-term signals. For the left plots, SACY data are shown with blue dots and the scaled well data for the

depth of maximum correlation are shown solid red lines. The optimal scaling factors are indicated by

. The red shaded

area represents the uncertainty in the scaling of the well data.

Water Resources Research

10.1029/2017WR021978

RIEL ET AL.

3566

2012 (caused by a change in annual pumping practices) can be observed in the short-term GPS signal. Addi-

tionally, the long-term decrease in water levels and two periods of lower pumping activity due to greater

water availability are manifested as vertical subsidence and uplift, respectively. We find that the seasonal

ground deformation is best correlated with hydraulic head variations at a depth of about 450 m. Interest-

ingly, the depth of best correlation for the long-term ground deformation is deeper (

580 m). In the next

section, we perform a similar comparison of hydraulic head to ground deformation using an 18 year InSAR

time series. The spatially dense observations provided by InSAR allow us to assess the time-dependent

ground deformation at every well location, which is an advantage over the sparse GPS network over the

coastal basins.

4. Central and Santa Ana Basin InSAR Time Series

We use 165 SAR acquisitions from the European Space Agency ERS (European Remote Sensing) and Envisat

satellites spanning from 1992 to 2011 to form 881 interferograms. The maximum perpendicular baseline

(spatial separation between orbits) is 480 m, which is reasonable for isolating ground deformation over met-

ropolitan areas for C-band SAR instruments. The temporal repeat times range from 35 to 210 days, although

after 1995, the repeat times are generally 35 or 70 days which is sufficient to model most deformation sig-

nals observed in the GPS data. We use a coherence threshold of 0.4 to mask poorly resolved areas (such as

over water) and any areas with unwrapping errors. Interferometric phase contributions due to topography

are removed using a digital elevation model (DEM) produced by the Shuttle Radar Topography Mission

(SRTM) with approximately 30 m spacing (Farr et al., 2007). We estimate and remove phase delays due to

atmospheric effects using global atmospheric reanalysis data from the European Center for Medium-Range

Weather Forecasts (ECMWF; Jolivet et al., 2014). We also remove long-wavelength signals due to orbital

errors by estimating a two-dimensional linear ramp for each interferogram. Finally, we reference the time

series to a 400 m

3

400 m window colocated with the SCIGN GPS station SNHS which is located in an area

of high coherence showing stable ground motion unaffected by groundwater withdrawal (Figure 3).

4.1. Seasonal Amplitude and Phase Maps

We first reconstruct the InSAR time series using the same time parameterization approach as equation (13)

where we assume the ground deformation can be described as a superposition of sinusoidal seasonal and

transient effects and estimate the sinusoidal and

-spline coefficients independently for each pixel (Hetland

et al., 2012). We limit our initial analysis of the InSAR time series to interferograms prior to 2008 to reduce

the effect of changes in annual pumping practices on the time series reconstruction. In supporting informa-

tion section S1 and section 5, we describe our method for performing a fully spatiotemporal time series

analysis for the full time series to account for nonsteady seasonal and transient deformation. At this point,

our primary goal is to examine the characteristics of the steady seasonal deformation prior to 2008, and our

experiments show that a pixel-by-pixel approach is suitable for estimating the coefficients of the sinusoidal

components. We can then generate maps of seasonal amplitude and phase using equation (15) for each

pixel.

Maps of the estimated amplitude and phase for the seasonal signal between 1992 and 2008 show that

most of the seasonal deformation is concentrated within the Los Angeles and Santa Ana Coastal Basins in

the region corresponding to the confined aquifer system (Figure 8). The maximum peak-to-peak amplitude

is 5 cm in the southern end of the Santa Ana Coastal Basin, which agrees with the results obtained by Baw-

den et al. (2001), Watson et al. (2002), and Lanari et al. (2004). We can also observe a smaller pair of high-

amplitude regions closer to Long Beach with amplitudes of 3 cm. Here withdrawal occurs primarily at older

sequences that have been folded up closer to the land surface. Low vertical permeabilities at these older,

deeper sequences result in large water-level variations at these pumping sites in response to pumping

stress. The seasonal amplitude decreases rapidly outside of the confined aquifer systems, particularly in the

western edge of the basin bounded by the Newport-Inglewood Fault (NIF) where the fault is an effective

barrier to across-fault fluid flow. This effect can also be observed in the map of the seasonal phase where

the ground east of the NIF has a peak signal in March whereas the ground west of the NIF has a peak signal

in July. The seasonal amplitude decreases rapidly to the north toward the forebay in Los Angeles where

aquifers are semiconfined/unconfined, reducing the water-level response to groundwater withdrawal.

Water Resources Research

10.1029/2017WR021978

RIEL ET AL.

3567

In general, the seasonal amplitude appears to be inversely correlated with the seasonal phase in the center

of the basin, i.e., higher amplitude areas peak earlier in the year, which suggests that groundwater dynamics

here follow a diffusion process. The central areas of the basin experience the highest amplitudes and earli-

est peak times, and we observe delays as one moves toward the margins of the basin. From section 3.2, we

observed a similar delay in hydraulic head from the well data, suggesting that the main driver of the delay

in ground deformation in the central regions of the basin is the time delay necessary for aquifer pressures

to equilibrate in the horizontal direction in response to pumping. Furthermore, in the northwest area of the

basin where aquifer thicknesses are relatively uniform, we can observe an exponential decline in seasonal

amplitude and linear variation in seasonal phase away from an amplitude peak (A–A

transect in Figure 8),

which agrees with the diffusion solution for periodic head variations in equation (11). As with the well data,

Figure 9.

Map of seasonal phase delay (left) with contour lines and labels corresponding to the depth of the bottom of

the principal aquifer system in meters. The dashed white line corresponds to the approximate boundary between the

forebay and confined areas as defined by OCWD. The solid blue line corresponds to the Santa Ana River. The right plot

shows the transects (C–C

) for the seasonal phase (solid blue) and aquifer depth (dashed black). While there appears to be

a general trend between depth of the principal aquifer system and seasonal phase delay, several other factors control

spatial propagation of hydraulic head (see section 6.3).

Figure 8.

Maps of seasonal peak-to-peak amplitude and phase delay. The dashed black lines show the location of the

transects used for the bottom plots. The dashed white line indicates the approximate boundary between the forebay and

confined areas. The arrow indices the satellite-to-ground line-of-sight (LOS) direction. The red circle in the map of sea-

sonal phase delay shows the location of a discontinuity in amplitude and phase within the basin. For the lower transect

plots, blue lines correspond to the phase delay while the dashed red lines correspond to the peak-to-peak amplitude. In

general, the deforming areas of the basin correspond to the confined aquifer system.

Water Resources Research

10.1029/2017WR021978

RIEL ET AL.

3568

the boundaries of the deforming areas become more difficult to interpret due to the transition from con-

fined aquifer units to unconfined/semiconfined units.

We also observe a sharp discontinuity in both the seasonal amplitude and phase maps on the eastern edge

of the basin where the Santa Ana River enters the forebay region indicating some form of impediment to

groundwater flow. Here the peak amplitude occurs in March on the west side of the discontinuity and in

May on the east side. In section 3.2, we discussed the possible influence of specific yield contrasts and faults

acting as impediments to fluid flow. While no known faults exist in this area, this area also corresponds to

the approximate boundary between the forebay and confined zones. Comparison of the seasonal phase

map to the thickness of the principal aquifer system (using the aquifer model developed by OCWD) shows

that the phase discontinuity is also roughly coincident with large changes in depth of the principal aquifer

system in the OCWD model across a short distance (Figure 9). In section 6.3, we explore dynamic models

investigating the impact of variations in aquifer thickness on propagation of hydraulic head.

4.2. Groundwater Withdrawal and Seasonal Ground Deformation

We expect that groundwater pumping practices would have a strong impact on the amplitude of seasonal

ground deformation. Hydraulic head in a confined aquifer system near a production well will experience a

drawdown during periods of groundwater withdrawal (Middleton & Wilcock, 1994; Todd & Mays, 1980). To

explore the impact of withdrawal on head levels, we use groundwater production time series for 250

OCWD production wells that measure total groundwater withdrawal on a monthly basis from 1996 to 2007.

We create a 50

3

50 uniform grid where the dimension of each grid cell is approximately 0.6

3

0.6 km and

compute the amplitude of total seasonal groundwater withdrawal in each grid cell. We perform this proce-

dure for three different depth ranges (0–200, 200–370, and 370–700 m depth) to explore the depth depen-

dency of withdrawal. We can clearly observe an association between the high ground deformation areas

and areas with high seasonal groundwater production at depths mostly associated with the principal aqui-

fer system (Figure 10). The pair of high-amplitude ground deformation regions near Long Beach directly

correspond to two regions of concentrated groundwater withdrawal between 200 and 370 m depth. How-

ever, not all areas of high groundwater withdrawal are colocated with high ground deformation, particularly

Figure 10.

Maps of (a) seasonal peak-to-peak LOS amplitude and seasonal groundwater withdrawal from OCWD produc-

tion wells located at depths of (b) 0–200 m, (c) 200–370 m, and (d) 370–700 m. See supporting information Figure S2.4 for

the corresponding seasonal phase delays. The locations with the highest seasonal withdrawal at 200–370 m depth corre-

spond to the locations of the highest seasonal ground deformation within the pressure area. The red outline encom-

passes the seasonal ground deformation as observed in the InSAR map of seasonal phase delay. The dashed white line

indicates the approximate boundary between the forebay and confined areas. In general, the areas of highest seasonal

ground deformation correspond to locations where confined aquifers are most heavily pumped.

Water Resources Research

10.1029/2017WR021978

RIEL ET AL.

3569