Hazard Implications of the 2016 Mw 5.0 Cushing, OK Earthquake from a Joint Analysis of Damage and InSAR Data

The Cushing Hub in Oklahoma, one of the largest oil storage facilities in the world, is federally designated as critical national infrastructure. In 2014, the formerly aseismic city of Cushing experienced a Mw 4.0 and 4.3 induced earthquake sequence due to wastewater injection. Since then, an M4+ earthquake sequence has occurred annually (October 2014, September 2015, November 2016). Thus far, damage to critical infrastructure has been minimal; however, a larger earthquake could pose significant risk to the Cushing Hub. In addition to inducing earthquakes, wastewater injection also threatens the Cushing Hub through gradual surface uplift. To characterize the impact of wastewater injection on critical infrastructure, we use Differential Interferometric Synthetic Aperture Radar (DInSAR), a satellite radar technique, to observe ground surface displacement in Cushing before and during the induced Mw 5.0 event. Here, we process interferograms of Single Look Complex (SLC) radar data from the European Space Agency (ESA) Sentinel-1A satellite. The preearthquake interferograms are used to create a time series of cumulative surface displacement, while the coseismic interferograms are used to invert for earthquake source characteristics. The time series of surface displacement reveals 4–5.5 cm of uplift across Cushing over 17 months. The coseismic interferogram inversion suggests that the 2016 Mw 5.0 earthquake is shallower than estimated from seismic inversions alone. This shallower source depth should be taken into account in future hazard assessments for regional infrastructure. In addition, monitoring of surface deformation near wastewater injection wells can be used to characterize the subsurface dynamics and implement measures to mitigate damage to critical installations.


Introduction
Cushing, Oklahoma, the "pipeline crossroads of the world," is home to the largest crude oil storage facility in the United States and the Keystone pipeline. Wastewater injection threatens this critical national energy infrastructure through both large, instantaneous ground acceleration from induced seismicity and through smaller, incremental, longer term movement from ground surface deformation [1].
Prior to 2009, seismicity in Oklahoma was concentrated along the Meers fault zone, southwest of the recently reactivated Nemaha and Wilzetta fault zones [2]. Since seismicity initiated in Cushing in 2014, three earthquake sequences have ruptured two fault zones ( Figure 1). The October 2014 sequence,

Synthetic Aperture Radar Data
To measure the surface displacement before and during the 2016 earthquake sequence, we process 32 European Satellite Agency (ESA) Copernicus Sentinel-1A (S1A) Synthetic Aperture Radar (SAR) Single Look Complex (SLC) images from the Alaska Satellite Facility (ASF). The SAR images spanning the Cushing area correspond to ascending relative orbit 34, frames 109-113 (north) and 114-118 (south), swaths 2 and 3 ( Figure 1). Table 1 lists the acquisition dates used for this study. The dates are separated by a dashed horizontal line, with preearthquake dates listed above and postearthquake dates listed below the horizontal line. We multilook the data, applying 12 looks in range and two looks in azimuth. For time series pairs, we apply a temporal baseline threshold of 12 days to minimize decorrelation due to vegetation (Table 2). We make an exception for a few pairs by increasing the temporal baseline up to 60 days to complete the network connectivity. For the source inversion, we pair images with a temporal baseline threshold of 36 days (Figure 2).

InSAR Processing
We perform interferometry on the SAR data using the Jet Propulsion Laboratory (JPL)/Caltech Interferometric SAR Scientific Computing Environment (ISCE) software [17]. The interferograms are corrected for topography using a 10-meter United States Geological Survey (USGS) National Elevation Dataset (NED) Digital Elevation Model (DEM). A Gaussian filter of strength 0.6 is applied to facilitate the unwrapping process. Interferograms are unwrapped using the Statistical-Cost, Network-Flow Algorithm for Phase Unwrapping (SNAPHU) algorithm [18][19][20] The time series of the cumulative surface displacement in the satellite line-of-sight (LOS) for the preseismic interferograms is constructed using the Generic InSAR Analysis Toolbox (GIAnT) software [21]. In GIAnT, we mask pixels whose phase coherence is less than 0.35, we correct for the atmospheric phase delay using the European Centre for Medium-Range Weather Forecasts (ECMWF) global atmospheric model, and we deramp the interferograms. After applying the corrections, we generate the time series using the New Small BAseline Subset (NSBAS) algorithm [22]. We select a reference region in the town of Carney, southwest of Cushing. This region was selected to ensure coherence over the time series.

Macroseismic Field
The EERI reconnaissance team specified a damage level from light to severe for structural, non-structural, and foundation damage for each photo based on post-earthquake site observations [7]. We assign a Modified Mercalli Intensity (MMI) to each set of photos and damage description representing a single building or site [23]. This process was based on matching qualitative descriptions of the MMI scale to the damage descriptions, similar to the technique used by the USGS to quantify damage from Did You Feel It? responses (Table 3) [24,25]. The MMI intensity that matched the observed damage then was assigned to that site. In the event that the observations and descriptions fell between two intensity levels, the higher intensity level was selected to avoid underestimating damage.

Surface Displacement Time Series
The surface displacement time series spans 5 June 2015 to 2 November 2016, ending five days before the Mw 5.0 earthquake. The cumulative displacement time series image reveals a broad region of 4 to 6 cm of uplift across the city of Cushing ( Figure 4) over 17 months. The broad deformation signal does not correlate with topography nor streams, rivers, or other waterbodies. These waterbodies, which are common sources of noise, do correlate with areas of decoherence ( Figure 5). Other broad regions of uplift are also observed in the cities of Stillwater and Perkins. We extract time series for pixels within two regions of interest, the Mw 5.0 earthquake ( Figure 6) and downtown Cushing ( Figure 7). We compare the displacement time series to seismicity using the Schoenball and Ellsworth [5] waveform-relocated earthquake catalog from May 2013 to November 2016, with a magnitude of completeness of 2.6. Note that the October 2014 Mw 4.3 earthquake featured in earlier figures was downgraded to a Mw 4.2 in the relocated earthquake catalog. Figures, henceforth, will reflect the downgraded magnitude. In addition to seismicity, we also compare the displacement data to monthly injection rates of of class II (oil and gas-related) Underground Injection Control (UIC) wells from the Oklahoma Corporation Commission (OCC) [27]. Wells are split into two types, type 2R for Enhanced Oil Recovery and type 2D for Salt-Water Disposal (SWD). Any well injecting above 300,000 bbl/month is classified as a high-rate injection well. A summary of the UIC wells within 6 km of our regions of interest are presented in Tables 4 and 5. Wells that have injected or are currently injecting into the Arbuckle have an additional column of information, injection period, where the table lists the time period that the well was injecting into the Arbuckle before the well was plugged and recompleted into a shallower formation, typically the Wilcox formation.  Although there is no immediate response of ground displacement to injection, there is an increase in displacement during the 2015 earthquake sequence. There was approximately 0.5 cm of uplift that occurred before the Mw 5.0 event. There are 15 UIC wells within 6 km of the Mw 5.0 earthquake ( Table 4). All but one well, 2R well API:3511923625, are 2D injection. 2D well API:3511923843, the only high-rate injection well, was high-rate in 2014, just two months before the first Cushing earthquake sequence of 2014.   For downtown Cushing, Figure 7 also shows steady uplift over time with a maximum displacement of 5.5 cm. The temporal signature of the displacement is strikingly similar to the signal for the Mw 5.0 earthquake region. The time series differ in the amplitude of displacement during the 2015 earthquakes and during the summer before the 2016 earthquake. During the 2015 earthquakes, the Mw 5.0 earthquake region, which is near the 2015 events, uplifts 2 cm more than downtown Cushing. Conversely, downtown Cushing experiences 2 cm of uplift more than the Mw 5.0 earthquake region a few months before the Mw 5.0 earthquake. There are nine low-rate UIC wells within 6 km of downtown Cushing, all but one is 2D ( Table 5). Six of the nine wells are within 6 km of the Mw 5.0 earthquake region. Unlike the area near the Mw 5.0 earthquake, none of the wells within 6 km of downtown Cushing are injecting into the Arbuckle formation. Thus, the column, "Injection Period," is omitted from Table 5.

Source Characterization
We construct six deramped coseismic interferograms ( Figure 8) ( Table 6). Interferograms 1 to 3 are the least saturated by an atmospheric signal. Interferogram 4, with a 5-day preseismic signal and a 7-day postseismic signal, best represents the earthquake deformation; however, it is heavily saturated by an atmospheric signal. We attribute the atmospheric signal to acquisition date 2 November 2016 due its presence in interferograms 4 through 6. Despite the atmospheric signal, we can see more than 1.6 cm of LOS displacement in opposite directions on both sides of the fault.  We average the six coseismic interferograms to increase the signal-to-noise ratio. Using the Pyrocko Kite software [28], we forward model LOS displacement by defining the event as a double-couple point source. We apply a genetic algorithm inversion scheme based on Covariance Matrix Adaptation (CMA-ES) for minimization between our observations and model in a least square sense (L2 norm). We perform 300 evaluations using the USGS National Earthquake Information Center (NEIC) parameters as a priori values. We list our best fit parameters and the the USGS NEIC parameters in Table 7. Our best fit parameters result in an underestimation of displacement (Figure 9), especially on the eastern side of the fault, suggesting the earthquake fault has a smaller dip, and thus has a greater oblique component than our almost purely strike-slip fault model. Forward modeling the USGS NEIC parameters results in an even greater underestimation of displacement on both sides of the fault, suggesting that the earthquake source is shallower than their estimated of 4.4 km ( Figure 10).

Damage Assessment
Earthquakes cause building and infrastructure damage through both dynamic shaking and permanent vertical or horizontal displacements of the ground surface. For large earthquakes, obvious manifestations of ground deformation effects include surface fault rupture, large lateral or vertical ground displacements associated with liquefaction, sand boils, landslides, and linear fissures, tilting or racking of homes, and foundation damage [29]. Indicators of more subtle impacts of ground deformation might include utility line breaks, damage to streets, curbs, sidewalks, plumbing, pipe bursts, and foundation cracks. Some examples of damage from ground acceleration are cracks in walls, spalling of exterior facade, permanent drift, and foundation damage. To assess damage from the 2016 Mw 5.0 Cushing earthquake, we remove duplicate data entries for the same damage site representing multiple photos of varying instances of damage. Thus, each point represents a unique damage site where a photo was taken ( Figure 11). The selected photos [7] of damage at American Legion, North Route 18 water main, Lion's Club, and Cimarron Tower from the 2016 Mw 5.0 earthquake, attest to the maximum MMI VII intensity ( Figure 12). For comparison with the modeled surface displacement, we select sites with damage indicators like pipe bursts, vertical offsets, and drops in elevation to indicate damage potentially caused or exacerbated by ground deformation (Table 8). Using the preferred parameters, we model the vertical and horizontal components of the surface displacement field ( Figure 13) and compare it to the ground failure damage sites. Sites 2 and 3, in the outskirts of town, occur in an area of either high vertical or horizontal displacement. The damage within downtown Cushing at sites 1,4, and 5 can potentially be attributed to a combination of differential settlement and ground acceleration. From Figure 13, we can see the majority of building damage from the earthquake is located in the Central Business District of downtown Cushing, OK. The buildings in this area were constructed between approximately 1910-1940 and are predominantly unreinforced masonry construction with wood framing. When the buildings were constructed, they were not designed to resist wind or seismic forces and thus the buildings may not have complete vertical load paths, or seismic detailing [7]. It is possible that some of the homes in more deformed regions have less damage because the buildings were constructed at a later date or because accelerations were locally lower there. Additionally, the varying geologic conditions between sites, including soil properties, have an impact on the intensity of dynamic shaking, which directly impacts the damage observed. There is a lack of strong motion instrumentation to quantify the differences in dynamic shaking intensity between sites [30].

Discussion
Since the late 2015 lift on the 40-year ban on U.S. oil exports, a surge in U.S. shale oil production has strained pipeline and storage capacity. Surface displacement through an induced earthquake or through wastewater injection can damage this critical infrastructure, significantly impacting the U.S. economy and the environment. Our time series analysis reveals 4 to 6 cm of positive LOS displacement (uplift) over ∼17 months (June 2015 to November 2016), a double in the uplift rate presented by Loesch and Sagan [13] of 4.4 cm over ∼49 months (December 2006 to January 2011). Loesch and Sagan [13] use an SBAS time series analysis of Japan Aerospace Exploration Agency (JAXA) Advanced Land Observing Satellite (ALOS) Phased Array L-band Synthetic Aperture Radar (PALSAR) images to show a broad uplift signal spanning 4 km. Our regions of uplift generally overlap; however, their uplift signal extends northeast by approximately 1 km. The JAXA ALOS PALSAR time series has greater coherence than our S1A time series due its longer, L-band wavelength, minimizing interference with vegetation.
We propose that the cumulative wastewater injection from multiple wells into shallow layers is causing pore pressure increase and in turn regional uplift. In downtown Cushing, none of the wastewater injection wells within 6 km are injecting into or deeper than the Arbuckle formation. In this same region, we are seeing the greatest localized regional uplift. Figure 14 shows that, although regional uplift is caused by the 2015 earthquakes, the ground continues to uplift into our final cumulative displacement image.
We attribute the 2014, 2015, and 2016 sequences to deep wastewater injection into the Arbuckle formation by local wastewater injection wells. When wastewater was injected into the Arbuckle, the pore pressure increase reduced the effective normal stress on an existing fault, inducing slip. Another possible explanation is the combined pore pressure increase from injection into the Arbuckle and shallower formations resulted in a pore pressure gradient, facilitating in situ fluid to flow into existing faults within the Arbuckle. The fluid flow reduced the frictional strength of the fault, inducing slip [31]. Figure 15 presents injection data for all UIC wells within 10 km of downtown Cushing. Data for the resulting 33 UIC injection wells dates back to 2011. Although we cannot comment on the source of the regional uplift from 2006 to 2011, we do see a significant increase in total injected volume per month starting in 2014 and peaking in 2015. From 2011 to 2018, only a single UIC well, API 3511923843, injects at high-rates. The period at which API 3511923843 injects at high-rates into the Arbuckle immediately precedes the 2014 earthquake sequence.  Our inversion suggests the depth of the 2016 Mw 5.0 earthquake is much shallower than that of the USGS model. The shallower source is in agreement with the Schoenball and Ellsworth [5] depth estimate of 3.625 km within the crystalline basement. Past studies [32][33][34] have attributed the discrepancy between seismic and geodetic source depths to the coarse approximation for the uppermost crust of the elastic halfspace model that does not consider possible rigidity contrasts. However, the observed MMI VII damage in downtown Cushing further supports a shallower source depth.

Conclusions
In our study, we characterize the impact of wastewater injection on the city of Cushing and its critical national infrastructure using DInSAR. In terms of long-term deformation, we observe 4 to 6 cm of regional uplift across the city over 17 months (Figure 16). We also successfully observe instantaneous, coseismic deformation from the 2016 Mw 5.0 earthquake. This is the best example of coseismic deformation associated with a moderate (5.0 ≤ M ≤ 5.9), shallow induced earthquake to date. Although the 2016 Mw 5.0 earthquake was too small to cause major infrastructure damage related to permanent displacement, a larger earthquake at this depth could have significant impacts. Following a magnitude 4.8-5.3 earthquake, the Oklahoma Department of Transportation (ODOT) mandates a 15-mile inspection radius from the epicenter of roads and bridges. Using critical infrastructure data in conjunction with high-resolution horizontal and vertical displacement maps derived from DInSAR can help identify those locations within such a large area that need to be prioritized for mitigation or remediation. Our source inversion demonstrates the potential for InSAR to facilitate rapid response efforts for shallow, moderate-sized earthquakes, especially in poorly seismically-instrumented areas. The phase coherence images overlaid with municipal boundaries ( Figure A1) show that the municipals are areas of consistent coherence. To avoid processing bias, we construct a time series using a reference region of 100 by 100 meters within the municipal boundary of Carney, OK, USA. The wrapped coseismic interferograms for the 2016 Mw 5.0 earthquake ( Figure A2), when compared to the unwrapped coseismic interferograms, can highlight unwrapping errors. No significant unwrapping errors are present in our coseismic unwrapped interferograms.