Hydrologically Induced Deformation in Long Valley Caldera and Adjacent Sierra Nevada

Vertical and horizontal components of GNSS displacements in the Long Valley Caldera and adjacent Sierra Nevada range show a clear correlation with hydrological trends at both multiyear and seasonal time scales. We observe a clear vertical and horizontal seasonal deformation pattern primarily attributable to the solid earth response to hydrological surface loading at large‐to‐regional (Sierra Nevada range) scales. Several GNSS sites, located at the eastern edge of the Sierra Nevada along the southwestern rim of Long Valley Caldera, also show significant horizontal deformation that cannot be explained by elastic deformation from surface loading. Due to the location of these sites and the strong correlation between their horizontal displacements and spring discharge, we hypothesize that this deformation reflects poroelastic processes related to snowmelt runoff water infiltrating into the Sierra Nevada slopes around Long Valley Caldera. Interestingly, this is also an area where water infiltrates to feed the local hydrothermal system, and where snowmelt‐induced earthquake swarms have been recently detected.


Introduction
Long Valley Caldera (LVC) is located at the eastern edge of the Sierra Nevada range (SNR) within the transition between the Eastern California Shear Zone (ECSZ) and the Walker Lane Belt (WLB) (Figure 1). This area is characterized by the coexistence of steady tectonic deformation with transient deformation due to hydrology, seismicity, and caldera inflation (Langbein, 2003;Savage & Clark, 1982). LVC is continuously monitored by different instruments, including a dense Global Navigation Satellite System (GNSS) stations network. This monitoring detected uplift episodes in 1997-1998, 2002-2003, and 2011-present, which are generally thought to be related to a volumetric source located beneath the LVC resurgent dome at ∼6-10 km depth (Feng & Newman, 2009;Hill, 2006;Langbein, 2003;Montgomery-Brown et al., 2015).
The SNR receives large amounts of precipitation concentrated in winter/spring, but with noticeable multiyear variability, alternating between periods of surplus precipitation and drought (Margulis et al., 2016). Several studies have highlighted the signature of these precipitation trends in deformation data, showing that observed subsidence/uplift trends at both seasonal and multiyear scales are driven mainly by surface loading/unloading from increasing/decreasing surface water and groundwater (Amos et al., 2014;Argus et al., 2017;Borsa et al., 2014;Kreemer & Zaliapin, 2018). These studies mostly focused on the vertical component of the deformation, since horizontal displacements induced by surface loading are typically much smaller (Wahr et al., 2013).
LVC is impacted by precipitation falling on the adjacent SNR. Its diffuse hydrothermal system is recharged mainly from snowmelt along the western/southwestern caldera rim (Peacock et al., 2016;Sorey et al., 1991), 10.1029 and SNR runoff has been linked to the occurrence of shallow earthquake swarms on the caldera margin (Montgomery-Brown et al., 2019). Furthermore, a recent study suggested a correlation between SNR unloading due to the recent 2012-2016 drought and the most recent inflation episode at LVC (Hammond et al., 2019).
Deformation records, such as those provided by GNSS data, are extensively used to study the different phases of LVC magmatic and seismic activity (Ji et al., 2013;Liu et al., 2011;Montgomery-Brown et al., 2015) and to analyze the effects of hydrological loading on the SNR (Argus et al., 2017;Amos et al., 2014;Borsa et al., 2014;Kreemer & Zaliapin, 2018). However, the interplay between the hydrological and solid earth systems and its effect on deformation in the LVC area have not been widely investigated. In this study, we describe and analyze both horizontal and vertical displacements from permanent GNSS networks in the context of independent hydrological and seismic data to investigate nontectonic deformation affecting the LVC region and the adjacent SNR. In the study of nontectonic transient deformation, observations of horizontal motion are particularly helpful in determining the spatial distribution of the load source and in identifying different processes causing hydrologically related deformation (e.g., poroelastic deformation) (Devoti et al., 2015;Larochelle et al., 2018;Wahr et al., 2013). Here, we identify anomalous transient horizontal deformation in the GNSS data southwest of LVC and determine that these displacements reflect the solid earth response to two different processes: regional hydrological load changes and local water infiltration into mountain slopes.

GNSS Data Sets and Observations
The continuous GNSS stations in our study area are operated by the U.S. Geological Survey (USGS) and the National Science Foundation's Network of the Americas (NOTA, formerly the Plate Boundary Observatory, PBO). Given the availability of a higher number of sites in the LVC area, we used publicly available daily positions processed by USGS (https://earthquake.usgs.gov/monitoring/gps/stations, Murray & Svarc, 2017) in the North America fixed reference frame. Since hydrological signals span a broad part of the network and region, we did not use the regionally filtered solutions produced by USGS because they may remove some of the signal of interest.
We selected 70 stations spanning the Central to Southern SNR (latitude 36 • N to 38.5 • N, longitude −120 • W to −117 • W, elevation 80 to 3,900 m above sea level, Figure 1 and Table S1 in the supporting information). While most of these stations started operation between 2005 and 2006, ∼30% date back to year ∼2000. We considered the time interval between 1 January 2000 and 31 October 2018 and we excluded seven stations (CASA, KRAK, JNPR, LCOV, P628, P723, and VINE), which lack data during this time span or contain substantial data gaps.
These GNSS stations record long-term trends related to interseismic plate motion, along with significant multiyear transients that have been attributed to LVC inflation activity (Liu et al., 2011;Montgomery-Brown et al., 2015) and variable hydrological loading, the latter of which impacts vertical displacements in particular (Argus et al., 2017). We detrended all time series by removing a steady-state velocity, which we estimated from position data spanning from 2005 to 2011 using robust least squares linear regression. We selected this time interval because of its relative stability in terms of both tectonic/volcanic (no strong inflation/deflation) and nontectonic (modest hydrological loading) deformation. Figures 2a, S2a, and S3a show examples of the so-obtained detrended horizontal and vertical time series (hereafter "multiyear time series"). For the horizontal component we estimated the main displacement direction for each GNSS site and projected the north and east time series along this direction (the direction angle, counterclockwise from horizontal, is indicated in brackets next to the site names in the figures and in Table S1 for all GNSS stations).
We then separated the seasonal deformation signal from the multiyear signal by applying a high-pass filter to each time series. We used a 2.5 year width Gaussian filter, whose weights are given by the Gaussian function and whose width is 6 times the conventional Gaussian sigma (filter1d function of the Generic Mapping Tools (GMT) software Wessel et al., 2013). The 2.5 year filtered GNSS time series are shown in Figures 2a, S2a, and S3a as black dashed lines superimposed to the daily position data.
We afterward used the variational Bayesian Independent Component Analysis (vbICA) decomposition technique by Choudrey and Roberts (2003), as modified for GNSS time series by Gualandi et al. (2016), to identify and remove the horizontal common-mode signal we noticed in the seasonal displacement data. Further details about GNSS data postprocessing are provided in section S1 of the supporting information.  Table S1 for all GNSS stations).
The GNSS time series show transient signals of variable amplitude over both multiyear and seasonal time scales (Figures 2, S2, and S3). To quantify the magnitude of these transients, we estimated the maximum peak-to-peak amplitudes at each site between 2010 and 2016 (time interval chosen for uniformity with loading models, refer to section 4.1), from the smoothed time series obtained with 2.5 and 0.2 year width Gaussian filters, respectively, for multiyear and seasonal time-scales (black dashed lines in Figures 2, S2, and S3). The horizontal amplitudes are calculated from the relative north (N) and east (E) amplitudes as Figure 3 displays the spatial distribution of the maximum amplitude and the main direction of horizontal motion at both the multiyear and seasonal scales.  (Figure 3b), vertical amplitudes are largest for stations on the SNR and decay with distance from the mountain range (see also the cross-section representation perpendicular to the SNR in Figure 4b, red circles). Horizontal amplitudes are generally only 20-40% of the vertical; however, several stations on the southwestern LVC rim show remarkably high horizontal amplitudes (∼10-15 mm), reaching 80-120% of their corresponding verticals. The seasonal horizontal motion of these high-amplitude sites is roughly perpendicular to the SNR, and it follows the change in the strike of the eastern SNR edge. In order to better understand the temporal trend of the seasonal motion, in Figure 4a we display the displacement distribution at two times of year 2011: late-winter/early-spring, when maximum precipitation is usually reached (hereafter designated as "snowfall" and indicated by red arrows); and late-summer/early-autumn, when maximum snowmelt and runoff is usually reached (hereafter "runoff" and orange arrows) (see also section 2.2). Seasonal displacements generally alternate between vertical-subsidence and horizontal-inward movements toward the SNR in late-winter/early-spring ("snowfall") and vertical-uplift and horizontal-outward movement in late-summer/early-autumn ("runoff"). The horizontal deformation direction generally follows the eastern edge of SNR. Except for station MINS, which is the only site on the western side of the SNR drainage divide, all stations on the southern rim of LVC move northeast during late-summer/early-autumn ("runoff").

Hydrological Data Sets and Observations
The SNR is characterized by a high-mountain Mediterranean climate, exhibiting large changes in seasonal temperature and precipitation. During the wet season, precipitation falls as rain and snow, respectively, at lower (below 1,000 m) and high elevations, with an annual average of 1,000 mm/year (1950-2005 interval). During winter, the mean daily temperature is <0 • C, increasing to 15-20 • C in summer. Consequently, runoff 10.1029/2020JB019495 is highly seasonal, with most of the winter precipitation transitioning to snowmelt-driven streamflow in summer (Ficklin et al., 2012(Ficklin et al., , 2013. The SNR and LVC are rich with continuous hydrological stations. We examined various observables (temperature, rainfall, snow water equivalent or SWE, stream discharge, and lake level) at sites within our study region, selected for their record completeness. Site locations, time series, and data sources are displayed in Figures 1, 2, S2, S3, and Table S2, and are summarized below.
Stream discharge data can be used as a proxy for groundwater recharge (Saar & Manga, 2003). Daily discharge data at Hot Creek Flume (HC) and the San Joaquin River (SJ) show seasonal maximum discharge in late-spring/early-summer, episodic discharge peaks immediately following high-precipitation events in winter, and multiyear fluctuations in response to successive years of decreasing/increasing precipitation.
Mono Lake (MN), at the eastern edge of the SNR, receives ∼82% of its annual inflow as runoff (mostly from SNR snowmelt) and 18% from precipitation (Ficklin et al., 2013). Seasonal and multiyear lake level fluctuations (up to 2 m) are clearly associated with precipitation fluctuations.
Temperature data at Mammoth Lakes (Figures S2 and S3) show seasonal peak-to-peak variations of up to ∼30 • C, but only small differences (few degrees) between wet and dry years. GNSS stations in the LVC region observe the superimposition of tectonic, volcanic, and hydrological processes. We focus the remainder of this paper on the seasonal time scale, since we have more information about seasonal dynamics driven by hydrology. We acknowledge, however, that there are intriguing phenomena occurring at longer periods, such as the possible temporal correlation between inflation episodes at LVC and periods of severe drought across the Western United States (Hammond et al., 2019), that warrant additional study.

Comparison Between GNSS and Hydrological Data
Vertical seasonal GNSS displacements are generally in-phase with precipitation patterns (subsidence in winter/spring months, uplift in late-summer/early-autumn), and their amplitudes are noticeably higher in years with higher precipitation (Figures 2b, S2b, and S3b). Horizontal seasonal displacements are ≤5-6 mm within our study area, except for the subset of stations located on the southwestern caldera rim (Figure 3b). Horizontal displacements from these stations noticeably mimic the shape of hydrographs, with a fast-rising limb peaking in late-spring/early-summer, followed by slower recession. The timing of maximum horizontal displacement is usually slightly delayed with respect to the spring discharge peak and advanced with respect to the vertical uplift peak (Figure 2b).
To examine the unusual displacement signal on the southwestern LVC rim in more detail, we selected seven representative GNSS stations with significant horizontal deformation (LINC, MINS, P630, P631, P642, SHRC, and TILC). Figure 5 shows the comparison of the horizontal components projected along the directions of maximum motion (refer to section 2.1), the vertical seasonal time series and hydrological data, focusing on three 3-years periods of heavy precipitation. Cross-correlation analysis (applied to 0.2-year width filtered time series-black dashed lines in Figure 5) between horizontal GNSS displacements and Hot Creek discharge time series (Table S3) results in high positive correlation coefficients (∼0.74) with positive movements (outward from SNR and towards LVC) delayed by ∼30 days (LINC, MINS, and P630) to ∼55 days (P631, P642, SHRC, and TILC) with respect to spring discharge. Lower and negative correlation values (∼ −0.60) result for the vertical components, with vertical subsidence preceding spring discharge by ∼80 days. We also cross-correlated GNSS displacements with time series of SWE over the SNR from Margulis et al. (2016). Vertical subsidence is well correlated with SWE (∼0.73) and slightly delayed (∼20 days) with respect to peak SWE, while positive horizontal displacement exhibits high correlation (∼0.77) and a delay of ∼5 months. Figure 6. Yearly time occurrence of "runoff" (a) and "snowfall" (b) peak epochs estimated from GNSS time series (red and orange lines) and hydrological and seismicity data. The latter are indicated with gray stars color coded by the peak prominence. The blue lines represent discharge at Hot Creek (dark blue) and San Joaquin (light blue). The gray line in (b) indicates SWE peak epochs from Margulis et al. (2016). Open symbols on the left of each data set represent (weighted) mean and standard deviations of the peak epochs over the entire period. Subscripts "r" and "s" refer respectively to runoff and snowfall phases. The table (c) summarizes the time difference between the means of the various data sets. Margulis et al. (2016) showed that the peak timing of SWE over the SNR varies significantly between years, ranging from mid-January to mid-May (here reproduced in Figure 6b, gray line). For each year between 2005 and 2018, we calculated the epoch of "runoff" and "snowfall" peaks estimated from the corresponding peaks in GNSS horizontal and vertical time series (specifically, the mean and standard deviation of peak times for each station) and the epoch of peak spring discharge for both Hot Creek Flume and the San Joaquin River. Discharge peaks generally fall between May and July and occur earlier in drought years than in high-precipitation years (blue lines in Figure 6a). GNSS horizontal "runoff" peaks (red line in Figure 6a) are correlated with interannual precipitation (delayed peaks in wet years), and exhibit a persistent temporal delay of ∼45 days with respect to annual peak spring discharge, in agreement with the cross-correlation results we obtained above using the entire time series. "Runoff" vertical displacement peaks (orange line in Figure 6a), while noisy, generally occur later (∼4 months) than the horizontal peaks. On the contrary, "snowfall" horizontal and vertical peaks occur at similar times in late winter and early spring, roughly in agreement with the SWE peak-time estimates from Margulis et al. (2016) (Figure 6b). The substantial temporal variability between snowfall and subsequent runoff, which can differ by up to ∼5 months annually, may partially explain why studies of seasonality in seismicity based on calendar time (Christiansen et al., 2005) find only small enhancements during calendar-defined winters.

Discussion
Different processes can contribute to the seasonal signal in GNSS time series (Xu et al., 2017, and references therein), such as surface loading due to mass redistribution, thermoelastic deformation from changing surface temperature, poroelastic deformation caused by subsurface hydrological variation, and instrumental and processing errors. Here we examine different possible causes of the observed signals.

Response to Surface Loading
As also observed at larger scales by Kreemer and Zaliapin (2018), the GNSS seasonal vertical subsidence/uplift and horizontal inward/outward motion around the SNR (Figure 4a) suggests that the solid earth elastic response to snow and water loading plays a significant role in driving the observed deformation. We demonstrated this by estimating the expected displacement response in our study area from hydrological loads at both regional scale (i.e., at the scale of Sierra Nevada) and large scale (i.e., at the scale of the contiguous USA, or CONUS). In both cases, we forward modelled the corresponding displacements at each GNSS station in our study using the SPOTL program (Agnew, 2012), assuming a spherical Earth with elastic rheology provided by the Gutenberg-Bullen Model A average-Earth model and displacements referenced to Earth's center of mass.
For our snow load data set over the SNR, we used the newly developed high-resolution SWE reanalysis from Margulis et al. (2016) ("Margulis SWE"), which is based on the assimilation of remotely sensed fractional snow-covered area data over the Landsat 5-8 record. The data set features daily temporal and 90 m spatial resolution, is updated through winter 2015/2016, and is available at the website https://margulis-group. github.io/data/. Since Margulis SWE data are only available through 2016 and because most of the GNSS time series start after 2007, we only considered the period between the wet winter years of 2010 and 2016 for our analysis. For our CONUS-scale data set, we used daily soil moisture (0-2 m depth) and SWE from NASA's North American Land Data Assimilation System (NLDAS) 0.125 • -gridded NOAH model.
We estimated the displacements associated with the load changes over the CONUS region (large-scale) and over just the SNR (regional-scale). For regional-scale loads, we used the sum of NLDAS soil moisture data and Margulis SWE. For large-scale loads, we used the sum of NLDAS soil moisture over all of CONUS, Margulis SWE over the SNR only, and NLDAS SWE for the rest of CONUS. The resulting seasonal time series for 2010-2016, obtained analogously to the observed GNSS time series (i.e., removing a 2.5-year-width Gaussian filter, refer to section 2.1 and step 3 in Text S1), are shown in Figures S4 to S6 as green (regional-scale) and blue (large-scale) lines.
The load model results are displayed in Figure 4. Figure 4a shows the displacements due to large-scale hydrological loading at the "snowfall" peak in April 2011 (blue arrows) superimposed to the corresponding observed displacements (red arrows). Figure 4b represents a cross-section view of the seasonal peak-to-peak amplitudes for the observed (red circles) and the modeled displacements (green and blue circles, respectively, for regional-scale and large-scale loads).
Focusing first on the vertical component, the peak-to-peak amplitudes in Figure 4b (right panel) show that regional-scale loads (green circles) are insufficient to explain observed displacements (red circles). However, considering loads over all of CONUS (blue circles) results in a greatly improved fit, indicating that both near-field and far-field hydrological loading drive the observed vertical displacements. The remaining disagreement with observations, particularly at high elevations, could be ascribed to the underestimation of the subsurface storage in the SNR, not fully characterized by NLDAS soil moisture data, as discussed in Argus et al. (2017) and Enzminger et al. (2019).
By way of contrast, Figure 4b (left panel) shows that neither near-field nor far-field loading results in significant horizontal displacements. Even if part of this disagreement could be also associated with the

10.1029/2020JB019495
underestimation of the subsurface storage, the large observed horizontal signals at the caldera southwestern rim are unlikely to be explained by surface loading. This is both because the horizontal amplitudes are too large (the expected horizontal response to loading is generally >2 times smaller than the vertical one (Wahr et al., 2013)) and because the shape and timing of the horizontal behavior does not align with that of the known SNR snow/water loading (refer to section 3). Furthermore, examining the residuals between the amplitudes of observed and modeled ("large-scale load") displacements (Figure 7a), we observe a strong (>10 mm) south-north variation across the south-western LVC in the horizontal component, which occurs over spatial scales of ∼10 km (Figure 7b). These horizontal residuals could potentially be related to the effects of lateral heterogeneities in crustal elastic parameters or to the presence of compliant zones related to extensive faulting in LVC (Flinders et al., 2018;McMechan et al., 1985), which are not considered in our laterally homogeneous model. These effects, however, would also be expected to impact the vertical residuals, which instead show no clear spatial pattern.
To be thorough, we also estimated the load effects of the water volume changes of some lakes and reservoirs located close to the GNSS sites (Cherry Lake, Hetch Hetchy Reservoir, Mono Lake, and Crowley Lake-data from https://maps.waterdata.usgs.gov/ and https://cdec.water.ca.gov/cdecstation2/?sta=BIS). The resulting deformation is however lower than the GNSS data noise level, with a maximum amplitude of ∼1.5 mm for the vertical component of the site DECH, which is located very close (∼3 km) to Mono Lake.
Given the evidence, we argue that one or more physical processes other than loading are driving observed horizontal deformation in the LVC. We examine two of the possibilities below.

Response to Temperature Variations
Temperature fluctuations of the Earth's surface may induce strain and tilting due to the thermal expansion/contraction of surface rocks. In areas with large topographic relief, thermoelastic deformation due to surface temperature variations can be significant (Ben-Zion & Allam, 2013;Prawirodirdjo et al., 2006). We estimated this contribution to the observed GNSS displacements using the model from Tsai (2011). Considering an elastic half-space covered by a thin unconsolidated layer, Tsai (2011) computed the equations for the displacement induced by a sinusoidal (in time and in space) surface forcing associated with temperature variations. We used these equations ((11a) and (11b) in Tsai, 2011) to estimate the expected horizontal and vertical surface displacements due to assigned surface temperature variations. We assumed peak-to-peak temperature variations ΔT = 30 • C (from Mammoth Lakes temperature data, Figures S2 and S3); a range of wavelengths for the spatial temperature variation (1, 5, 10, and 50 km) to take into account the short-and long-range features of local topography; annual periodicity (10 −7 s −1 ); varying values of upper layer thickness (0-10 m); and typical values of the Poisson's ratio (0.25), the coefficient of linear thermal expansion (10 −5 • C −1 ), and the thermal diffusivity coefficient (10 −6 m 2 /s) (Kraner et al., 2018;Tsai, 2011). The resulting maximum peak-to-peak displacement amplitude amounts to ∼2 and 1 mm, respectively, for the horizontal and vertical components.
The magnitude of the thermoelastic vertical deformation could be increased by up to 3 times at the near surface if the GNSS stations are anchored within the thermal boundary layer (Tsai, 2011). Most (66%) of the GNSS sites used in this paper are drill-braced monuments, which perform 30% to 50% better than other monument types on the basis of rate errors and seasonal noise (Langbein & Svarc, 2019). Specifically, 41% are short-drilled (up to 2 m) braced monuments (SDBM), usually used where competent rock is at the surface, and 27% sites are deep-drilled (up to 12 m) braced monuments (DDBM). Thirteen percent of the GNSS sites are shallow foundation masts (∼0.5 m), while the remaining are other (e.g., bedrock-bolted mast) or unspecified kinds of monuments. Monuments details of GNSS sites can be found in the stations log-files (e.g., https://earthquake.usgs.gov/monitoring/gps/LongValley/p642/logs).
Our results show that estimated thermoelastic deformation (∼1-3 mm) represents only a small fraction of the observed ∼10-15 mm amplitude of seasonal horizontal displacements at GNSS stations along the southwestern rim of LVC. The Tsai (2011) model does not consider lateral variations in Earth's elastic properties that, at local scales, could cause millimeter-level horizontal deformation (Fleitout & Chanard, 2018). A more realistic model would therefore increase the thermoelastic deformation at sites on the caldera rim, but we expect that it still would not entirely explain the significant magnitude, interannual variation, asymmetric shape, and peak timing of the observed horizontal signal (occurring before the temperature peak in summer).

Response to Pore Pressure Variations
In several regions characterized by high precipitation and fractured-rock outcrops, water infiltration and related subsurface pore pressure changes have been suggested as the cause of considerable (several millimeters) local-scale surface displacements, especially in the horizontal direction (Hansmann et al., 2012;Oestreicher, 2018;Serpelloni et al., 2018;Silverii et al., 2019). Water infiltration in the SNR during winter is limited due to its accumulation in the solid state as snow. Snowmelt in spring, however, produces large volumes of mobile water that can fill and expand fracture voids. Rapidly rising water levels in these fractures simultaneously increase pore-fluid pressure and cause an overall rock volumetric expansion, while subsequent draining results in longer-term contraction. This mechanism could be responsible for the observed alternating pattern of rapid outward and slower inward motion at GNSS stations on the southwestern rim of LVC. In this interpretation, the amplitude of the horizontal deformation would correlate with the amplitude of seasonal groundwater-level fluctuations (<2 m in dry years and >10 m in wet years, Montgomery-Brown et al., 2019) and the deformation would preferentially occur in typically highly fractured areas around the caldera (Acocella, 2010;Walter & Motagh, 2014;Wang et al., 2017). Interestingly, the largest horizontal amplitude occurs at station MINS, which moves in southwest direction in spring/summer and lies on the west side of Mammoth Mountain, one of the main infiltration areas of LVC hydrothermal system (Pribnow et al., 2003;Peacock et al., 2016).
Recently, Montgomery-Brown et al. (2019) showed that shallow seismicity (<4 km below surface) in the SNR south of LVC is strongly modulated by snowmelt: pressure diffusion from water infiltration accelerates seismicity rates, causing nonvolcanic earthquake swarms along steeply dipping strata that act as high-permeability pathways and faulting planes. In particular, the analysis of 33 years of seismicity in the Convict Lake area revealed a time lag of 24 days between peak streamflow and peak earthquake rate, suggesting a rapid downward propagation of the pore pressure perturbation. The location of Convict Lake is indicated in Figure 1, while seismic data (weekly earthquake counts) are displayed in Figure 5. We refer to Figure 3A in Montgomery-Brown et al. (2019) for a representation of the entire seismicity time series.
As previously described (refer to section 3 and Figures 5 and 6), our data show a delay between horizontal deformation on the southwestern LVC rim and spring discharge, particularly at stations located at higher elevations in the SNR (LINC, MINS, P630, P631, and P642). We analogously calculated the peak epochs of Convict Lake seismic swarms. Since the swarm was not always activated, also in some years when the runoff was still quite high (Montgomery-Brown et al., 2019), we show results for the years with a fairly pronounced peak in seismicity. In particular, we calculated the peak prominence (measure of how much the peak stands out due to its intrinsic height and its location relative to other peaks, as in the findpeaks MATLAB function), and we selected the years with a peak prominence >5 earthquakes (color-coded stars in Figure 6a). As also visible from the seismicity time series in Figure 5, 2017 is the year with the most prominent peak, followed by 2005 and 2006. Resulting peak epochs show a delay of ∼30 days relative to spring discharge. This is comparable to the time lag between horizontal deformation and spring discharge, especially for LINC, MINS, and P630 (Table S3), suggesting that pore pressure diffusion in fracture networks (associated with superficial water infiltration) is the driver of both the shallow seismic swarms and the (horizontal) deformation in this area.
To precisely model the observed horizontal deformation around LVC, we would need detailed information about the spatial distribution and the mechanical properties of local rock fractures (not accessible at our spatial scales) and diffusive (anisotropic) subsurface water flow. Here we simply simulated the overall response to pore pressure variations using an analytic model for axisymmetric reservoirs in an isotropic, poroelastic half-space. This model, from Segall (1992), is a generalized solution giving displacements and stresses generated in rocks surrounding a reservoir undergoing contraction or expansion associated with a known pressure distribution (Geertsma, 1973). The model considers an axisymmetric reservoir (defined by its radius r and thickness z) embedded in an isotropic poroelastic half-space, which requires specifying the bulk modulus (K), Poisson's ratio ( ), and effective stress coefficient ( ) of the medium, and assigning a pore pressure variation inside the reservoir (Δp).
In applying this model we assumed a uniform, constant pressure increase in the reservoir to simulate the amplitude and the spatial distribution of the residual horizontal seasonal displacement amplitudes at their maximum (Figure 7). We projected the residual seasonal amplitudes in Figure 7a on a profile perpendicular to the LVC south rim (red circles in Figure 7b) and compare them with the deformation produced by the model. We used a source domain located on the SNR rim, south of the GNSS stations, tangent to the cross-section line (black straight line in Figure 7a) and extending from the surface to a depth equal to the source thickness (z). Figures 7c.1 and 7c.2 depict a sketch of the model setup. Figure 7c.3 shows the comparison between observations and results for a model with fixed radius and thickness (r = 6 km, z = 7 km). The figure also shows two other results for reference (z∕3, dashed line, and 3z, dotted line) to show the sensitivity of the results to the source thickness: a larger source thickness does not match the profile of observations well at small radial distances, while a more shallow source would require a much larger (∼50-100%) pore pressure increase to match the observations.
Fixing the values of Poisson's ratio ( = 0.25) and effective stress coefficient ( = 1), we evaluated the effect of the bulk modulus and pore pressure change. We varied K between 10 and 60 GPa-the latter being compatible with the Gutenberg-Bullen Model A average Earth model used for the load model in section 4.1-and we varied Δp between 100 and 150 kPa, corresponding to upper-limit values of water level fluctuations in this area (Montgomery-Brown et al., 2019). For each K-Δp combination, we inverted for the source geometry (radius r and thickness z) that best fits the observations. We calculated reduced 2 misfit over a range of radii and thicknesses. We found a strong trade-off between source radius and thickness, producing a large range of statistically equivalent models-that is, ones with a reduced Chi square misfit of 2 ≤ 1. We therefore randomly sampled (using the uniform distribution) radius/thickness pairs falling within the 2 < 1 contour, and forward calculated the poroelastic model prediction. Afterward, we calculated the 99%, 75%, and 50% (median) quantiles of the predicted displacement distribution as a function of the distance from the source. The results are depicted in Figure S7: 2 values and random sampling within the 2 = 1 contour on the left, and predicted distributions on the right. The predicted displacements obtained for K = 60 GPa and Δp = 150 kPa provide a good agreement with the observations amplitude and spatial distribution. The distribution of predicted displacements in this case is represented as gray areas and black line in Figure 7b. It is important to note that, since we chose a value of 1 for the Biot coefficient ( ) to represent the most sensitive fluid-rock system possible, this value of pressure change is likely to represent a lower bound. If were ∼0.4, as would be the case with unfractured crystalline rock (Wang, 2000), the pressure change would need to be 2.5 times larger.
Topographic features could also affect amplitude and direction of infiltration-induced displacements at stations located on steep slopes edges (Devoti et al., 2016;Serpelloni et al., 2018). In particular, we expect that the topography aspect angle (the direction that a slope faces) could influence the motion direction, while the slope angle (degree of inclination of the slope relative to the horizontal plane) could affect the motion amplitude. In order to evaluate this possible effect, we implemented a simple analysis based on the estimation of the topographic gradient vector (defined by the aspect and slope angles, as in Materna, 2014), using increasingly coarser topography data, obtained by spatially filtering the high-resolution Shuttle Radar Topography Mission digital elevation model (DEM) data (https://www2.jpl.nasa.gov/srtm). Further details of this analysis are described in section S2 and related Figures S8 and S9 in the supporting information. Our results show that the orientation of the observed horizontal seasonal displacements, which are roughly perpendicular to the eastern boundary of SNR, correlates with topographic aspect at scales of several km (∼15-40 km), while it is minimally influenced by smaller-scale topographic features. This suggests that the forcing mechanism for observed horizontal deformation is related to features varying at the scale and with the shape of the SNR itself, as the spatial distribution of snow/groundwater and of rock fractures. The volumetric changes within fractures will be aligned with the direction of smallest confining pressure (here, the free surfaces associated with the eastern SNR slopes).
Given the compatibility between model results and observed GNSS displacements, we infer that poroelastic effects are able to generate the magnitude and spatial distribution of the seasonal horizontal displacement in LVC and environs. Since the source dimensions are difficult to constrain with available data, additional GNSS observations to the southwest of the cross-section line in Figure 7a would be beneficial. Incorporating material anisotropy and heterogeneity and layered structure in the model is also possibly needed to fully account for the observed horizontal deformation.
As an alternative to approximating the overall behavior of a dense network of fractures, previous works have implemented linear fault models using a single, equivalent tensile dislocation embedded in an elastic half-space (Devoti et al., 2015;Silverii et al., 2016;Serpelloni et al., 2018). Instead of a linear fault model 10.1029/2020JB019495 expanding at unrealistic depths, we prefer the spatially extensive poroelastic model which permits linking the pore pressure variation over a certain area (ideally the one covered by distributed fractures) to the observed displacements. The agreement between the displacement direction and topographic features at scales of several kilometers provides support for a model based on an extended source.
Because the process responsible for the observed surface displacements must be reversible in order to generate the observed cyclic behavior, we can also exclude the possibility that the observed deformation might be due to slip on shear dislocations (such as on the fault related to the Convict Lake seismic swarms) facilitated by normal stress reduction due to increased pore pressure. Slip on shear dislocations is usually not a reversible process, unlike the opening of a tensile fracture, which would be able to expand and contract to produce the observed cyclical displacements due to variations of stress associated with changing pore-pressure.

Conclusions
Combining displacement data from permanent GNSS networks with various hydrological records, we analyzed the nontectonic deformation affecting Long Valley Caldera and the adjacent Sierra Nevada range and highlight strong correlations between surface displacements and hydrology at both seasonal and multiyear scales. We observed a clear deformation pattern with spatiotemporal variability that is likely related to surface loading acting at large (CONUS) to regional (Sierra Nevada range) scales. We highlighted anomalous seasonal horizontal deformation not explainable by the solid earth elastic response to hydrological loading. This deformation is most apparent at GNSS sites along the eastern SNR front and on the southwestern rim of LVC. The timing is strongly correlated with the discharge profile from local streams. In this area, pressure diffusion from water infiltration into subvertical rock fractures dramatically accelerates shallow seismicity rates. We suggest that the anomalously large observed GNSS horizontal displacements are similarly related to pore pressure changes due to groundwater level increases and subsequent draining within fractured rock.
This study shows that GNSS networks can provide detailed information about groundwater features on local (km) scales. An enhanced understanding of the causes and implications of the deformation at LVC is of primary interest for a better comprehension of possible interactions between tectonic/volcanic processes and hydrological processes such as aquifer charging and draining. Increased awareness of the effects of nontectonic processes in LVC area would also improve the observation and study of the volcanic system itself.
In this paper we focused mainly on seasonal deformation and simple modeling. In future work, we plan to implement more realistic 3-D models to evaluate the effects of heterogeneity and topography on measured deformation and the load variations impact on the LVC volcanic system.