Abstract Observations that unequivocally link seismicity and wastewater injection are scarce. Here we show that wastewater injection in eastern Texas causes uplift, detectable in radar interferometric data up to >8 kilometers from the wells. Using measurements of uplift, reported injection data, and a poroelastic model, we computed the crustal strain and pore pressure. We infer that an increase of >1 megapascal in pore pressure in rocks with low compressibility triggers earthquakes, including the 4.8–moment magnitude event that occurred on 17 May 2012, the largest earthquake recorded in eastern Texas. Seismic activity increased even while injection rates declined, owing to diffusion of pore pressure from earlier periods with higher injection rates. Induced seismicity potential is suppressed where tight confining formations prevent pore pressure from propagating into crystalline basement rocks.

In recent years, the eastern and central United States have experienced a sharp increase in the number of earthquakes, with more than 1570 events with moment magnitudes (M w ) ≥ 3 between 2009 and 2015 (1–3). Many of these events occurred near injection disposal wells and were preceded by a high rate of fluid injection over a period of months to years, suggesting a link between seismicity and injection operations (1, 3–8). In general, earthquake hazard is proportional to the seismic rate; thus, the current increase in the seismic rate implies an elevated hazard in the central and eastern United States (9).

On 17 May 2012, the city of Timpson, Texas, experienced a M w 4.8 earthquake, the largest recorded event in the region (Fig. 1A). This event was preceded and followed by several earthquakes located on an inferred northwest-southeast–trending basement fault, including three with M w ≥ 4.0 over the following 16 months. Focal depths were shallow, ranging from 1.6 to 5.0 km, with the majority of the strain released between 3.5 and 5 km (5). Four high-volume class II disposal wells are located within ~10 km, and two lie directly above the site of the earthquakes (table S1). They dispose coproduced saline formation water from oil and gas production operations in the area by injecting into Lower Cretaceous limestones within the Sabine Uplift of eastern Texas (10). There are other injection wells in the vicinity, but none is closer than 7 km to the four studied wells (table S8 and fig. S8) or is a high-volume injector. The immediate area near the disposal wells and the earthquakes has limited oil and gas production. The four disposal wells began injection between 2005 and 2007 at a net average rate of 890,000 m3/year until mid-2012, when injection dropped to 720,000 m3/year (5).

Fig. 1 Study area and data sets. (A) Three overlapping frames from ALOS in ascending orbit (heading = 350°, incidence = 34.5°). Locations of the seismicity and focal mechanism of the Timpson sequence are shown (white dots). Latitude and longitude are shown on the y and x axes, respectively. (B) Contours of the LOS deformation velocity field obtained from multitemporal processing of the overlapped ALOS SAR data set (yr, year). The west wells (W1 and W2) and the east wells (E1 and E2) are indicated. (C) Time series of the volume of injected fluid for each of the wells shown in (B) (Mo, month).

The proximity of the earthquake clusters to the injection wells suggests a link between them (5, 11). As wastewater is injected into the disposal formation, it increases pore pressure within the connected hydrogeologic system. Over time, the pressure perturbation can spread to distances of many kilometers (12, 13). The increase in pore pressure caused by the injection of fluids decreases the effective normal stress on faults, bringing them closer to failure (14–16), as well as locally changing differential stress within the reservoir and surrounding rocks (17–19). Moreover, increased pore pressure can cause surface deformation (19), measurable using geodetic tools (20), which provides the possibility of documenting subsurface evolution from the surface. Among geodetic tools, interferometric synthetic aperture radar (InSAR) and the Global Positioning System (GPS) have been widely used to monitor surface deformation resulting from natural and anthropogenic processes.

We applied a multitemporal InSAR approach (21) to three overlapping sets of L-band SAR images (tables S2 and S3) acquired by the Advanced Land Observing Satellite (ALOS) over the Timpson area between 6 May 2007 and 14 November 2010 (Fig. 1A). High-quality interferograms were generated from these L-band data (fig. S1). We improved the signal-to-noise ratio of the measurements by estimating the linear velocity from time series obtained by inverting a large number of interferograms (22). Using velocity maps for each individual data set and the combined map (22), we found a line-of-sight (LOS) uplift rate of up to 3 mm/year over the area between the injection wells (Fig. 1B). We estimated that the rate of volume increase under the LOS velocity surface is 800,000 to 1,000,000 m3/year, assuming an elastic material with a Poisson’s ratio of 0.25 to 0.33 estimated from seismic velocities profiles, which is consistent with the net injected volume rate at the injection wells. To validate ALOS observations and evaluate the current state of ground deformation, we also applied the same processing scheme to eight C-band images acquired by the RADARSAT-2 satellite between 6 March and 21 August 2014 (fig. S2 and tables S4 and S5). Although the location of the zone of maximum deformation is consistent with the velocity map obtained from the L-band data, its spatial extent is broader, with a maximum uplift of ~5 mm over a ~6-month interval.

The two western wells (W1 and W2) inject (Fig. 1C) at a depth of 1800 m into Trinity Group formations (fig. S6), which are porous and permeable limestones overlain by the regionally extensive and much less permeable Ferry Lake Anhydrite (10). The east wells (E1 and E2) inject (Fig. 1C) into carbonate formations of the Washita Group at a depth of 900 m (fig. S6), which is stratigraphically above the Ferry Lake Anhydrite (10). We applied an inverse modeling scheme (23) to characterize the rate of volume strain. We discretized the volume beneath the half-space into rectangular prisms, 3 by 3 km in area and 0.2 km high, between the surface and a depth of 5 km. We assumed constant volume strain within each prism (22). Our optimum strain model accurately reproduced the observed deformation data (fig. S3) and calculated a total volume change of 700,000 ± 1600 m3/year, slightly lower than the injection rate. This discrepancy is likely due to diffusion of injected fluids into the surrounding rocks without generating any measurable deformation. We also found a maximum volume strain rate of ~5 × 10−6 yr−1, located at a depth of 0.8 to 1.1 km adjacent to wells E1 and E2 (Fig. 2A).

Fig. 2 Volumetric strain and pore pressure. (A) Distribution of the estimated volume strain rate (shaded grid cells). Colored circles show the timing of earthquakes with respect to the first event (day 0). The injection wells are shown by green bars. Contour lines show the surface deformation rate between 6 May 2007 and 14 November 2010. (B) Distribution of the cumulative pore pressure (shaded grid cells) between 2006 and 2013. Colored circles show the pore pressure increase at the location of the earthquakes.

The availability of geological profiles (fig. S6) and the distribution of hydraulic conductivity and Poisson’s ratios allowed us to characterize parameters of a poroelastic layered Earth model (table S6). Using this Earth model and the time series of injected fluid volume, we solved for the evolution of the pore pressure in the crust (22) (Fig. 2B). We identified two zones of maximum pore pressure at depths of ~0.85 and ~1.85 km near the east and west wells, respectively. The shallower zone of elevated pore pressure coincides with the zone of maximum volume strain. Higher pressures occur around the two west wells, where uplift has been negligible. The contribution from other injection wells to the south (figs. S9 and S10) is small and of second-order importance. The pore pressure associated with wells to the north is unlikely to influence the seismicity because these wells lie on the downthrown block of the Mount Enterprise fault system (10), which probably acts as a barrier to southward migration of the fluids.

To investigate the relationship between the pore pressure distribution associated with injection and the observed seismicity, we estimated the increase in pore pressure at the location of the 2012 earthquake sequence, in which the main events nucleated between 3.5 and 4.5 km depth (5). Overall, fluid injection caused a pore pressure increase of 0.5 to 1.5 MPa at the hypocentral depth. Pressure changes of this magnitude are known to trigger earthquakes elsewhere (24). In the context of Mohr circle stress analysis, a localized increase in pore pressure shifts the circle to the left (i.e., reduces the effective normal stress) and changes its radius because of poroelastic strain (i.e., increases the differential stress), whereas a spatially uniform pore pressure increase only shifts the circle to the left until it touches the failure envelope (25). Given the historical lack of large earthquakes in the region and the 5-year delay between the initiation of injection and the first large event, we suggest that a decrease in effective normal stress (due to a homogeneous increase in pore pressure) triggered seismicity. The initiating seismicity and associated stress change may have transiently enhanced the permeability (26–29), increasing the pore pressure at the location of the deeper events, which in turn promoted additional earthquakes.

We investigated injection from two pairs of wells that began injecting at about the same time, disposing approximately the same volumes of wastewater. The main differences between them are the depth of injection and the presence of an impermeable barrier below the shallower east wells that blocks fluid and pressure from reaching deeper formations. The deeper west wells are associated with the 2012 Timpson earthquake sequence, whereas no detected seismicity occurred near the east wells. This observation highlights the importance of hydrogeology for the consequences of fluid injection.

The extent to which induced pore pressure change occurs is an important parameter for accurately estimating seismic hazard. From a regulatory perspective, however, constraining this parameter is not trivial, given its complex relationship with local hydrogeology (30). Using deformation data, we can put a lower bound on the extent of the rock volume change caused by a pore pressure increase. Measurable uplift more than 8 km from the east wells demonstrates the long reach of pressure perturbations, which has been inferred in other studies (12, 13, 15).

Studies of potentially induced earthquakes suggest that the majority of seismicity occurs within basement rocks, even though most of the injection is done in more shallow sedimentary layers (2–5, 7, 31, 32). Although pore pressure has increased adjacent to both the east and west wells, little surface deformation was detected in the vicinity of the west wells, where the seismicity occurred (Fig. 2). This is likely due to lower rock compressibility near the west wells compared with near the east wells, a feature we did not account for in solving for the pore pressure evolution, in addition to the greater injection depth. Moreover, injection in the east wells is done in a shallow layer, typically characterized by velocity-strengthening frictional properties (33); thus, pore pressure changes are less likely to initiate seismic rupture at this location. The uplift signal also has an asymmetric pattern, which we cannot explain with standard models of radial pore pressure diffusion in a homogeneous medium. Our model shows that pore pressure propagates downward below the west wells with a delay and that the period of elevated seismicity between 2010 and 2014 coincides with the period of the maximum pore pressure change of 1 to 2.5 MPa at the average depth of 3 km (Fig. 3A). Though all events coincide with the period of pore pressure increase in the focal zone, the onset of the main sequence in May 2012 corresponds to pressures of about 1 MPa reaching the focal zone (Fig. 3B).

Fig. 3 Pore pressure time series. (A) Colored lines show time series of pore pressure at various depths. The error bars show 1-standard-deviation uncertainties and were obtained through bootstrapping. Vertical black lines show the time and magnitude of earthquakes (Eq.) in the Timpson sequence (table S7). (B) Close-up view of the period 2010–2015 from (A). The time series of the rate of injected fluid at the west wells is superimposed in purple.

The frequency and magnitude of seismic events reached a climax between May 2012 and September 2013, when more than 80% of events occurred, including four M w ≥ 4 earthquakes. This period of elevated seismic activity followed a rapid decline in injection at the west wells (Fig. 3B). Thus, the timing of seismicity may be determined by pore pressure diffusion to the depths of the earthquakes. Additional contributions to the stresses may also promote seismicity, such as from the poroelastic stresses near the well that accompany the decrease in injection. For an optimally oriented strike-slip fault with the fault-normal direction radial to the injection site, a sudden decline in injection rate relaxes the compressive stress (34). All of the seismic events discussed here occurred on a fault with the fault-normal direction oriented at N60°E, radial to the west injection wells.

Our coupled flow and poroelastic model allows us to predict the future pore pressure distribution after injection ends (Fig. 3A). As injection ends in the west wells by about 2016, the model shows the pressure decreasing approximately exponentially. However, the decay rate is fastest in the most permeable formations and at the depths where injection occurs. At the east wells, 10 years after shut-off, the pore pressure remains close to its maximum level. For the west wells, only 5 years after a hypothesized shut-off, pore pressure drops to less than 10% of its maximum value. This observation has direct implications for future injection operations and seismic hazard estimations. Changes in the seismicity rate are a function of changes in Coulomb stress and background stress (35). Our study shows that the background stress is characterized by a relaxation time that depends on both the injection history and hydrogeological properties. Therefore, injection history at a given site may modify future estimates of the seismic hazard.

Better quantification of the evolution of the stress and pore pressure in the crust is vital to forecasting fault reactivation (9, 36). Despite improvements in seismic monitoring capacity and the resulting decrease in the magnitude detection threshold (37, 38), observations of the in situ pore pressure and stress field remain elusive because of the scarcity of deformation observations and limited integration of those observations with physical models. This work highlights the value of monitoring surface deformation, in particular by using advanced remote sensing techniques, to understand the evolution of pore pressure and stress at depth. The ability to measure crustal stress evolution affords a proactive approach to managing the hazard associated with fluid injection. Observation of the time-dependent stress field permits the construction of temporally variable statistical frameworks (34), which are useful for earthquake operational forecasting (39). The key to successful operational earthquake hazard assessment is being able to continuously update information about the probability of a future earthquake, which can be achieved by using data and models such as those presented in this study. Geodetic monitoring and modeling schemes are valuable components of efforts to mitigate induced seismic hazard.

SUPPLEMENTARY MATERIALS www.sciencemag.org/content/353/6306/1416/suppl/DC1 Materials and Methods Figs. S1 to S10 Tables S1 to S8 References (40–57)