Exploring the InSAR signature associated with river-sourced recharge in California's San Joaquin Valley

InSAR time series of vertical displacement over our study region (figure 2(b)), were generated by TRE Altamira and made publicly available by the California Department of Water Resources [27]. We used the vertical displacement time series provided as input into our InSAR time series decomposition (section 4.1). While a full description of processing and calibration/validation workflows are outlined in supporting data documentation [27], we provide a summary here. Data spanning 1 January 2015–1 January 2023 from the European Space Agency's Sentinel-1 missions was processed using the SqueeSAR method [28] over groundwater basins in California. The resulting line-of-sight (LOS) displacement estimates were then calibrated using a network of 536 continuous Global Navigation Satellite Systems (GNSS) stations. Projecting the 3D GNSS measurements into the satellite 1-D LOS, the average displacement velocity (linear trend) of each station was compared against the displacement velocity of the average time series determined from InSAR pixels within a 100 m radius. At regional scales, pairwise differences between GNSS and InSAR velocities were used to estimate and remove a first order plane from all InSAR pixels. InSAR time series were then shifted by a common time series of residuals, with respect to the GNSS LOS time series, to force the relative InSAR measurements into the absolute reference frame of the GNSS network. This has the additional benefit of minimising non-geophysical signals in the time series, such as tropospheric delays or changes in scattering properties [29, 30]. To create a consistent dataset from multiple orbital geometries and tracks, full resolution InSAR LOS results were subsampled to a common spatial grid (100 × 100 m). Similarly, the time series for each pixel on the common spatial grid are interpolated and resampled to a common temporal grid with an approximate 7 d interval, resulting in 481 measurement times. Where multiple orbital geometries (ascending and descending) were available, these re-gridded LOS results were combined to estimate the component of vertical displacement. Otherwise, the LOS results were projected into vertical direction using the satellite look angle.

Using an independent set of 181 continuous GNSS stations, the Root Mean Square Error (RMSE) of time series for each station and calibrated InSAR dataset resulted in a state-wide RMSE value of 8.86 mm with individual stations ranging from 29.37 to 1.33 mm. The quantitative validation analysis determined a vertical accuracy assessment of 18 mm (95% confidence interval) for the InSAR time series [28].

AEM data were previously acquired in the study region using SkyTEM 508, 312 and 304 systems [21–23, 31], and inverted to obtain electrical resistivity models of the subsurface. Shown in figure 2(b) are the flightlines along which the AEM data, corresponding to the resistivity models used in this study, were acquired. While data were collected continuously, measurements were stacked at ∼35 m-intervals (referred to as AEM soundings) along the flightlines.

For all three previous surveys, a spatially constrained inversion was used to obtain models of the subsurface resistivity from the AEM data. This approach prioritised a resistivity model that fit the data while favoring smooth transitions of resistivity values in vertical and lateral directions [22, 32]. A resistivity model is composed of 1D resistivity profiles at each AEM sounding. While measurement errors and inversion decisions introduce uncertainty in the resistivity profiles, the penalisation of resistivity variations not supported by the data during the inversion process provide strong confidence in the large-scale structure imaged [22]. These 1D resistivity profiles consisted of 30–40 vertical cells. The thickness of the cells ranged from 1–3 m at the surface, then increased with depth with a constant geometric factor (e.g. 1.07) to compensate for the degrading resolution of the AEM data [33]. The depth of investigation (DOI), which describes the maximum depth at which there is measurable change in the AEM data due to changes in subsurface resistivity, was on average 300 m below ground surface (mbgs), ranging from 150–400 mbgs, in the three surveys.

To analyse the InSAR data, we adopted the curve fitting approach of Neely et al [20] using linear and sinusoidal terms to approximate the InSAR time series:

$$\begin{align}Y\left( t \right) & = vt + A\cos \left( {2\pi \left( {t - \tau } \right)} \right) + {Y_0}\end{align} \tag{ 1 }$$

where $Y\left( t \right)$ is the modeled displacement time series (mm), $t$ is time (year), $v$ is the linear term estimating the rate of displacement (mm yr⁻¹), $A$ is the peak amplitude (mm) for the seasonal term, $\tau$ is the timing to peak amplitude for the seasonal term (year, where $\tau$ = 0 represents a timing in the peak amplitude term on 1 October, the start of the California water year and suggesting groundwater storage peaks around this time), and ${Y_0}$ is a constant displacement shift (mm). Terms in equation (1), conceptually shown in figure 1, were estimated using least-squares minimisation techniques applied to the vertical displacement time series and computed for each water year. The average to the year-by-year standard deviation to peak amplitude and timing to peak amplitude are shown in figure S1. Examples of model fits to the time series are shown in figures S2–6.

Neely et al [20] showed that these terms differed between a dry year (WY2016) and a wet year (WY2017), with the determination of dry vs. wet based on drought conditions and water availability relative to the average [34, 35]. We applied the same methodology to observe the behavior on a year-by-year basis. While the previous study [20] was limited to two water years, we assessed the consistency of terms for 7 years. These include 5 dry years (WY2016, WY2018, WY2020–2022) and 2 wet years (WY2017 and WY2019).

To relate displacement values to recharge via storage changes e.g. [11], we assumed InSAR is observing primarily poroelastic deformation in our study region [4, 36–38]. As pixels with low peak amplitude values are expected to give higher uncertainties for the timing to peak amplitudes, we excluded $\tau$ values for pixels that had $A$ values less than the accuracy of the InSAR data (18 mm peak-to-trough) [27]. For visualisation purposes, we filled in regions of missing data using a median filter (5 km diameter) on our derived maps of rate of displacement, peak amplitude, and timing to peak amplitude.

We used the available 1D resistivity profiles to create a laterally continuous 3D model of subsurface resistivity for the study region. We first aligned each sounding to a common reference frame by vertically shifting each by their ground elevation. As the 1D resistivity profiles have variable thickness of cells increasing with depth, we then resampled resistivity values at 5 m intervals in the vertical direction with depths below the DOI excluded from the analysis. A nearest-neighbor interpolation scheme with a 12 km search radius was applied to each vertical interval using Generic Mapping Tools [39] to create 400 m × 400 m resolution horizontal grids.

To interpret sediment type from resistivity in our study region, we applied the method presented in Knight et al [21]. Using sediment type reported in driller's logs located within 100 m of a 1D resistivity profile (349 well logs), we determined the distribution of resistivity values that corresponded to coarse-grained dominated sediments and fine-grained dominated sediments (hereafter referred to as coarse- and fine-dominated sediments respectively). We limited this analysis to the saturated zone, using water table levels reported at the time of AEM data acquisition (10–20 mbgs at the Fresno site and 20–40 mbgs at the Visalia site). We presumed that the resistivity distributions in the relatively small region above the water table would broaden (due to the saturation state varying from dry to saturated) and shift to higher resistivity values as shown in a detailed study of the impact of saturation on AEM-derived resistivity values [40].

The InSAR time series decomposition using a curve fitting approach resulted in three parameter terms, (1) the average rate of displacement, $v$, (2) the peak amplitude, $A$, and (3) the timing to peak amplitude, $\tau$, on a year-by-year basis. Although the focus of this study was to understand how InSAR data might capture river-sourced recharge during the wet years, we also analysed the signals observed in the dry years for completeness.

5.1.1. Dry water years

Our curve-fitting approach for the dry water years (WY2016, WY2018, WY2020–2022) displayed similar spatial patterns in the linear rate of displacement, peak amplitude, and timing to peak amplitude with respect to each other dry year. The rate of displacement term featured a large subsidence bowl with surface elevations falling more than 250 mm yr⁻¹ (purple in figure 3 column 1). The magnitude and spatial pattern of subsidence, often attributed to groundwater pumping, was consistent with previous interferometric studies conducted over the San Joaquin Valley during drought [7, 20, 41, 42]. We found that in regions of higher subsidence rates, peak amplitude was generally higher (up to 80 mm peak-to-trough motion and indicated by yellow in figure 3 column 2). The timing to peak amplitude occurred between late winter (pink in figure 3 column 3) and early summer (yellow-orange in figure 3 column 3).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For these dry years, we did not observe the general spatial pattern in the timing—progressing from east to west—that had been interpreted as river-sourced recharge in the study by Neely et al [20]. Rather, the timing to peak amplitude for these years were in-phase with expected seasonal oscillations of groundwater storage which, on average, peaks around April for this region [16, 37, 43]. The expected seasonal oscillations are aligned with local precipitation primarily falling between winter through early spring [16, 43] and with lowest groundwater elevations due to groundwater pumping occurring in autumn. The overall drier conditions and increased reliance on groundwater, as evidenced by decreases in groundwater elevations [25], over the study site during these dry years limited our ability to assess InSAR signatures related to river-sourced recharge.

5.1.2. Wet water years

Given our hypothesis, we focused our analysis on the wet water years, when we expected to see river-sourced recharge. The wet water years (WY2017 and WY2019), when compared to each other and to the results in Neely et al [20] for WY2017, had consistent spatial patterns in the linear rate of displacement, peak amplitude, and timing to peak amplitude. The rate of displacement for both years showed broad subsidence with values up to 150 mm yr⁻¹ (dark orange in figure 3 column 1). While these wet years exhibited increases in groundwater levels (upwards of 20 m for both study sites) [25], subsidence due to the residual compaction of clay, can occur on decadal to centurial time scales [44]. We found peak amplitude values up to 70 mm peak-to-trough (indicated by yellow-orange in figure 3 column 2) at various locations, including the Visalia site. In locations where peak amplitude values were greater than the reported accuracy of the InSAR displacement time series (18 mm), we had confidence in our ability to quantify timing to peak amplitude—the critical feature we wanted to isolate in the InSAR data. Further providing confidence in our analysis, the average standard deviation for the timing to peak amplitude term, shown in Figure S1b, was less than 30 d for our Visalia site. As the peak amplitude values at our Fresno site were less than 18 mm, there was effectively no InSAR signature at this site that could be considered as related to recharge.

The timing to peak amplitude at the Visalia site exhibited a distinct spatial pattern in both WY2017 and WY2019. Near the eastern edge of the valley, peak amplitudes occurred in late winter (pinks and reds in figure 3 column 3). Moving towards the valley center, peak amplitude occurred at progressively later times in the year (orange—blue in figure 3 column 3). Taking a transect across this feature, the progression of the timing to peak amplitude is exemplified in figure 4 with the time series and models shown in figures S2–6. This InSAR signature was focused around Cross Creek and Mills Creek; on the order of 10 s km in scale. This differs from the broad, regional behavior in the timing to peak amplitude observed in the dry years. Although the most prominent example of this spatial pattern in timing to peak amplitude was at the Visalia site, this pattern was present elsewhere in the study region; for example along Tule River and Deer Creek. Due to the absence of high density AEM resistivity data, we did not assess these other locations of interest.

Zoom In Zoom Out Reset image size

Given the wet conditions [34, 35], locally focused seasonal deformation patterns, proximity to rivers and streams, and evidence of increased groundwater levels in the area [25], we interpreted the InSAR signature near Visalia to be the expression of river-sourced recharge.

The AEM resistivity values found at the study sites ranged from 3–1000 Ohm·m. The highest resistivity values (300–1000 Ohm·m) in the region were interpreted as the bedrock of the Sierra Nevada, as suggested by earlier studies [21, 23].

Using the relationship derived with the methodology presented in section 4.3, we found the distribution of resistivity values corresponding to fine-dominated sediments ranged from 11–18 Ohm·m while for coarse-dominated sediments, the distribution ranged from 25–50 Ohm·m, as shown in figure 5. Resistivity values below 18 Ohm·m were interpreted as clay-rich sediments, as clay (rather than silt) was the primary sediment type within the fine-dominated classification of the driller's logs; resistivity values above 25 Ohm·m were interpreted as sand and gravel. Resistivity values between 18 and 25 Ohm·m were interpreted as mixtures of clay, sand, and gravel, with the amount of sand and gravel increasing as resistivity increases.

Zoom In Zoom Out Reset image size

At the Fresno site, the presence of the perennial King's River to the east is a known source for recharge [26]. Increases in local groundwater elevations during wet years [25] and previous isotope studies [18] supported our assumption that river-source recharge occurred. When looking at the resistivity structure of this study site, shown in figure 6(a), we interpreted three key features. On the eastern edge of the valley, the high resistivity values (red in colour) corresponded to the granitic bedrock of the Sierra Nevada. Moving into the valley, the base of the resistivity model showed a low resistivity package (blue—green in colour), interpreted as clay-rich sediments, with the top of the package reaching near the ground surface at the eastern edge and dipping to the greater depths to the west. This clay-rich sediment package would suggest low permeability and be unlikely to provide a recharge pathway into the deeper aquifer. In the remainder of the displayed section there is a package, extending from the ground surface to the top of the clay-rich sediments, dominated by the resistivity values corresponding to sand and gravel (orange—red in colour) which likely provide pathways for recharge. The absence of evidence for significant clay within this package is consistent with reported unconfined conditions in the aquifer system across the entirety of the Kings subbasin [26].

Zoom In Zoom Out Reset image size

While there is evidence of river-sourced recharge in the Fresno site [18, 25], the low clay content associated with the permeable pathways limited the potential for poroelastic deformation to occur or be observed. This, combined with the absence of an InSAR signature, supported our hypothesis that clay is essential for observing an InSAR signature.

At the Visalia site, where an InSAR signature interpreted as being related to river-sourced recharge is present, there are the ephemeral Cross Creek and Mills Creek that connect to the Kaweah River system during wet years. As with the Fresno site, there is a distinct isotopic signature attributed to river-sourced recharge and local groundwater elevations rose during WY2017 and WY2019, supporting assumptions that river-sourced recharge occurred [18, 25]. The underlying resistivity model shown in figure 6(b) has a more complex structure than the Fresno site. Along the eastern edge, we observe high resistivity values corresponding to the Sierra Nevada bedrock (red in colour). This is juxtaposed by a low resistivity package (blue—green in colour) interpreted as clay-rich sediments. The top section of the resistivity model, overlying the clay-rich sediments, has resistivity values primarily in the range of ∼18–35 Ohm·m (yellow—orange in colour) with thickness of <10 m to the east and thicken to the west. In contrast with the interpreted sand and gravel package near the surface of the Fresno site, the top sections of the model at the Visalia site displayed lower resistivity values. We interpreted this as sand and gravel with relativity high amounts of interbedded clay. Intersecting this package in the west, and underlying the InSAR signature, is a ∼30 m thick lower resistivity unit (blue—green in colour) dipping the southwest with depth ranging from 50 to 150 m. We interpreted this unit as a clay-rich confining layer, likely corresponding to the Corcoran Clay.

We interpreted, shown on figure 6(b), the sand and gravel as recharge pathways for river-sourced water to migrate from the surface to greater depths in the aquifer system. The high clay content from interbedded clays and confining layer creates the semi-confined and confined conditions needed for observable poroelastic deformation responses. As such, an InSAR signature is observed during times of recharge. Our hypothesis of the essential role high clay content has in determining whether or not we can observe an InSAR signature due to river-sourced recharge is supported by the underlying hydrogeological structure and sediment type in the Visalia site.

This study explored the InSAR signature associated with river-sourced recharge in California's San Joaquin Valley. Interpretations of hydrogeologic structure and sediment type from AEM resistivity data supported claims in Neely et al [20] and our hypothesis necessitating there to be relatively high clay content associated with recharge pathways. Below are considerations for the extension this work to other aquifers systems and InSAR signatures.

We made the reasonable assumption, as demonstrated by a number of studies with focus in the San Joaquin Valley [4, 36–38], that the vertical displacement measurements are predominately attributed to poroelastic deformation. This may not be the case for other aquifer systems. For example, elastic loading due to changes in the mass of ice and snow in the Sierra Nevada have been shown to be the largest process driving seasonal changes in surface displacement in the Sacramento Valley to the north [36]. Thus, InSAR signatures found in other aquifer systems, where poroelastic deformation is not the dominant driver of seasonal surface displacements, may not reflect recharge or may need to be interpreted differently.

This study followed Neely et al [20] by fitting the time series with a sinusoidal function with an annual period. While this is simple and easily implemented, it cannot provide a sense of the dominant oscillation for deformation. For example, a time series where only seasonal groundwater extraction occurred would still result in a model that suggests both seasonal subsidence and uplift of equal magnitudes. Thus, this modeling approach describes the relative seasonal displacements. However, increases in groundwater elevations for both study sites (up to 20 m) during the wet water years [25] as well as true uplift observed in the InSAR time series (shown in figures S2–6) provided confidence that recharge occurred rather than a halting in groundwater storage declines. The use of more sophisticated methods such as polynomial fitting and principal component decomposition may improve the estimation for the timing to peak amplitude and be useful for investigating non-seasonal recharge events. Further, we focused on addressing the question of the necessary hydrogeological conditions in order to observe an InSAR signature in the San Joaquin Valley. Future detailed analysis of these InSAR signatures with the other components of the InSAR time series, the annual rate and peak amplitude, may provide deeper insight into the driving forces of deformation and the underlying heterogeneity of the aquifer system.

This investigation was limited by the availability of dense AEM surveys. Extension to other aquifer systems would likely face similar constraints. Serendipitously, our two sites demonstrated ideal endmembers of clay content and confining conditions with the Visalia site coinciding with an InSAR signature. The work accomplished here provides the foundational framework for the interpretation of InSAR signatures attributed to recharge. Collection of additional AEM data over other regions and aquifer systems, with and without InSAR signatures, may prove invaluable for further understanding the hydrogeological controls on the surface expression of river-source recharge and recharge in general.