Advancing the Limits of InSAR to Detect Crustal Displacement from Low-Magnitude Earthquakes through Deep Learning

: Detecting surface deformation associated with low-magnitude ( M w ≤ 5) seismicity using interferometric synthetic aperture radar (InSAR) is challenging due to the subtlety of the signal and the often challenging imaging environments. However, low-magnitude earthquakes are potential precursors to larger seismic events, and thus characterizing the crustal displacement associated with them is crucial for regional seismic hazard assessment. We combine InSAR time-series techniques with a Deep Learning (DL) autoencoder denoiser to detect the magnitude and extent of crustal deformation from the M w = 3.4 Gallina, New Mexico earthquake that occurred on 30 July 2020. Although InSAR alone cannot detect event-related deformation from such a low-magnitude seismic event, application of the DL method reveals maximum displacements as small as ( ± 2.5 mm) in the vicinity of both the fault and earthquake epicenter without prior knowledge of the fault system. This finding improves small-scale displacement discernment with InSAR by an order of magnitude relative to previous studies. We additionally estimate best-fitting fault parameters associated with the observed deformation. The application of the DL technique unlocks the potential for low-magnitude earthquake studies, providing new insights into local fault geometries and potential risks from higher-magnitude earthquakes. This technique also permits low-magnitude event monitoring in areas where seismic networks are sparse, allowing for the possibility of global fault deformation monitoring.


Introduction
Interferometric synthetic aperture radar (InSAR) is a geodetic technique that uses phase interferometry to produce two-dimensional fields of surface displacement along the line-of-sight (LOS) of the radar.These displacement maps, also referred to as interferograms, typically have footprints on the order of tens to hundreds of kilometers, and the technique itself can resolve deformation at fine resolution (e.g., centimeter scale).These characteristics make InSAR an invaluable geodetic tool for observing geophysical phenomena, including seismic hazards.Since its initial application to measuring crustal deformation from the Landers, California earthquake in 1992 [1], this technique has been applied to well over 100 earthquake studies [2] and this number is currently growing by 20-30 events per year [3].The actual signal being measured by InSAR, the apparent change in the distance between ground and satellite difference from two synthetic aperture radar (SAR) images, represents a combination of crustal deformation plus noise potentially stemming from changes in atmospheric conditions, topography, soil moisture, or vegetation [4][5][6].Thus, care must be taken to properly partition out the signal of interest from other non-seismic factors.In particular, the detection of small and transient seismic deformation has remained a challenge, often requiring interpretation of noisy InSAR time series using a priori knowledge of the possible sources of deformation (e.g., [7][8][9]).
Time-series analysis is a common approach for mitigating atmospheric and topographic noise in InSAR datasets as well as effects from soil moisture and vegetation conditions.Such analysis looks for pixel-wise trends in displacement over the time span of an InSAR dataset, smoothing out transient signals and improving the signal-to-noise ratio of the data.Several noise mitigation techniques combined with time-series analysis have been used to extract deformation signals from low-magnitude (M w < 6) events.Such approaches include reducing atmospheric interference through time-series analysis [10] or accounting for topographic-related tropospheric delays, which was achieved while analyzing a M w 5.9 earthquake in Tibet [11].Auxiliary information including Global Positioning System (GPS) data and data from multiple SAR satellites, as well as advanced coregistration techniques, were used to investigate displacement from a M w 3.9 earthquake in southern Italy [12].Finally, in at least one case, InSAR was used to investigate four low-magnitude earthquake events (M w ≤ 5.1) in Turkey where the use of advanced techniques was unnecessary, and therefore absent, due to the events occurring over ideal conditions for InSAR processing [13].For example, the extent of deformation in this region was small with good coherence due to minimal tropospheric delay.Additionally, most of the study area was flat and the spatial extent of the fringe was large with a regular pattern of deformation.
In the aforementioned studies, InSAR processing techniques were applied to estimate crustal displacement at the centimeter scale.To increase signal detection still further to subcentimeter scales, other techniques to improve the signal-to-noise ratio of the data must be used.Persistent Scatter (PS) InSAR is one such technique that has been used due to its ability to minimize the effects of temporal and geometric decorrelation, as well as the influence of the atmospheric phase screen (APS) [14].PS works by considering consistent scatterers, which are generated from the satellite radar signal bouncing off stationary objects [14][15][16].However, PS performs optimally only when large data stacks are available, high PS density is present, and crustal deformation is minimal.The use of PS is thus problematic in highly vegetated areas or steep terrain that leads to low scattering densities, which often prevent the successful application of this technique.Moreover, PS analysis is computationally very expensive [16] and relies on knowledge of deformation smoothness in time and/or space, making it most suitable for monitoring stable deformation.Accounting for distributed scatterers (DS) (e.g., desert areas, crop fields, etc.) is a related technique that captures typically weaker signals than PS, but still may improve signal detection in InSAR analyses when smaller levels of displacement are expected.By selecting pairs using a minimum spanning tree framework weighted by a measure of image quality, one can obtain similar results using methods that consider both PS and DS (e.g., SqueeSAR; [17]) [18].However, this method has not been tested using data from ScanSAR/TOPSAR imaging modes (e.g., Sentinel-1), nor has it been applied to fault activity, which may have a more complicated signal than geothermal signatures.
More recent work has shown promise for the use of a deep learning (DL)-based method, which was used to automatically detect and extract ground deformation signals from InSAR time series [18].This technique relies on a fully convolutional deep auto-encoder to recognize and extract a signal from noisy InSAR time series.Using this technique, deformation on faults can still in principle be recovered from an InSAR signal for signal-tonoise ratios (SNRs) as low as 1%.Therefore, when combined with standard InSAR output, this method can potentially be used to detect millimeter-scale displacement.
In this study, we examine the potential of the aforementioned DL method for detecting sub-centimeter crustal displacement from the 30 July 2020 M w 3.4 earthquake in Gallina, NM in the vicinity of the Nacimiento -Gallina fault arch.This earthquake occurred less than 70 km from several northern New Mexico communities, as well as the Valles Caldera National Preserve and Bandelier National Monument, which together hold geological and cultural significance.A better understanding of seismicity along this fault system is critical, as a recent study provided evidence that the Gallina fault is currently active [19] and supported earlier, poorly constrained estimates that this fault may be capable of accommodating an M w 6-7 event in the future [20].

Site Background
The M w 3.4 Gallina, NM earthquake occurred on 30 July 2020 near the southern end of the Gallina fault and the northern end of the Nacimiento fault at 36 • 13.09 ′ N, 106 • 50.87 ′ W at a depth of 11.4 km (location determined from Los Alamos Seismic Network (LASN) data, as shown in Figure 1 with faults from the United States Geological Survey (USGS) Geologic Map Database) [21,22].The Gallina and Nacimiento faults are largely Laramide structures (80-50 Ma) [23,24], although there is evidence for some post-Laramide, Miocene slips [25].
It is not yet resolved whether the Gallina and Nacimiento faults are currently active (i.e., have Quaternary offset).In two probabilistic seismic hazard analyses for the Los Alamos National Laboratory [20,26], and a study on the nearby El Vado Dam [27], both the Gallina and Nacimiento faults were conservatively interpreted to be active, despite no clear evidence of Quaternary offset, largely due to microseismicity observed within and east of these faults [28,29].More recently, the first clear geologic evidence for Quaternary deposit offset was provided across the northern part of the Gallina fault and also extended the fault 14 km further north [19].No similar studies have been performed near the southern end of the Gallina fault or the northern end of the Nacimiento fault, near the location of the 30 July 2020 earthquake.
The Nacimiento Fault has been interpreted as a predominantly east-dipping thrust and reverse fault, transitioning northward along strike through vertical to a west-dipping normal fault near its northern end [25].Its dip is thought to be steep at depth, flattening upwards, with the shallow section accommodating a westward slip over the San Juan Basin sediments [30].The Gallina fault is mapped as being nearly vertical and has both predominant down-to-the-west and minor down-to-the-east offsets [25].Detection of ground deformation from the July 2020 event, if it aligns with the known faults, would bolster both the evidence for the Gallina and/or Nacimiento faults being active and provide a means to constrain fault geometry.Caldera National Preserve and Bandelier National Monument, which together hold geological and cultural significance.A better understanding of seismicity along this fault system is critical, as a recent study provided evidence that the Gallina fault is currently active [19] and supported earlier, poorly constrained estimates that this fault may be capable of accommodating an Mw 6-7 event in the future [20].

Site Background
The Mw 3.4 Gallina, NM earthquake occurred on 30 July 2020 near the southern end of the Gallina fault and the northern end of the Nacimiento fault at 36°13.09′N, 106°50.87′Wat a depth of 11.4 km (location determined from Los Alamos Seismic Network (LASN) data, as shown in Figure 1 with faults from the United States Geological Survey (USGS) Geologic Map Database) [21,22].The Gallina and Nacimiento faults are largely Laramide structures (80-50 Ma) [23,24], although there is evidence for some post-Laramide, Miocene slips [25].
It is not yet resolved whether the Gallina and Nacimiento faults are currently active (i.e., have Quaternary offset).In two probabilistic seismic hazard analyses for the Los Alamos National Laboratory [20,26], and a study on the nearby El Vado Dam [27], both the Gallina and Nacimiento faults were conservatively interpreted to be active, despite no clear evidence of Quaternary offset, largely due to microseismicity observed within and east of these faults [28,29].More recently, the first clear geologic evidence for Quaternary deposit offset was provided across the northern part of the Gallina fault and also extended the fault 14 km further north [19].No similar studies have been performed near the southern end of the Gallina fault or the northern end of the Nacimiento fault, near the location of the 30 July 2020 earthquake.
The Nacimiento Fault has been interpreted as a predominantly east-dipping thrust and reverse fault, transitioning northward along strike through vertical to a west-dipping normal fault near its northern end [25].Its dip is thought to be steep at depth, flattening upwards, with the shallow section accommodating a westward slip over the San Juan Basin sediments [30].The Gallina fault is mapped as being nearly vertical and has both predominant down-to-the-west and minor down-to-the-east offsets [25].Detection of ground deformation from the July 2020 event, if it aligns with the known faults, would bolster both the evidence for the Gallina and/or Nacimiento faults being active and provide a means to constrain fault geometry.[21].White lines are faults from the USGS Geologic Map Database [22].Map inset in upper right corner depicts the location of the study area in blue relative to the state of New Mexico.Map was created using Generic Mapping Tools v6 [31].

Data
We utilize data from the European Space Agency's C-band Sentinel-1 constellation due to its repeat global coverage and open-source status.Sentinel-1 has an average cadence of 12 days.We work with single-look complex (SLC) data from the Interferometric Wide (IW) mode, which is collected at roughly 5 by 20 m resolution over three subswaths, forming a full image with an overall swath size of 250 km [32].For our analysis, we use the third subswath of ascending track T49 and the second subswath of descending track T56.The corresponding unit look vectors, considered from a pixel on the ground to the satellite in the global Cartesian coordinate system, are ŝT49 = [0.70, 0.12, −0.70] and ŝT56 = [−0.63, 0.11, −0.77].We focus on data acquired two to three months pre-and post-earthquake event (March to October 2020), excluding months outside of this timeframe to avoid processing scenes affected by snow cover.
Pairs are chosen to maximize the coherence of the dataset by setting thresholds for both the orbital baseline (120 m) and temporal baseline (90 days).In total, we processed 45 interferometric pairs from 16 out of 19 different scenes (i.e., epochs) for ascending track T49 (3 scenes were deemed unfit for interferometric analysis due to poor correlation) and 57 pairs from 16 scenes spanning descending track T56 (Figures A1 and A2 in Appendix A) using the open-source Generic Mapping Tool Synthetic Aperture Radar (GMTSAR) InSAR processing software [33].We apply a 30 m digital elevation model (DEM) provided by the Shuttle Radar Topographic Mission (SRTM) to remove topographic contributions from each interferogram [34].Pairs are processed in batch using built-in capabilities from GMTSAR with ionospheric phase correction included via the split-spectrum method (e.g., [35]).The pairs are unwrapped using GMTSAR's implementation of Statistical-cost, Network-flow Algorithm for Phase Unwrapping (SNAPHU) [36].This results in pairwise measurements of unwrapped range change or LOS deformation referenced to the radar.In this convention, positive values of range change indicate subsidence while negative values indicate uplift.We note that a common convention for time-series analysis is to reference deformation to the ground, such that negative values of LOS deformation indicate subsidence.For clarity, we hereby refer to ground-referenced LOS deformation as LOS displacement (following time-series analysis convention) and radar-referenced LOS deformation as range change.
A comparison to global positioning system (GPS) data from the Geodetic Facility for the Advancement of Geoscience (GAGE) network suggests our InSAR data have an accuracy on the order of millimeters (Appendix B).

Methods
The overall workflow for our methods is shown in Figure 2. We start by processing our InSAR data into interferometric pairs, the results of which we refer to as observed range changes.We then apply time-series analysis to these pairs to improve the signal-to-noise ratio of our dataset.This process results in cumulative displacements, which we will refer to as TS displacements.We feed these TS displacements to our DL method, which enhances signals of interest.The resulting displacement maps from this method are referred to as DL displacements.Lastly, we apply spatial deformation modeling to these DL displacements using an analytical displacement model for seismic applications.The resulting deformation fields from this model are referred to as modeled slips, which we further temporally analyze in terms of seismic slip.

2D Time-Series Analysis
After pair formation, the 2D average velocity field and cumulative deformation fields associated with each epoch in the dataset are estimated using pixel-wise time-series analysis.The open-source Python package MintPy [37] is used to perform a weighted least squares inversion informed by coherence (e.g., [38]), with additional tropospheric correction via atmospheric models supplied by PyAPS [39].Also included are DEM error estimation and orbital phase ramp estimation.

Deep Learning
The TS displacements created using the SBAS time-series analysis procedure described above functioned as input for the DL technique.The DL method applied in this study detects and extracts transient ground deformation signals from the time-series displacements using a fully convolutional autoencoder [18].More specifically, it has been trained to isolate spatial and temporal patterns of deformation resembling those found along faults following a seismic event from signals due to other modes of surface displacement, such as aquifer depletion or magma intrusion.At the same time, it removes noise associated with variability in atmospheric conditions.The autoencoder uses a 9-timestep window to parse the input TS displacement fields, outputting the cumulative deformation occurring during the time window (excluding the first and last steps, which are used as references by the model).The autoencoder has previously been trained on synthetic data, in which signals of low-magnitude displacement in fault zones were overlain by fields carefully designed to mimic atmospheric noise.The trained tool has shown remarkable success in earthquake detection.Within its training and validation datasets, it recovered displacement from interferograms with SNRs as low as 1%, and, in a case study along the North Anatolian fault in Turkey, it was successfully employed to detect millimeter-scale displacement from a slow slip event [18].We use our TS displacements and DEM as inputs to this DL model.This results in displacement fields containing refined signals of interest, which we refer to as DL displacements.

Slip Assessment via Spatial Deformation Modeling and Temporal Adjustment
In addition to identifying surface deformation signals amidst background noise in a seismically active area, we aim to identify whether these signals can be reasonably interpreted as seismically related.To do this, we apply spatial deformation modeling by fitting a seismic model to identify signatures of interest found in our DL displacements spanning the event window.We use an open-source package called the General Inversion of Phase Technique (GIPhT) [40].When supplied with a deformation model and accompanying initial estimates of parameters for such a model, GIPhT uses simulated annealing to solve the nonlinear inversion problem of estimating the model parameters that best describe the deformation observed in the InSAR data.Since we are particularly interested in whether observed signals may be explained seismically, we employ the Okada model [41].This model calculates the surface deformation generated by the slip of a solid over a rectangular patch at a given depth below the free surface of a homogeneous elastic half-space.We assume a standard Poisson's ratio of ν = 0.25 following similar studies in the area [20,42].Parameters of interest that are optimized in this formulation include the dimensions and spatial orientation of the rectangular patch, its depth, and the direction and magnitude of the slip.We first fit the average field of the DL displacements for both descending and ascending tracks to find the best-fitting structural parameters for the slip event.We then fit each individual DL displacement field with the optimized structural parameters constrained to estimate the amount of slip over each temporal pairing.
To analyze the relationship between estimated slip and seismic activity, we apply a time-series analysis technique called temporal adjustment.This technique converts the pairwise estimates of slip to cumulative slip over each epoch of the dataset [43].Assuming separable temporal and spatial dependence, the displacement field u is then defined by the product of a distinct temporal function f (t) and spatial function G(x): In our case, G(x) is represented by a single Okada or collection of Okada source models.We examine several different temporal representations f (t) to best describe temporal trends in the slip.First, we consider a constant rate parameterization with rate a for an epoch t i given starting epoch t 0 : We also consider piecewise-linear parameterizations: We examine a two-segment (m = 2) parameterization f 2 (t i ) with a break at the event on 30 July.We also examine a three-segment (m = 3) parameterization f 3 (t i ) with breaks on 1 July and 30 July.
While we also consider nonlinear representations such as exponential decay and polynomial parameterizations, initial analysis indicates these functions are not well suited to characterizing the temporal nature of slip in this study and are thus not discussed further.

InSAR Time Series
The TS cumulative displacement fields derived from time-series analysis on descending track T56 are shown in Figure 3, with the corresponding ascending track T49 analysis having similar results.The sequence of TS displacement fields shows cumulative deformation from the initial data point (t 0 ), where the date of the final timestamp for each deformation field is shown on the bottom right of each image.There is no apparent deformation surrounding the epicenter of the earthquake (denoted by the outlined star) in either the ascending or descending track.Furthermore, both results are seemingly dominated by an atmospheric gradient effect present from southwest to northeast in TS displacement fields.Consequently, the TS displacement fields derived from our InSAR time-series analysis do not on their own reveal any meaningful deformation present at or in the vicinity of the earthquake epicenter.

Signal Refinement via Deep Learning
To further extract any subtle signals occurring in the data, we apply the DL autoencoder to our TS displacements. Figure 4 shows the average DL displacement rate fields derived for both tracks along with the LASN-located epicenter and fault traces obtained from the USGS Geologic Map Database [22].As opposed to the fields in Figure 3, we see noticeable displacement in several sections of the field.We focus our attention on deformation signatures that fall along fault lines and that have paired uplift-subsidence patterns following typical seismic-related deformation.This leads to a site along the Nacimiento Fault system (Figure 4), which we denote as the Nacimiento Fault (NF) site.

Signal Refinement via Deep Learning
To further extract any subtle signals occurring in the data, we apply the DL autoencoder to our TS displacements. Figure 4 shows the average DL displacement rate fields derived for both tracks along with the LASN-located epicenter and fault traces obtained from the USGS Geologic Map Database [22].As opposed to the fields in Figure 3, we see noticeable displacement in several sections of the field.We focus our attention on deformation signatures that fall along fault lines and that have paired uplift-subsidence patterns following typical seismic-related deformation.This leads to a site along the Nacimiento Fault system (Figure 4), which we denote as the Nacimiento Fault (NF) site.When further analyzing the DL results at the Nacimiento Fault site, we find that the model detects surface deformation between the Nacimiento and Gallina faults in both the ascending and descending tracks (Figure A3 in Appendix C).Cumulative displacement measurements within 9-step windows occurring prior to the July 30th event show paired subsidence and uplift signatures that mirror the configuration of the Gallina and Nacimiento faults at this site.Displacement reaches a maximum uplift of ~1.5 mm and maximum subsidence of ~2.5 mm during this time.While cumulative displacements from 9-step windows that encompass the July 30th event maintain the characteristic shape seen in pre-event displacements, they consistently diminish in magnitude as time progresses until cross-sections reveal almost no discernible displacement in the final estimates.
We note that the model has no a priori knowledge of the faults, so its detection of displacement at the Nacimiento Fault site that mirrors the known configuration of the fault system increases our confidence in attributing the signal to the earthquake.When further analyzing the DL results at the Nacimiento Fault site, we find that the model detects surface deformation between the Nacimiento and Gallina faults in both the ascending and descending tracks (Figure A3 in Appendix C).Cumulative displacement measurements within 9-step windows occurring prior to the July 30th event show paired subsidence and uplift signatures that mirror the configuration of the Gallina and Nacimiento faults at this site.Displacement reaches a maximum uplift of ~1.5 mm and maximum subsidence of ~2.5 mm during this time.While cumulative displacements from 9step windows that encompass the July 30th event maintain the characteristic shape seen in pre-event displacements, they consistently diminish in magnitude as time progresses until cross-sections reveal almost no discernible displacement in the final estimates.
We note that the model has no a priori knowledge of the faults, so its detection of displacement at the Nacimiento Fault site that mirrors the known configuration of the fault system increases our confidence in attributing the signal to the earthquake.

Spatial Deformation Modeling and Slip Analysis
We further investigate potential causes of the DL displacement signals identified at the Nacimiento Fault site by considering whether the observed displacement could have been caused by co-seismic slip occurring at the time of the event.We investigate this using the Okada model within GIPhT, as described in Section 2.3.3.We start by finding the average DL displacement rate field for each track.We then use simulated annealing via

Spatial Deformation Modeling and Slip Analysis
We further investigate potential causes of the DL displacement signals identified at the Nacimiento Fault site by considering whether the observed displacement could have been caused by co-seismic slip occurring at the time of the event.We investigate this using the Okada model within GIPhT, as described in Section 2.3.3.We start by finding the average DL displacement rate field for each track.We then use simulated annealing via GIPhT to find Okada model parameters that best explain the deformation measured in both tracks' average displacement rate fields.Once these parameters are identified, we then apply the model to each individual DL displacement field, holding constant the structural parameters of the model (e.g., location and geometry) and inverting only for slip.The similarity between ascending and descending DL displacement fields allowed for one model to describe the structure of the deformation signatures in both tracks.
Based on the deformation pattern at the NF site, we fit two Okada models.The best-fitting model parameters are shown in Table 1.We define a dimensionless misfit of the model to the DL displacement data using the square root of the reduced χ 2 -test statistic [44].We find the misfit to be χ = 1.6 for both ascending and descending average DL displacement rate fields.Figure 5 shows the corresponding model compared to the descending average DL displacement field.then apply the model to each individual DL displacement field, holding constant the structural parameters of the model (e.g., location and geometry) and inverting only for slip.The similarity between ascending and descending DL displacement fields allowed for one model to describe the structure of the deformation signatures in both tracks.
Based on the deformation pattern at the NF site, we fit two Okada models.The bestfitting model parameters are shown in Table 1.We define a dimensionless misfit of the model to the DL displacement data using the square root of the reduced  -test statistic [44].We find the misfit to be  1.6 for both ascending and descending average DL displacement rate fields.Figure 5 shows the corresponding model compared to the descending average DL displacement field.

Table 1.
Best-fitting Okada model parameters for the NF Site.Following standard Okada conventions [40,41], dip is measured downward relative to the top edge of the fault plane, strike is measured clockwise from North, and rake is measured in-plane and counterclockwise from the strike.

Parameter
Estimate Uncertainty   Once we have estimates of slip for each pairwise DL displacement field at each site, we then perform temporal adjustment to model the resulting epoch-wise slip and determine temporal trends in the slip over time according to our f 1 (t), f 2 (t), and f 3 (t) parameterizations.We model each Okada source separately per site.Our results for both sites are described in Table 2.We find that the piecewise-linear parameterization using three breaks ( f 3 (t)) provides the best fit for both Okada sources (Figures 6 and 7), although the first Okada source (O1) is also reasonably described using the two-segment parameterization f 2 (t).To determine whether the added complexity of including a break at the start of July is meaningful, we perform an F-test for model complexity at 95 percent confidence [44].We find an F-test score of F calc = 7.8 as compared to a critical value of F α=0.05 = 7.8 with degrees of freedom do f 1 = 8 and do f 2 = 7.We thus reject the null hypothesis that f 2 (t) and f 3 (t) provide similar measures of fit, indicating that f 3 (t) is significantly better at explaining the temporal slip pattern for Okada 1.
Once we have estimates of slip for each pairwise DL displacement field at each site, we then perform temporal adjustment to model the resulting epoch-wise slip and determine temporal trends in the slip over time according to our   ,   , and   parameterizations.We model each Okada source separately per site.Our results for both sites are described in Table 2.
Table 2. Results from temporal adjustment on slip estimated at NF site from the best-fitting Okada 1 (O1) and Okada 2 (O2) sources.O1 is associated with slip at a centroid depth of ~207 m and O2 is associated with slip at a centroid depth of ~370 m.

Parameterization
Planar We find that the piecewise-linear parameterization using three breaks (  ) provides the best fit for both Okada sources (Figures 6 and 7), although the first Okada source (O1) is also reasonably described using the two-segment parameterization   .To determine whether the added complexity of including a break at the start of July is meaningful, we perform an F-test for model complexity at 95 percent confidence [44].We find an F-test score of  7.8 as compared to a critical value of  .7.8 with degrees of freedom  8 and  7 .We thus reject the null hypothesis that   and   provide similar measures of fit, indicating that   is significantly better at explaining the temporal slip pattern for Okada 1.
We also consider similarities in temporal trends between the Okada 1 and Okada 2 sources.While there is no significant difference in slip rate estimates for   based on uncertainty bounds, the best-fitting parameters for the   and   parameterizations are markedly different between sources, suggesting that there are more complicated dynamics at this site.

Discussion
The detection of displacement associated with the Gallina earthquake is the smallest magnitude earthquake and displacement identified yet using InSAR, with detected displacements as small as 2.5 mm in the area between the Nacimiento and Gallina faults (e.g., Nacimiento Fault site), where deformation signatures were consistently captured amongst co-seismic pairs.Although the DL technique applied does not definitively link the estimated displacement directly to the earthquake, our confidence in this interpretation is supported by the fact that the denoising autoencoder identified maximum deformation in an area between the sections of the Nacimiento and Gallina fault systems that are just West of the reported epicenter, despite having no prior knowledge of the fault location.Moreover, the timing of the displacement we observed correlates with the timing of the earthquake, with deformation peaking just prior to the seismic event, and then gradually relaxing throughout the co-seismic time intervals.
Taking further into consideration only those DL-identified displacements that occur along faults in both the ascending and descending tracks, we focus on the region between the southern end of the Gallina Fault and the Nacimiento (e.g., Figure 4).This area shows complex spatial deformation patterns that suggest two areas of localized slip, which our inverse spatial deformation modeling suggests are related to movement along two normal faults.While both of these slip sources are aligned with the SW-NE trend of the Nacimiento fault bend, the western slip source (Okada 1) has a steeper dip (~80°) as compared to its eastern counterpart (Okada 2, with a dip of ~50°).Additionally, both sources are estimated at shallow depths occurring between 200 and 400 m.Taking the strikes of both sources to be roughly equal, this suggests that the sources' corresponding faults intersect at a depth of ~2 km roughly 300 m NE of Okada 1's location.Temporal adjustment on Figure 7. Results of temporal adjustment on NF Okada 2 slip estimates occurring at ~370 m depth using f 3 (t) parameterization.Plotting conventions as in Figure 6.
We also consider similarities in temporal trends between the Okada 1 and Okada 2 sources.While there is no significant difference in slip rate estimates for f 1 (t) based on uncertainty bounds, the best-fitting parameters for the f 2 (t) and f 3 (t) parameterizations are markedly different between sources, suggesting that there are more complicated dynamics at this site.

Discussion
The detection of displacement associated with the Gallina earthquake is the smallest magnitude earthquake and displacement identified yet using InSAR, with detected displacements as small as 2.5 mm in the area between the Nacimiento and Gallina faults (e.g., Nacimiento Fault site), where deformation signatures were consistently captured amongst co-seismic pairs.Although the DL technique applied does not definitively link the estimated displacement directly to the earthquake, our confidence in this interpretation is supported by the fact that the denoising autoencoder identified maximum deformation in an area between the sections of the Nacimiento and Gallina fault systems that are just West of the reported epicenter, despite having no prior knowledge of the fault location.Moreover, the timing of the displacement we observed correlates with the timing of the earthquake, with deformation peaking just prior to the seismic event, and then gradually relaxing throughout the co-seismic time intervals.
Taking further into consideration only those DL-identified displacements that occur along faults in both the ascending and descending tracks, we focus on the region between the southern end of the Gallina Fault and the Nacimiento (e.g., Figure 4).This area shows complex spatial deformation patterns that suggest two areas of localized slip, which our inverse spatial deformation modeling suggests are related to movement along two normal faults.While both of these slip sources are aligned with the SW-NE trend of the Nacimiento fault bend, the western slip source (Okada 1) has a steeper dip (~80 • ) as compared to its eastern counterpart (Okada 2, with a dip of ~50 • ).Additionally, both sources are estimated at shallow depths occurring between 200 and 400 m.Taking the strikes of both sources to be roughly equal, this suggests that the sources' corresponding faults intersect at a depth of ~2 km roughly 300 m NE of Okada 1's location.Temporal adjustment on pairwise slips estimated during the span of our dataset suggests that both sources have temporal trends associated with the July 30th event; we see a change in slip for both slip sources starting a month prior to the event and immediately following the event.While both sources show an increase in slip at the start of our dataset and a marked decrease in slip post-event, their slip responses in the month preceding the event differ.For the western deformation (Okada 1), slip starts to decrease, while the eastern deformation (Okada 2) shows a drastic increase in slip.We note here that while we select July 1st as the start of our second time interval for the f 3 (t) parameterization, the lack of temporal data coverage in June limits us from distinguishing between any significance of a temporal break in June as compared to the start of July.
The timing of the deformation, which starts prior to the seismic event, makes it unlikely that the observed deformation is a direct result of the earthquake.However, the close correlation in both the time and space of the observed deformation and the earthquake suggests that both are related to slip on the same fault system.If so, then we observe subtle precursor slip, leading to the intriguing possibility that this method could potentially be used to identify regions of aseismic slip precursors to earthquakes.This would require additional validation at other sites and time periods.
While we focus our analysis on activity along known fault systems, we note that there are several signals in the NE portion of Figure 4 that occur in both the ascending and descending tracks (further distinguished in Appendix D, Figure A4).While these signals do not occur along a known fault system, their distinguishable presence in both tracks suggests that activity may be occurring here as well.The faulting systems in this area are not well constrained [45], and there is the possibility that our method is detecting movement along faults that have not yet been mapped.While outside the scope of this paper, examining the signatures and geology of this region would better inform on the potential for our method to also help distinguish undetected faults.
One possible approach for constraining the relationship between deformation signatures and geology in this region is to further decompose the InSAR-observed LOS displacements into their vertical and horizontal components.Considering that the primary goal of this study was to assess the potential for sub-centimeter, seismic-related surface deformation detection using the DL technique in ref. [18] and that this technique utilizes LOS deformation, displacement decomposition analysis was not explored in this study.However, such an analysis may be warranted in future studies of this area if there are ascending and descending scenes available with near-coincident temporal pairings.
A comparison of the region of displacement with the highest-resolution geologic map for the area [46] (Figure A6 in Appendix E) indicates that the westernmost boundary between uplift and subsidence in both the ascending and descending track deformation areas follow a line at the headscarp of a mapped landslide.The co-location of the displacement region within a landslide formation raises the possibility that the mechanism of displacement was not seismically related and could perhaps be associated with the instability of the landslide.In order to claim detection of deformation associated with the smallest magnitude earthquake, it is necessary to distinguish between landslide and fault-related deformation.
Given this ambiguity in geologic relations, we performed a field check of the area of observed InSAR deformation on 10 November 2022 and found that the mapped landslide cannot be modern, as the source area for the deposits has been eroded away.In addition, the area of uplift is a topographically high area relative to the western area of subsidence.Neither landsliding nor other geomorphic processes can explain the uplift of a topographically high area, while this is readily explainable by faulting.In addition, the geometry of deformation observed in this study, with 2 uplift-subsidence pairs, is more consistent with fault displacement.Finally, the feature is interpreted to be a landslide scarp on the map in ref. [46] is continuous with a feature interpreted to be the Gallina fault along the strike further east.Thus, the existing geologic data support the interpretation of this InSAR deformation presented here to be related to slip on the Gallina fault, and not landslide-related.
The amount of atmospheric attenuation of a SAR signal is inversely proportional to the wavelength of the satellite [47].This means that the application of the DL method described here would be most useful for InSAR processing in conjunction with smaller wavelength X-or C-band missions where this technique is likely to be more effective at filtering noise related to atmospheric artifacts from the InSAR image.The utility of the method here is most obvious for the ongoing Sentinel-1 mission, as well as the recent C-band Radarsat Constellation Mission and the X-band family of COSMO-SkyMed, KOMPSAT-5, TerraSAR-X, TanDEM-X, or PAZ [2], where atmospheric delays are often greater than potential signals of interest.Although these other wavelengths were not utilized in this study due to either a lack of availability of open-source data or a lack of data coincident with our study area, they may be useful for follow-on analysis or investigations of other low-magnitude earthquake events where the DL approach is warranted.In addition, there is potential to utilize data from the upcoming L-and S-band NASA-ISRO SAR mission, the joint Earthobserving satellite mission between NASA and the Indian Space Agency [48], which may suffer less from atmospheric delays due to the longer wavelength of the radar.Given the ability of the DL model to extract deformation as a function of the signal-to-noise ratio, its application to less noisy data even from longer wavelength SAR missions may yield even smaller deformation detection capability, thereby still demonstrating its usefulness.We do note that the removal of ionospheric contributions through other techniques (e.g., split spectrum method [35]) would need to be given special consideration for longer wavelength data sources, as the DL technique applied in this study does not explicitly address ionospheric effects.
While the standard InSAR analysis techniques have not found any deformation located at the earthquake site location, the DL technique shows promise for future applications to identify deformation related to low-magnitude earthquakes.This study is at the forefront of low-magnitude, small-displacement earthquake studies from InSAR.To demonstrate this, Figure 8 shows the measured LOS crustal displacement for a number of earthquakes and the associated earthquake magnitudes from previous InSAR processing of Sentinel-1 studies conducted from 2015 until recently.The vast majority of these studies have been on earthquakes with a magnitude greater than 6.The Gallina earthquake that was targeted for this study is not only the smallest magnitude earthquake represented (M w = 3.4), but the application of the DL technique has improved the LOS displacement detection by more than an order of magnitude relative to other studies.
Through further applications to similar earthquake settings, the technique applied here promises to improve the understanding of earthquake hazard potential, should a higher magnitude earthquake occur, as well as provide information on the local fault geometries through additional modeling.Moreover, the application of this technique to higher-magnitude events holds the potential to reveal more about the extent of crustal displacement by highlighting new signals of displacement that were previously undetected.Such a feat was already demonstrated by the application of this technique to the North Anatolian Fault where the model recovered 2 mm of previously undetected slip along a fault [18].These new findings of displacement extent uncovered by this method led to the determination that the length of the fault that accommodated slip was twice what was previously thought and new fault lineaments were identified.Further studies on the application of this method to other small earthquakes that have occurred in the Gallina-Nacimiento fault system could be used to investigate similar features.demonstrate this, Figure 8 shows the measured LOS crustal displacement for a number of earthquakes and the associated earthquake magnitudes from previous InSAR processing of Sentinel-1 studies conducted from 2015 until recently.The vast majority of these studies have been on earthquakes with a magnitude greater than 6.The Gallina earthquake that was targeted for this study is not only the smallest magnitude earthquake represented (Mw = 3.4), but the application of the DL technique has improved the LOS displacement detection by more than an order of magnitude relative to other studies.InSAR pairs.The red circle represents the Gallina earthquake targeted for this study, which is an order of magnitude smaller in LOS displacement than typical studies.Associated references used to obtain this information are provided in Appendix F (Table A1).

Conclusions
We have demonstrated the effectiveness of combining InSAR time-series analysis with a DL autoencoder denoiser to detect surface deformation associated with low-magnitude (M w ≤ 5) seismicity.Applying this technique to Sentinel-1 InSAR pairs spanning the M w = 3.4 Gallina, New Mexico earthquake that occurred on 30 July 2020, we can detect surface deformation signals as subtle as 2.5 mm in the area between the Nacimiento and Gallina faults that are otherwise undetectable with standard time-series analysis techniques.Additionally, the spatial and temporal patterns of the displacement coincide with the seismic event, suggesting a relationship between the observed subsidence and the July 30th event itself.Additionally, the DL autoencoder was able to identify several locations of displacement in both the ascending and descending tracks that were not aligned with known faults, suggesting that this could be a method for fault identification as well.Additional analysis across other low-magnitude events will further determine the extent of this technique for low-magnitude seismic risk characterization and fault identification.
accessed on 8 April 2024), as well as through the Sentinel Scientific Data Hub.SRTM DEM data for InSAR processing was accessed free of charge through GMTSAR's DEM retrieval interface (https://topex.ucsd.edu/gmtsar/demgen/,last accessed on 1 October 2021).The GAGE GPS time-series data used for displacement validation is freely available from the Nevada Geodetic Laboratory (http://geodesy.unr.edu/index.php,last accessed on 20 May 2024).The DEM used for the k-means cluster was obtained from the Earth Data Analysis Center of the University of New Mexico (https://edac.unm.edu/, last accessed on 20 May 2024).Fault information is made freely available by the United States Geological Survey (USGS) via the ScienceBase-Catalog (https: //www.sciencebase.gov/catalog/item/4f4e496ee4b07f02db5a354e,last accessed on 6 March 2024).

Appendix B. InSAR-Measured Displacement Accuracy Assessment Using GPS
We use GPS data to perform a calibration check on our InSAR measurements.There are two GAGE network continuous GPS stations, RG09 (36.305°N, 107.056°W) and RG10 (36.451°N, 106.539°W), within the coverage of our ascending T49 and descending T56 tracks [49,50].Both stations have been analyzed according to standard procedures, with their resulting time series of relative positions made freely available [51].However, both stations only provide temporal coverage until mid-2019, which is outside the range of our 2020 analysis.To perform our validation, we instead use InSAR pairs spanning the same months used in our 2020 analysis but from T49 and T56 scenes acquired in 2018.This resulted in 46 pairs derived from 15 ascending track T49 scenes spanning from 4 March 2018 to 24 September 2018 and 36 pairs from 12 descending track T56 scenes spanning from 10 March 2018 to 1 August 2018.
We perform our analysis similar to previous studies [52][53][54].In order to account for the relative nature of InSAR measurements, we assign GPS station RG09 as a reference station and treat RG10 as being within a deforming region.To arrive at estimates of InSAR-observed range change within the deforming region of RG10 with respect to RG09 for each InSAR pair, we compare InSAR-observed range change averaged over a 1 km 2 area centered around station RG10 to InSAR-observed range change averaged over a 1 km 2 area centered around station RG09.To compare our findings with the corresponding GPS measurements, we first convert GPS-observed displacements to range change using the unit pointing vectors described in Section 2.2.The GPS measurements (now in terms of range) are then compared at each station across time, corresponding with the epochs of each InSAR pair, to calculate each GPS station's measure of relative range change over the span of each InSAR pair.We are then able to find the difference of the GPS-measured range change estimates between the stations to derive GPS-measured range changes from RG10 with respect to RG09, giving a direct comparison to our InSAR-observed range changes at RG10 with respect to RG09.We find a mean difference of 5 3 mm in the range changes derived from GPS and InSAR measurements from both our T49 pairs and T56 pairs.Moreover, using uncertainties for our InSAR measurements calculated using standard deviation and

Appendix B. InSAR-Measured Displacement Accuracy Assessment Using GPS
We use GPS data to perform a calibration check on our InSAR measurements.There are two GAGE network continuous GPS stations, RG09 (36.305 • N, 107.056 • W) and RG10 (36.451 • N, 106.539 • W), within the coverage of our ascending T49 and descending T56 tracks [49,50].Both stations have been analyzed according to standard procedures, with their resulting time series of relative positions made freely available [51].However, both stations only provide temporal coverage until mid-2019, which is outside the range of our 2020 analysis.To perform our validation, we instead use InSAR pairs spanning the same months used in our 2020 analysis but from T49 and T56 scenes acquired in 2018.This resulted in 46 pairs derived from 15 ascending track T49 scenes spanning from 4 March 2018 to 24 September 2018 and 36 pairs from 12 descending track T56 scenes spanning from 10 March 2018 to 1 August 2018.
We perform our analysis similar to previous studies [52][53][54].In order to account for the relative nature of InSAR measurements, we assign GPS station RG09 as a reference station and treat RG10 as being within a deforming region.To arrive at estimates of InSARobserved range change within the deforming region of RG10 with respect to RG09 for each InSAR pair, we compare InSAR-observed range change averaged over a 1 km 2 area centered around station RG10 to InSAR-observed range change averaged over a 1 km 2 area centered around station RG09.To compare our findings with the corresponding GPS measurements, we first convert GPS-observed displacements to range change using the unit pointing vectors described in Section 2.2.The GPS measurements (now in terms of range) are then compared at each station across time, corresponding with the epochs of each InSAR pair, to calculate each GPS station's measure of relative range change over the span of each InSAR pair.We are then able to find the difference of the GPS-measured range change estimates between the stations to derive GPS-measured range changes from RG10 with respect to RG09, giving a direct comparison to our InSAR-observed range changes at RG10 with respect to RG09.We find a mean difference of (5 ± 3) mm in the range changes derived from GPS and InSAR measurements from both our T49 pairs and T56 pairs.Moreover, using uncertainties for our InSAR measurements calculated using standard deviation and uncertainties for our GPS measurements calculated using error propagation, we find overlap between µ ± 1σ ranges between all GPS-derived and InSAR-derived range change measurements per pair, suggesting there is no significant difference in the measurements obtained from each modality.This analysis was repeated using 250 m 2 and 500 m 2 areas for InSAR-observed range change sampling, with no difference in the results.

Appendix C. DL Displacement Results at Nacimiento Fault Site
Remote Sens. 2024, 16, x FOR PEER REVIEW 17 of 25 uncertainties for our GPS measurements calculated using error propagation, we find overlap between  1 ranges between all GPS-derived and InSAR-derived range change measurements per pair, suggesting there is no significant difference in the measurements obtained from each modality.This analysis was repeated using 250 m 2 and 500 m 2 areas for InSAR-observed range change sampling, with no difference in the results.

Appendix C. DL Displacement Results at Nacimiento Fault Site
Figure A3.Pairwise LOS displacement results from DL analysis on ascending track T49 and descending track T56 for the Nacimiento Fault site (e.g., region denoted by the square in Figure 4).[46].Also shown are fault locations from the USGS Quaternary Fault Database [56] and the New Mexico USGS Fault Database [22]

igure 1 .
Satellite image (Google Earth Pro v7.3, image taken on 19 September 2017 courtesy of Landsat/Copernicus) showing the epicenter of the 30 July 2020 event from LASN, denoted with a yellow circle

Figure 2 .
Figure 2. Flowchart showing the major steps for our displacement analysis and the names of LOS displacement resulting from each step.2.3.1.2D Time-Series Analysis After pair formation, the 2D average velocity field and cumulative deformation fields associated with each epoch in the dataset are estimated using pixel-wise time-series analysis.The open-source Python package MintPy [37] is used to perform a weighted least

Figure 2 .
Figure 2. Flowchart showing the major steps for our displacement analysis and the names of LOS displacement resulting from each step.

Figure 3 .
Figure 3. Example of TS displacement fields from time-series analysis on InSAR pairs from descending track T56.(a) LOS cumulative displacement from the start of the descending track T56 dataset (11 March 2020) to 2 August 2020, just after the event on 30 July (LASN-located epicenter denoted with a triangle).The TS deformation field is overlaid on topography from an SRTM DEM.(b) Full TS displacement fields for the descending track analysis.Geographic coverage and color scale as in (a).

Figure 3 .
Figure 3. Example of TS displacement fields from time-series analysis on InSAR pairs from descending track T56.(a) LOS cumulative displacement from the start of the descending track T56 dataset (11 March 2020) to 2 August 2020, just after the event on 30 July (LASN-located epicenter denoted with a triangle).The TS deformation field is overlaid on topography from an SRTM DEM.(b) Full TS displacement fields for the descending track analysis.Geographic coverage and color scale as in (a).

Figure 4 .
Figure 4. Average DL-derived LOS displacement rate fields for ascending track T49 and descending track T56.The Nacimiento Fault site is denoted with a black square, the LASN-located epicenter is denoted with a black star, and faults from the USGS are denoted with black lines.

Figure 4 .
Figure 4. Average DL-derived LOS displacement rate fields for ascending track T49 and descending track T56.The Nacimiento Fault site is denoted with a black square, the LASN-located epicenter is denoted with a black star, and faults from the USGS are denoted with black lines.

Figure 5 .
Figure 5. Best-fitting model for NF Site using two Okada sources, as fitted on DL-derived average displacement rate fields from ascending track T49 and descending track T56.Pictured is the observed average displacement rate field from descending track T56, the modeled displacement rate field derived from the best-fitting model parameters for both Okada sources (O1 and O2), and the residual field (observed minus modeled).

Figure 5 .
Figure 5. Best-fitting model for NF Site using two Okada sources, as fitted on DL-derived average displacement rate fields from ascending track T49 and descending track T56.Pictured is the observed average displacement rate field from descending track T56, the modeled displacement rate field derived from the best-fitting model parameters for both Okada sources (O1 and O2), and the residual field (observed minus modeled).

Figure 6 .
Figure 6.Results of temporal adjustment on NF Okada 1 slip estimates occurring at ~207 m depth using f 3 (t) parameterization.Modeled cumulative slip is shown with the black line, with corre- sponding 68% confidence intervals indicated by the dashed black lines.Temporal breaks in the parameterization are indicated with dashed, vertical green lines.Red segments represent pairwise measurements of slip with black dots distinguishing individual epochs per pair.Vertical blue bars indicate 1σ uncertainty in the pairwise slip measurements post-scaling using the square root of the variance.

Figure 6 .
Figure 6.Results of temporal adjustment on NF Okada 1 slip estimates occurring at ~207 m depth using   parameterization.Modeled cumulative slip is shown with the black line, with corresponding 68% confidence intervals indicated by the dashed black lines.Temporal breaks in the parameterization are indicated with dashed, vertical green lines.Red segments represent pairwise measurements of slip with black dots distinguishing individual epochs per pair.Vertical blue bars indicate 1 uncertainty in the pairwise slip measurements post-scaling using the square root of the variance.

Figure 7 .
Figure 7. Results of temporal adjustment on NF Okada 2 slip estimates occurring at ~370 m depth using   parameterization.Plotting conventions as in Figure 6.

Figure 8 .
Figure 8.Comparison of LOS crustal displacement and moment magnitudes studied using Sentinel-1InSAR pairs.The red circle represents the Gallina earthquake targeted for this study, which is an order of magnitude smaller in LOS displacement than typical studies.Associated references used to obtain this information are provided in Appendix F (TableA1).

Figure A3 . 25 Figure A4 .
Figure A3.Pairwise LOS displacement results from DL analysis on ascending track T49 and descending track T56 for the Nacimiento Fault site (e.g., region denoted by the square in Figure4).Each plot shows resulting cumulative displacement per 9-step window as estimated after applying our DL method to TS displacements within said time window.Star denotes the fault junction.

Figure A5 .
Figure A5.Elbow plot used to determine the optimum number of clusters [55].Appendix E. Geologic Relationship of Observed DisplacementFigureA6shows that the observed deformation is more consistent with faulting than landslides or other surficial deformation.

Figure A5 . 25 Figure A4 .
Figure A5.Elbow plot used to determine the optimum number of clusters [55].Appendix E. Geologic Relationship of Observed DisplacementFigureA6shows that the observed deformation is more consistent with faulting than landslides or other surficial deformation.

Figure A5 .
Figure A5.Elbow plot used to determine the optimum number of clusters [55].Appendix E. Geologic Relationship of Observed DisplacementFigureA6shows that the observed deformation is more consistent with faulting than landslides or other surficial deformation.

Figure A6 .
Figure A6.Overlay of deformation between 4 and 15 to 7 and 20 near the Nacimiento Fault site, with geology from the New Mexico Bureau of Mines & Mineral Resources[46].Also shown are fault locations from the USGS Quaternary Fault Database[56] and the New Mexico USGS Fault Database[22].

Table 1 .
[40,41]tting Okada model parameters for the NF Site.Following standard Okada conventions[40,41], dip is measured downward relative to the top edge of the fault plane, strike is measured clockwise from North, and rake is measured in-plane and counterclockwise from the strike.

Table 2 .
Results from temporal adjustment on slip estimated at NF site from the best-fitting Okada 1 (O1) and Okada 2 (O2) sources.O1 is associated with slip at a centroid depth of ~207 m and O2 is associated with slip at a centroid depth of ~370 m.