Aquifer deformation and active faulting in Salt Lake Valley, Utah, USA

Abstract Aquifers and fault zones may interact through groundwater flow and stress redistribution, yet their spatiotemporal relationship remains enigmatic. Here we quantify changes in water storage and associated stress along the Wasatch Fault Zone in Salt Lake Valley, recently shaken by a M5.7 earthquake on March 18th, 2020. Ground deformation mapped by Sentinel-1 SAR imagery (2014-2019) reveals an elongated area with ∼50-mm seasonal uplift corresponding to 0.03-0.06-km3 water storage cycles. Phase shifts across active faults in both water level and deformation suggest control by the low-permeability structures. The seasonal stress changes on the adjoining faults from poroelastic volume strain are two orders of magnitude larger than those from hydrological surface loading, but both are small compared to the annual increase of tectonic loading at seismogenic depths. Historic seismic events, limited in number, do not exhibit statistically significant annual periodicity and hydrological modulation of microseismicity or triggering of the recent M5.7 event is not evident.


Introduction
Natural water discharge (e.g., evaporation and drainage) and recharge (e.g., rainfall and snowmelt infiltration) maintain a sustainable hydrosphere and ecosystem. In particular, aquifers help regulate the water balance by storing and releasing the groundwater as needed. Such natural subsurface reservoirs are invaluable in arid regions where freshwater resources are limited. Human extraction of groundwater is sustainable, if net extraction is balanced by recharge and water levels can be maintained at stable levels.
Land subsidence is often observed over sedimentary basins due to water level decline and gradual consolidation of the confining units and fine-grained silts and clays that constitute the interbeds. Subsidence may be large (up to several meters), permanent and unrecoverable (inelastic), if the water head drops below previously achieved lowest levels and the stress exceeds preconsolidation conditions (Galloway and Burbey, 2011;Ojha et al., 2018Ojha et al., , 2019. Cyclic seasonal subsidence and uplift by millimeters to centimeters are typically associated with water discharge and recharge producing poroelastic deformation (e.g., Amelung et al., 1999;Lu and Danskin, 2001;Chaussard et al., 2014;Hu et al., 2018;Carlson et al., 2020). Horizontal movements also exist and generally occur in the vicinity of operating wells, near fault zones traversing aquifers, and along the margins of aquifer basins (e.g., Chaussard et al., * Corresponding author. E-mail address: xiehu@berkeley.edu (X. Hu). 2014; Helm, 1994;Fu et al., 2013). Consideration of the horizontal movements can improve our ability to quantify the properties and geometry of subsurface aquifer systems (Burbey, 2008).
Hydrological loading and unloading may regulate seismicity through elastic stresses in the seismogenic zone (e.g., González et al., 2012;Johnson et al., 2017;Craig et al., 2017). In addition, poroelastic stresses due to subsurface pore-fluid pressure diffusion driven by precipitation and/or groundwater variations may also contribute to modulating seismicity, at least at shallow depths and in especially permeable rocks (e.g., Hainzl et al., 2006;Montgomery-Brown et al., 2019;Wetzler et al., 2019). Anthropogenic oil and gas production and fluid injection may also trigger earthquakes through pore pressure redistribution (e.g., Ellsworth, 2013;Shirzaei et al., 2016;Goebel and Brodsky, 2018). How natural groundwater processes in smaller sedimentary basins can affect seismic hazards remains an open question.
Salt Lake Valley, Utah is a sedimentary basin that hosts the commercial, industrial and financial state capital Salt Lake City. Three-fourths of the state's population (∼3 million) is concentrated within a 160-km radius of the city. The valley is bounded by the generally NS-trending Oquirrh Mountains to the west, the Wasatch Range to the east, and the EW-trending Traverse Mountains to the south. The 70-km-long Jordan River traverses the central axis of the valley, connecting two remnants of prehistoric Lake Bonneville (30,000-14,000 yr BP) -Great Salt Lake and Utah Lake. The basins are composed of three distinct hydrological units ( Fig. 1): the water discharge area with an upward hydraulic gradient in the lower-elevation northern part of the confined basin and a narrow  (Herring et al., 2016). The error ellipses represent 95% confidence intervals. Black lines are the Quaternary faults. (b) The horizontal strain-rate field determined from the GPS velocities. Arrows represent the direction of the principal strains. Dilatational strain (blue) governs most parts of the eastern Basin-Range. Our study area, Salt Lake Valley (SLV), is highlighted by a dashed box in the center of panels a and b. (c) A close-up view of SLV. Dashed black lines delineate the boundary of basin-fill deposits. Blue line shows the Jordan River. Solid black lines show the faults. Faults in the SLV include the West Valley fault (WVF), the East Great Salt Lake fault (EGSLF), and the Wasatch fault zone (WFZ), including three major Salt Lake City segments -the Warm Springs fault (WSF), East Bench fault (EBF), and Cottonwood fault (CF). The epicenter of the M5.7 Magna earthquake west of the WVFZ is shown by its normal-faulting focal mechanism. (For interpretation of the colors in the figure(s), the reader is referred to the web version of this article.) unconfined zone bounding the Jordan River; the primary recharge area at the foot of the mountains where the hydraulic head gradient is downward; and the secondary recharge area in between where the confined and unconfined layers are not clearly distinguished (Thiros et al., 2010).
The alluvial basins also host the parallel and sub-parallel N20 • W trending Wasatch fault zone (WFZ) along the front of the Wasatch Range and the inner-valley West Valley fault zone (WVFZ), which make Salt Lake County one of the most seismically hazardous metropolitan areas in the interior of the western U.S. (Wong et al., 2002;Valentini et al., 2020). The 390-km-long WFZ extends from Malad City, Idaho, to Fayette, Utah along the western flank of the Wasatch Range, and separates the stable Rocky Mountains and Colorado Plateau to the east and the extending crust of the Basin and Range Province to the west (Fig. 1). The regression of Lake Bonneville and the deglaciation of mountain ranges around the WFZ during the Late Pleistocene to Early Holocene epochs caused lithospheric rebound and accelerated the slip rates to ∼1 mm/yr, about twice as high as the average geologic slip rate on a 10 5 years time scale (Friedrich et al., 2003;Hetzel and Hampel, 2005;Hampel et al., 2010).
Three en-echelon fault segments of the WFZ surrounding Salt Lake City include the Warm Springs fault (WSF), the East Bench fault (EBF), and the Cottonwood fault (CF) (Moschetti et al., 2017). The Salt Lake City segment of the WFZ is believed to produce large earthquakes (M 7.0+) every 1,300 to 1,500 years, and the last one occurred about 1,300 ± 200 years ago (DuRoss and Hylland, 2015).
The Utah Geological Survey and U.S. Geological Survey (2016) forecast a 93% likelihood of one or more moderate earthquakes of magnitude 5 or greater striking the Salt Lake Valley (SLV) in the next 50 years. Thus, the recent M5.7 Magna, Utah earthquake on March 18th, 2020 ( Fig. 1c) was not a complete surprise. Earthquake hazard stems not only from the shaking, but also the potential liquefaction in lowland areas, and tsunami and seiches in Great Salt Lake if extensive ground subsidence were to occur due to rupture along the East Great Salt Lake fault (EGSLF) (Earthquake Engineering Research Institute, 2015). The 1997The continuous and 1992The -2003 campaign GPS observations of horizontal motions in a stable North America reference frame indicate ≤∼1.6 mm/yr of extension across the WFZ (Chang et al., 2006). Other GPS-based studies found somewhat different rates, depending on the dataset and approach used to determine deformation and slip rates along the WFZ (e.g., Niemi et al., 2004;Puskas and Smith, 2009). The dilatational strain rates calculated from Plate Boundary Observatory network under the North America reference frame (Herring et al., 2016) show an accumulation of extension at 0.1 μstrain/yr (Fig. 1). However, limited by the sparse distribution and inconsistent surveying time among the stations, GPS measurements alone are insufficient for the basin-wide characterization of deformation. Interferometric synthetic aperture radar (InSAR) provides complementary geodetic observations to monitor the spatially continuous crustal deformation with weekly to monthly updates, though the measurements are limited to one-dimensional line-of-sight (LOS). Here, we compile ascending (AT122) and descending (DT100) Sentinel-1 imagery (2014-2019) (Fig. S1), and continuous GPS observations (Blewitt et al., 2018;Fig. S2) to decipher the multi-annual and seasonal vertical and horizontal motions in the SLV.
Wells provide a direct window into the subsurface hydrology. We use 1931-2019 well data from the U.S. Geological Survey (https://waterdata .usgs .gov /usa /nwis/) (Table S1). Earthquake catalogues help us assess the potential effects from spatiotemporally variable stressing patterns. To assess spatio-temporal variations in seismicity, we draw on the decadal earthquake catalog from 1981 to 2018 provided by the University of Utah (see supplement for details). We use a joint analysis of geodetic displacement measurements, water levels ( Fig. S3; Table S1) and earthquake information ( Fig. S4; Table S2) to quantify the seasonal variation in water storage, estimate the commensurate stress changes on nearby faults, and explore the potential coupling between the hydrological and tectonic processes in the SLV.

Separation of vertical and horizontal motions from temporal behaviors
We use ascending (AT122) and descending (DT100) Sentinel-1 tracks to resolve the multiannual and seasonal temporal behaviors of ground deformation during 2014-2019. The time-series analysis is performed on persistent scatterers in small baseline subset interferograms (see supplement for details). Different timings of seasonal motion make it challenging to separate the horizontal and vertical components from the characteristic seasonal amplitudes, yet this also provides an opportunity. We ignore NS motions to which spaceborne SAR systems with near-polar orbits are insensitive as the heading angle is ∼10 • from due north/south (see supplement for details). The right-looking ascending and descending observations have similar sensitivity to vertical motions and contrasting representations of EW motions from opposing flight directions, i.e., ascending orbits look down to the east and descending orbits look down to the west. If the seasonal motion is dominated by the EW component, the peak motions in ascending and descending time series will be out of phase. Therefore, we can identify the targets that mainly move seasonally in the vertical component, as their peak motions in ascending and descending time series occur at the same time of the year (here within 30 days).

Analytic modeling of the seasonally deforming aquifer system
We use an analytical solution of finite strain sources in a halfspace (Barbot et al., 2017) to fit the observed InSAR displacement fields during seasonal groundwater-level changes. Strong vertical seasonal movements are mainly contained in the asymmetrical elongated area of interest (referred to as AOI hereafter), shown by white dashed outlines in Figs. 2c, d. We mesh the AOI using an arrangement of 481 grid cells with individual dimensions of 500-by-500 m. The NW trending cuboid boundaries strike N20 • W, sub-parallel to the basin and surrounding fault strands. Here we consider isotropic volume-strain sources reaching up to the surface, assuming that any shallow confining layer is thin. As a firstorder approximation of the isopach map of unconsolidated and semi-consolidated sediments over this AOI (Mattick, 1970), we assume a bulk aquifer thickness of up to 600 m and apply a Poisson's ratio of 0.25.
There is a strong trade-off between the volume strain and thickness of the model cuboids. We thus consider two end-member scenarios of constant volume strain and variable thickness, and variable strain and constant thickness of the cuboid elements to generate best-fit LOS displacement fields, which capture both the vertical and horizontal motions. We focus on the skeleton expansion during the wintertime phase of peak uplift. In the first model, assuming that the vertical motion linearly correlates with the water level and thus the aquifer bulk thickness, we use the InSARresolved vertical seasonal amplitudes to obtain aquifer thicknesses ranging from 0 to 600 m. In the second end-member model, we consider a constant aquifer thickness of 500 m for all the cuboids and compute the spatial distribution of volume-strain to fit the displacement fields. We compare results from these two end-member parameterizations and their consequent bulk volume changes are consistent (see section 4.3 for details).

Elastic deformation due to seasonal hydrological loading and unloading
The storage coefficient describes the amount of water drained from the aquifer per unit decline in water level, and it can be re-solved by a linear correlation between the vertical displacements and water level changes (Chaussard et al., 2014, and references cited there). Based on this concept, Hu et al. (2018) obtained a storage coefficient map of SLV through interpolation from the estimates at four daily sampled water wells using the earlier SAR data sets (2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011). This allows us to estimate the seasonal groundwater level changes using the Sentinel-1 derived vertical displacements. Confined aquifers commonly have a porosity of 0.2-0.4 (e.g., Chaussard et al., 2014). Therefore, the annual changes in water storage amount to ∼0.03-0.06 km 3 for areas with equivalentwater-thickness change larger than 0.5 m, where the estimated seasonal variation in the groundwater level is up to 3 m. We simplify the subsurface scenario to surface loading and unloading in a Cartesian elastic half-space (Becker and Bevis, 2004). We use a Poisson's ratio of 0.25 and a shear modulus of 3 GPa for the weakly lithified sedimentary rocks.
We apply a simple 2D solution to estimate the stress change due to an equivalent distributed line load distributed across the ∼8-km-wide aquifer. Then we infer the equivalent distributed line load for the seasonal groundwater mass change (see section 4.4). The stress components at a given depth in an elastic half-space can be represented by the angular distance from the loading source (Jaeger et al., 2007;Amos et al., 2014;Kundu et al., 2015). We consider a characteristic fault dip of 55 • (Chang et al., 2006) and frictional coefficient of 0.4 to estimate the Coulomb stress change along the EW profile traversing the loading source. Note that the stress values associated with volume strain and the elastic loading depend directly on the applied elastic moduli.

Regional seasonal and multi-annual deformation in space and time from InSAR
We extract targets whose seasonal movements are predominantly vertical by correlating the timing of the LOS motions from ascending and descending data. Annual uplift peaks in the timeseries displacements are generally around March to April due to abundant recharge from mountain snowmelt. Our measurements can be validated by GPS station SUR1, the only site located inside the AOI in its southern part (Fig. 2). Although SUR1 has less than 2 years of data available (1997)(1998), the seasonal peak-to-peak motions are well resolved with amplitudes in the EW, NS and UD components of about 5.2, 2.9, and 27.4 mm, respectively (Fig. S2). For comparison, the values for stations ZLC1 and SLCU at the margins of the AOI are about 11.4, 3.5, and 12.8 mm, and 14.1, 7.7, 12.8 mm, respectively (Fig. 4), and the EW and UD components have similar amplitudes that are much larger than in the NS direction. We project the 3D GPS displacements of ZLC1 into the Sentinel-1 LOS ascending and descending directions for comparison of time-dependent motions during 2014 to 2019 (Figs. 2e, f). We apply a constant offset adjustment between them for a best fit to account for their different starting times and reference systems. The GPS and InSAR time series match well with residuals of 2.58 and 1.00 mm for the AT122 and DT100 tracks, respectively. The time of year of peak uplift (∼4 mm) at GPS sites on the surrounding ranges (i.e., COON and RBUT) occurs in summer/fall, while that for GPS sites within the basin-fill deposits occurs in winter/spring with larger amplitude (Fig. S2).
We retrieve the 2D displacement maps (EW and vertical; Fig. 3) for seasonal amplitudes and multi-annual velocities from ascending and descending Sentinel-1 InSAR, assuming that the NS displacements are negligible (see supplement). The AOI presents pronounced seasonal motions in both horizontal and vertical components with sharp margins (Figs. 3c, d). The AOI uplifts by ∼50 mm from fall to spring, accompanied by EW extension with a for targets whose seasonal amplitude is larger than 1 mm. Colored circles (with station labels in c and d) represent the amplitude and phase information obtained from the time series of the vertical GPS component; unfilled circles are stations whose seasonal uplift was not resolved due to short time spans. White dotted lines highlight the area with seasonal motions dominated by the vertical component. (e) and (f) compare Sentinel-1 LOS and GPS time series during 2014-2019. GPS station ZLC1 is the only station within the basin that overlaps in time with the Sentinel-1 observing period (Fig. S2). The 3D GPS displacement time series at ZLC1 in gray symbols have been projected into radar LOS directions for comparison with ascending track AT122 (red circles) and descending DT100 (blue circles) LOS time series at the same location. We apply a constant offset adjustment between InSAR and GPS results due to their different start times and reference systems. The yellow star denotes the reference area.
net horizontal motion of ∼30 mm across the uplift zone. The displacements reverse for the other half of the year from spring to fall, with subsidence and EW shortening of the same magnitude. This N20 • W oriented zone of hydrological deformation has a larger (∼600 m) sediment thickness than the surrounding areas (Mattick, 1970). The seasonal deformation zone is bounded by the WVFZ and EBF in the north, while the southern end without such bounding structures appears more diffuse in its deformation pattern. Hydrogeologically, the AOI is part of the water discharge unit. The Jordan River cuts longitudinally through the central AOI and divides the horizontal displacement field into several smaller, isolated patches.
The long-term displacement map for 2014-2019 reveals that the eastern half of the valley is subsiding at ∼1-2 mm/yr relative to the western SLV. The spatial distribution of longer-term ground subsidence coincides with the areas of largest water level decline and S5), consistent with the rates determined from Envisat ASAR spanning 2004-2010 (Hu et al., 2018). The seasonal displacement field highlights a local area experiencing highly variable aquifer storage, whereas the multi-annual displacement field presents a regional long-wavelength signal correlated with prolonged water drawdown.

Temporal variations of water levels and 3D GPS observations over the basin
While the temporal sampling of water-level measurements is sparse, we are able to determine the phase and amplitude of average annual variations for some of the wells in the SLV region. The timing of the seasonal water level fluctuations varies among wells at different locations with phase shifts of several months (Figs. 4b and S6). Wells located on either side of the EBF represent remarkably contrasting patterns in time. Artesian wells 23301 and 30901 in the water discharge area to the west of the EBF have the lowest water level from June to August, likely due to summer pumping. In contrast, this time period features the highest water levels at wells 94001 and 03901 on the east side of the fault and in the water recharge area at the foot of the ranges (Fig. 4b). Other wells distributed across the basin have varying temporal patterns that depend on their location with respect to the principal recharge and discharge zones and faults (Figs. S3 and S6).
To further investigate the controls of the orientation and timing of seasonal displacements, we focus on GPS time series from three stations close to the target aquifer that overlap in time (Fig. 4). Stations ZLC1 and SLCU are ∼1.5 km apart and located east of the WVFZ in the northeast portion of the AOI and within the water discharge area (confined aquifer), while UTCR is located on the southwestern edge of the AOI in the secondary recharge area just west of the WVFZ (undistinguished confined-unconfined aquifer). Seasonal uplift of UTCR is accompanied by southwesterly motion, whereas the uplift of ZLC1 and SLCU is accompanied by northeasterly motion, as expected for the expansion of a finite elastic porous medium (e.g., Burbey, 2008). Interestingly, the seasonal displacements observed in those two groups are shifted by ∼4 months: ZLC1 and SLCU have the largest subsidence in fall, in contrast to UTCR with peak subsidence in spring-summer. As for the 3D displacements of UTCR, the smallest horizontal motion (most southwesterly position) occurs up to 4 months earlier than that of the vertical component, while no evident difference in phase between the vertical and horizontal motions exists for the other two sites. Overall, the time-series GPS observations illuminate phase differences in 3D seasonal motions depending on the location in the groundwater basin, but the small number of stations limits our ability to make out systematic patterns in this behavior.

Relationship between seasonal water levels and GPS-/InSAR-derived displacements
We attribute seasonal deformation patterns captured by the GPS and InSAR time series to annual variations in water storage in the SLV groundwater system, which is also reflected in the changing well water levels. Well 75901, southwest of the AOI and within the secondary recharge area (Fig. 2), is the only one that has daily sampled water levels during our observation period. The seasonal LOS displacements are modest at this site (Fig. 5). The peak displacements measured by both tracks are a few weeks prior to that of the water level. This may be because this well taps water at a depth of 242 m, above which there may be additional deforming layers whose water levels change earlier than the deeper aquifers.
The colocated GPS-derived ground motions and well water levels are correlated. For example, UTCR and its closest well #75901 reach their minima around May to June (Fig. 4), and the UTCR phase for 2011-2014 is consistent with that for Sentinel-1 in 2014-2019 (Fig. 5). The regional storage coefficient at SLV is between 0.002 and 0.07, and ∼0.024 near downtown Salt Lake City (Hu et al., 2018). Referring to the seasonal vertical displacement of the AOI, we estimate that the principal aquifer experiences up to ∼3 m of seasonal water-level variations, corresponding to up to ∼1 m of equivalent-water-thickness and seasonal change in water storage of ∼0.03-0.06 km 3 , considering a porosity of 0.2-0.4. Such hydrological loading can produce up to 6 mm elastic subsidence in the spring (reversed for the unloading scenario; Fig. S7) (Becker and Bevis, 2004), which is negligible compared to the 50-mm vertical motion due to the expansion and contraction of the poroelastic aquifer skeleton (Fig. 3c). The load model predicts up to 2 mm of convergence across the aquifer during peak spring loading (Fig.   S7), compared to ∼30 mm of extension produced by the aquifer strain (Fig. 3d). Note that the direction of the vertical motions from these two physical processes associated with seasonal water storage change (i.e., poroelastic volume strain and elastic loading), are opposite of one another.

Deformation from elastic loading vs. poroelastic aquifer strain
The timing difference in cyclic ground motions between the mountain ranges and the adjacent unconsolidated alluvial basins is well understood as a consequence of their distinct controlling mechanisms (e.g., Amos et al., 2014). Elastic loading and unloading by snow and water result in instantaneous ground subsidence and uplift as illustrated by the GPS stations located on the ranges (RBUT and COON in Figs. 2 and S2). On the other hand, groundwater inflow and outflow in the basin environment cause poroelastic uplift and subsidence, respectively, generally with delays due to diffusion of water into and out of the aquifer and/or inelastic compaction processes. When hydraulic head declines, groundwa-ter outflows from pore spaces in the fine-grained interbeds and confining units, and thus the compressible materials elastically compact and the land subsides. The opposite phenomenon occurs when hydraulic head increases, raising pore fluid pressure and decreasing the effective elastic stress on the granular skeleton supporting the vertical load (e.g., Rice and Cleary, 1976;Chaussard et al., 2014). Therefore, land surface elevations above the aquifer reach maxima during snowmelt runoff from the mountains and reach minima when groundwater levels are depleted by surface and subsurface flow, pumping, and evaporation.

Role of fault-aquifer interaction
Multiple lines of evidence suggest that faults may act as physical boundaries, defining and perhaps controlling groundwater re- distribution. In the spatial domain, the margins of the AOI agree with the extent of the confined water discharge area and nearby active fault traces. The Jordan River cutting through the AOI longitudinally also affects the groundwater system and complicates the displacement field in the center of the AOI. In the temporal domain, different sides of the fault splays have distinct phase patterns in their seasonal motions and also in water level (Figs. 4 and S6). The faults and fractures at depth may act as low-permeability barriers to horizontal flow, so the groundwater flow is regulated but not completely obstructed. This may be the reason for the observed phase shift by several months of the water levels on either side of the EBF (Fig. S6), and phase differences in ground motions between the two sides of the WVFZ (Fig. 2); similar phenomena have been observed in the Los Angeles basin (Bawden et al., 2001) and the Santa Clara Valley (Chaussard et al., 2014) in California.

Estimating water storage and volume strain changes
To estimate the water storage changes and quantify the stress contribution from the seasonal deformation of the aquifer system, we rely on an analytical solution of finite strain volumes in a halfspace for cuboid sources (Barbot et al., 2017). We investigate two end member scenarios in which we fix either the thickness or the volume strain rate of each strain source. In the first model with variable thickness ranging from 0 to 600 m, a homogeneous isotropic strain of 9.1 × 10 −5 yields 3D displacements that best fit the ascending and descending InSAR observations. In the second model of constant thickness of 500 m, we invert for the distribution of volume strain that produces a displacement field that best fits the InSAR results, and the resulting strains range from ∼2-12 × 10 −5 .
Both models can fit the observations well as shown in crosssectional profiles (Fig. 6), but the variable-strain model captures the details better with an overall residual of 3.49 mm, compared to 3.74 mm for the variable-thickness model. We thus prefer the model with the variable strain and focus on that to compute the Coulomb stress change assuming pure normal faulting on the principal fault planes of the WFZ near Salt Lake City, i.e., the Warm Springs fault, East Bench fault, and Cottonwood fault (Moschetti et al., 2017). Ascending-orbit results can be better recovered in both models, and the displacements in the western half of the AOI in the descending results are underestimated. This may imply more complexity and spatial heterogeneity in this subsection in both the source strain and thickness, as well as anisotropic volume strain, which seem likely in the context of multiple parallel faults and their southern ends over this area.
The distribution of the uplift and thus the bulk thickness in the first model and that of the strain in the second model are very similar, suggesting consistent vertical integration of the strain sources. The consequent seasonal bulk volume changes for these two models are estimated to be 3.3 × 10 −3 and 3.5 × 10 −3 km 3 , respectively, similar to the product of the previously estimated representative storage coefficient (∼0.024) and the volumetric variation of the water-bearing unit (0.15 km 3 ).

Stress changes from volume strain and elastic surface loading
Using the volume-strain sources from the variable-strain model inverted from the seasonal deformation data, we can forward model the seasonal changes in stress on nearby faults, assuming a shear modulus of 3 GPa for the young basement. Accompanying annual surface uplift of ∼50 mm and water storage increase of 0.03-0.06 km 3 , the estimated Coulomb stresses on the dipping fault planes (WSF, EBF and CF) change by about −450 to 50 kPa at shallow depth (<∼600 m) during the wintertime (peaking in March); the stress changes reverse for the summertime (Figs. 7c and S9). In the normal-faulting regime, larger earthquakes tend to nucleate near the brittle-ductile transition zone (>10 km) and propagate upwards (Roten et al., 2011). At these depths, the stress perturbations from the nontectonic aquifer strain are about −10 to 2 kPa during the wintertime. Overall, the seasonal stress changes in the spring are dominated by negative normal stress changes (clamping) underlying the aquifer and larger positive normal stress changes (unclamping) at the sides on dipping faults (Fig.  S9). The seasonal stress changes at seismogenic depths due to shallow aquifer processes generally lie below estimates of the annual background loading rate on the WFZ (∼15 kPa/yr, Bagge et al., 2019; ∼3.6 kPa/yr, Verdecchia et al., 2019). Note that a wide range of elastic moduli in the natural Basin and Range setting brings uncertainty to the absolute values of our stress-change estimates.
In addition to aquifer deformation, seasonal stress variations also result from other hydroclimatic periodic sources, including elastic water loads, atmospheric pressure, temperature, and Earth pole tides (e.g., Johnson et al., 2017). In California, the largest regional source of seasonal stressing comes from elastic water loads in the form of snow, lakes and groundwater and may periodically increase seismicity rates by nearly 10% (Johnson et al., 2017).
For a first-order estimate of stress changes at depth due to elastic loading, we model deformation and stress from the Salt Lake Valley aquifer storage changes by applying an equivalent distributed line load rate distributed across the width of the deforming aquifer (Jaeger et al., 2007;Amos et al., 2014;Kundu et al., 2015). The volume mass across the distributed line load is given by N 0 = ρ g w n hl c w c , where the water density ρ = 997 kg m −3 ; the gravitational acceleration on the Earth g = 9.8 m s −2 ; the endto-end EW width of the aquifer w is 12.5 km; notably, the EW width of the aquifer model load is 8 km (Fig. 8); each cuboid has dimension l c by w c of 500 by 500 m in NS and EW directions; the porosity n ranges between 0.2 and 0.4; the spatially variable groundwater level changes are labeled as h, and Fig. 8a shows the equivalent water height change nh. Therefore, the seasonal groundwater mass change amounts to (3.04∼6.09) × 10 10 kg and the corresponding distributed line load is (2.38∼4.75) × 10 7 N m −1 .
We find that the Coulomb stress changes during peak spring loading on a fault plane dipping 55 • (Chang et al., 2006) are up to only 0.7 kPa concentrated at depths of 0.4-1.2 km underneath the aquifer (Fig. 8). It is worth mentioning that a west dip of 55 • does not represent the geometry of all seismogenic structures in the SLV. The dip and strike both control the distribution and the sign of Coulomb stress changes at a given location, but their magnitude is always very low compared to background stress and stressing-rate levels (Figs. S11, S12). Overall, the stress change from the elastic volume strain source at shallow depth is more than two orders of magnitude larger than that from the surface loading.

Seismicity analysis to assess role of annual and multi-year stress perturbations
The 1981-2018 earthquake catalog for the SLV contains a total of 635 seismic events, up to M4.16 (Fig. 7a). After declustering the catalogue, we are left with 512 events (see supple- Month-of-year histograms of the full and declustered earthquake catalogue in boxes a and b and their combined areas. The bottom plot shows the seasonal LOS displacements at the center of box a (red and blue shades represent AT122 and DT100 results, respectively). In contrast, the displacements in box b exhibit no seasonality (Fig. S10). The March 18th, 2020 M5.7 earthquake (focal mechanism in A and hypocenter shown as a white star in C). (c) The Coulomb stress change on the WFZ fault planes adopted from Moschetti et al. (2017), and the second invariant of stress at each hypocenter (1981-2018) due to volume strain during peak water levels in the spring. ment; Wiemer, 2001). The major faults of the WFZ do not host a significant number of events. Instead, the northwest SLV contains two major clusters (Fig. 7a). Cluster a is bounded by splays of the WVFZ. Cluster b is separated by the WVFZ and lies ∼7 km west of a and at a greater depth (∼8 versus ∼5 km), in the hanging wall of the deep extension of the WSF and EBF. The time series displacements over cluster a indicate regular seasonal variations with a peak around May (Fig. 7b), whereas motions above cluster b, near the recent M5.7 earthquake and next to a large compacting tailings impoundment (Hu et al., 2017), are fairly stochastic (Fig. S10). The second invariant of the seasonal stress changes ( √ | I 2 |) from the aquifer strain √ | I 2 | at the hypocenters reaches up to ∼20 kPa in cluster a while stress changes are low (<3 kPa) in the more distant cluster b. The March 18th, 2020 M5.7 earthquake (https://earthquake .usgs .gov / earthquakes /eventpage /uu60363602 /origin /detail) is located within cluster b and the springtime Coulomb stress changes on 30 • westdipping or 70 • east-dipping normal faults near the hypocenter are 0.1 kPa and −0.03 kPa due to the surface loading, respectively; while that resulting from the volume strain is 0.36 kPa. Unlike the apparent seasonal variation in seismicity rates due to regional hydrological load cycles in the Nepal Himalayas (Bettinelli et al., 2008), California (Amos et al., 2014;Johnson et al., 2017), and the New Madrid Seismic Zone (Craig et al., 2017), there is no clear indication of annually cyclic seismicity in the SLV (Fig. 7b). While the volume expansion and loading of the principal aquifer peak in spring (Figs. 2b,d), with decreased Coulomb stress concentrated at ∼0-1 km to discourage failure and with increased Coulomb stress at depth of ∼1-4 km to promote failure on the WSF and EBF (Fig. 7c), we are not able to resolve corresponding seasonal changes in seismicity rates in clusters a and b that would support a direct triggering relationship.
On a multi-decadal timescale, while there are temporal variations in both precipitation (proxy for groundwater level) and the number of earthquakes, there does not appear to be a signifi- It is worth mentioning that a west dip of 55 • does not represent the geometry of all seismogenic structures in the SLV. The results for different west and east dips are shown in Figs. S11 and S12. cant correlation (Fig. S4). The limited number of events during four decades over the ∼700-km 2 SLV basin may simply be insufficient to decipher the code of nature with confidence, compared to the significant seasonality seen in orders-of-magnitude larger seismicity catalogs in Nepal, California and New Madrid (Bettinelli et al., 2008;Johnson et al., 2017;Craig et al., 2017). In future work, we hope to explore the role of regional hydrological loading and unloading across the larger Wasatch Range front area, including contributions of regional seasonal snow loads and highly variable levels of the Great Salt Lake.

Conclusions
To sum up, we map out a multi-annual subsidence coinciding with prolonged water level decline in the eastern SLV along the front of the Wasatch Range. We also identify an elongated aquifer following a regular peak-to-peak seasonal uplift (50 mm) and extension (30 mm) during wintertime (reversed for summertime), revealing a seasonal variation in water storage by ∼0.03-0.06 km 3 .
The spatial association of the seasonally deforming area, hydrological discharge units and fault splays, as well as phase shifts in the displacement time series and water levels in areas separated by active faults, indicate that the faults modulate the groundwater flow and poroelastic strain field. The seasonal groundwater breathing of the aquifer exerts up to a few kPa Coulomb stress from the poroelastic volume strain and elastic loading at seismogenic depth of nearby fault zones, generally below the annual increase of tectonic stress. There is currently no evidence to suggest that earthquakes in the SLV, including the March 18th, 2020, M5.7 Magna earthquake, are directly related to the seasonal or multi-year aquifer deformation processes.

CRediT authorship contribution statement
X.H. and R.B. designed the study. X.H. analyzed SAR images and performed the modeling. X.H. drafted the paper with important inputs from R.B.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.