WEBVTT Kind: captions Language: en 00:00:02.120 --> 00:00:06.320 [background conversations] 00:00:07.120 --> 00:00:08.680 Okay, folks. I’m going to get started. 00:00:08.680 --> 00:00:12.860 I’m filling in for Jack today, as he will be our presenter. 00:00:12.860 --> 00:00:14.500 Quick announcement. For next week, 00:00:14.500 --> 00:00:18.260 we’ll have Lise Retailleau from Stanford, and she’s going to 00:00:18.260 --> 00:00:20.880 talk about core mantle boundary imaging. 00:00:21.420 --> 00:00:23.760 But today, it’s my pleasure to introduce Jack. 00:00:23.769 --> 00:00:28.029 He received his undergraduate degree from CU-Boulder, followed by 00:00:28.029 --> 00:00:31.739 a master’s at Colorado School of Mines, and then he joined Roland Horne’s 00:00:31.739 --> 00:00:38.300 group at Stanford for a Ph.D. in geothermal energy-related topics. 00:00:38.300 --> 00:00:40.550 He only joined the USGS as a Mendenhall fellow 00:00:40.550 --> 00:00:45.660 earlier this year and proposed to look at improving the quantification 00:00:45.660 --> 00:00:48.800 of the induced seismicity hazard. 00:00:48.800 --> 00:00:51.440 And this is actually exactly what he’s going to propose – 00:00:51.440 --> 00:00:54.500 what he’s going to talk about today. 00:00:54.500 --> 00:00:57.340 And it may make him the first Mendenhall to actually deliver 00:00:57.340 --> 00:01:01.530 on the proposed things exactly, and that within the first year. 00:01:01.530 --> 00:01:04.380 So without further ado, Jack. 00:01:06.340 --> 00:01:09.039 - Okay, thanks for the introduction, Ole. 00:01:09.040 --> 00:01:11.600 And good morning, everyone. 00:01:12.780 --> 00:01:17.660 So as Ole mentioned, today I’m going to be speaking about induced seismicity. 00:01:19.000 --> 00:01:22.690 In fact, I did get introduced to the topic of induced seismicity 00:01:22.690 --> 00:01:25.210 through the geothermal community. 00:01:25.210 --> 00:01:28.450 There are actually quite a few interesting case studies of 00:01:28.450 --> 00:01:32.590 induced seismicity at geothermal sites. Probably the most famous one 00:01:32.590 --> 00:01:39.240 is the Basel EGS project, which, after experiencing a few 00:01:39.240 --> 00:01:45.100 magnitude 3 earthquakes during stimulation, was actually cancelled. 00:01:45.100 --> 00:01:50.920 And it’s really interesting to see the differences between, you know, 00:01:50.920 --> 00:01:53.990 how that has actually influenced the geothermal community 00:01:53.990 --> 00:01:58.210 over the last few years, where now every major geothermal project 00:01:58.210 --> 00:02:04.939 really has a strong focus on understanding induced seismicity. 00:02:04.939 --> 00:02:08.209 But later on, I got introduced to the Stanford Center for Induced 00:02:08.209 --> 00:02:10.009 and Triggered Seismicity. And of course, 00:02:10.009 --> 00:02:14.440 there, over the last few years, the main issue is to try and 00:02:14.440 --> 00:02:19.160 understand Oklahoma’s earthquake problem. 00:02:19.160 --> 00:02:23.450 So in fact, in 2015, there were nearly 900 magnitude 3 and larger 00:02:23.450 --> 00:02:28.010 earthquakes in Oklahoma. And last year, there were actually 00:02:28.010 --> 00:02:31.700 the most magnitude 5 and larger earthquakes – there were three of them – 00:02:31.700 --> 00:02:37.320 signifying that the hazard is still likely elevated in that area of the country. 00:02:37.960 --> 00:02:44.720 Now, since the earthquake rates in Oklahoma were changing really rapidly, 00:02:44.720 --> 00:02:50.680 that also gave rise to, like I said, an elevated hazard. 00:02:50.680 --> 00:02:57.260 And because things were changing so quickly, the USGS 00:02:57.260 --> 00:03:01.040 wanted to try and get on top of that by introducing kind of these 00:03:01.040 --> 00:03:06.900 one-year hazard models, which are different from 00:03:06.900 --> 00:03:12.100 the hazard models that have been introduced previously, 00:03:12.100 --> 00:03:15.740 where they’re kind of looking at longer-term 50-year forecasts. 00:03:15.740 --> 00:03:18.110 Things were changing quickly, so they opted to do these 00:03:18.110 --> 00:03:23.340 kind of one-year models. But up this point, although they 00:03:23.340 --> 00:03:27.430 included the impact of induced earthquakes in the sense of the statistical 00:03:27.430 --> 00:03:33.650 observations of the seismicity, they have yet to incorporate the operational 00:03:33.650 --> 00:03:37.860 parameters – the injection data, the locations of where the wells are. 00:03:37.860 --> 00:03:44.870 And since the hazard is likely really being driven by the saltwater disposal wells, 00:03:44.870 --> 00:03:47.660 it makes sense to try and incorporate that type of information. 00:03:47.660 --> 00:03:51.340 And really, the goal of my research has been to try and develop 00:03:51.340 --> 00:03:55.370 a physics-based understanding of the link between fluid injection 00:03:55.370 --> 00:03:58.980 and the mechanical behavior of fractures and faults 00:03:58.980 --> 00:04:03.870 in the subsurface that are capable of hosting these earthquakes. 00:04:03.870 --> 00:04:06.770 So it fits in well with the goal of trying to 00:04:06.770 --> 00:04:10.660 incorporate the operational data into a hazard framework. 00:04:11.640 --> 00:04:14.960 So in my talk today, I’m going to be speaking about 00:04:14.960 --> 00:04:20.040 two studies that I’ve performed. The first is some work 00:04:20.040 --> 00:04:24.560 I did during my Ph.D. And it’s kind of a case study 00:04:24.560 --> 00:04:30.320 looking at the magnitude 5.7 Prague earthquake that occurred in 2011. 00:04:30.320 --> 00:04:34.340 And hopefully I’ll get through that relatively quickly because 00:04:34.340 --> 00:04:37.490 I’m really excited to share some of the more recent work 00:04:37.490 --> 00:04:41.000 that I’ve been doing here since starting at the USGS back in January. 00:04:41.000 --> 00:04:43.990 And that’ll be looking at forecasting rates of 00:04:43.990 --> 00:04:47.220 induced seismicity using injection data. 00:04:47.760 --> 00:04:50.300 All right, so we’ll jump right into it. 00:04:50.310 --> 00:04:55.650 The motivation for the study of the Prague sequence was that, for one, 00:04:55.650 --> 00:04:58.430 it was the largest earthquake that had occurred in Oklahoma 00:04:58.430 --> 00:05:00.410 at the time I was doing the study, so it was really important 00:05:00.410 --> 00:05:04.060 to understand the processes involved in that sequence. 00:05:04.060 --> 00:05:08.430 But also, more broadly speaking, you know, the hazard is really 00:05:08.430 --> 00:05:13.120 dominated by these large magnitude 5 and greater earthquakes. 00:05:13.120 --> 00:05:18.060 And when you think about these larger earthquakes, you really have to consider 00:05:18.070 --> 00:05:23.600 finite source considerations. So things like the geometry of the fault, 00:05:23.600 --> 00:05:28.030 the state of stress in that particular location, and then also, in the context of 00:05:28.030 --> 00:05:31.000 induced seismicity, we have to think about how flow is occurring 00:05:31.000 --> 00:05:35.310 through the aquifer and down and along these large basement faults 00:05:35.310 --> 00:05:38.500 that can be on the order of 5 or 10 kilometers long. 00:05:38.500 --> 00:05:42.610 So really the goal of this study was to try and understand 00:05:42.610 --> 00:05:47.140 that fault zone hydrology part of the problem. 00:05:48.090 --> 00:05:52.500 So the Prague sequence is interesting for a lot of reasons. 00:05:52.500 --> 00:05:58.389 But I think it’s really well-suited for this study where we’re trying to 00:05:58.389 --> 00:06:03.880 look at hydraulic properties of these large basement faults. 00:06:03.880 --> 00:06:10.840 And really, the – so the sequence occurred in early November of 2011. 00:06:10.840 --> 00:06:15.500 And what actually happened was, there was a magnitude 4.8 foreshock 00:06:15.500 --> 00:06:19.760 that occurred kind of along this main section of the Wilzetta Fault. 00:06:19.760 --> 00:06:24.071 And one day later, the magnitude 5.7 earthquake 00:06:24.080 --> 00:06:28.620 occurred kind of along this splay fault. 00:06:28.620 --> 00:06:31.360 And so really the key observation driving this study was that 00:06:31.360 --> 00:06:35.430 there was this one-day delay between the foreshock and the main shock. 00:06:35.430 --> 00:06:39.850 And really, based on some prior work that had been done, 00:06:39.850 --> 00:06:43.620 it had been shown that the geometry of these faults were such that the static 00:06:43.620 --> 00:06:51.600 stress change caused by that foreshock may have triggered the main earthquake. 00:06:51.600 --> 00:06:54.400 But the fact that there was this one-day delay 00:06:54.400 --> 00:06:57.540 implies some sort of transient process. 00:06:57.540 --> 00:07:01.560 And so what I proposed is that it could have been influenced by a 00:07:01.560 --> 00:07:07.440 hydromechanical response, but we also looked at rate-and-state friction effects. 00:07:09.150 --> 00:07:15.000 So fortunately, there had been a lot of work done to learn about this sequence, 00:07:15.000 --> 00:07:17.010 and so we have a fairly good understanding of 00:07:17.010 --> 00:07:19.919 what the fault geometry actually looked like, 00:07:19.919 --> 00:07:22.370 which is – which is really important for this problem. 00:07:22.370 --> 00:07:27.890 And so we were able to kind of identify the fault locations and their geometry 00:07:27.890 --> 00:07:32.780 based on focal mechanisms and the distribution of aftershock locations. 00:07:32.780 --> 00:07:37.780 So really, the two earthquakes that we’re considering in this study 00:07:37.780 --> 00:07:42.020 occurred on more or less vertical strike-slip faults, 00:07:42.020 --> 00:07:45.040 and they were located at about 5 kilometers’ depth. 00:07:46.240 --> 00:07:50.419 We also had access to some studies that had been done through the 00:07:50.419 --> 00:07:52.610 Stanford Center for Induced and Triggered Seismicity. 00:07:52.610 --> 00:07:58.530 So Rall Walsh and Mark Zoback performed a stress inversion study 00:07:58.530 --> 00:08:01.450 using focal mechanisms from nearby earthquakes. 00:08:01.450 --> 00:08:04.570 And so that gave us good information about the local state of stress, 00:08:04.570 --> 00:08:07.720 both in terms of the stress orientations and their magnitudes at depth. 00:08:07.720 --> 00:08:11.800 So we took all of that into account during the modeling. 00:08:12.660 --> 00:08:16.720 So for this study, we used a modeling software called CFRAC, 00:08:16.720 --> 00:08:19.370 which was introduced originally by Mark McClure, 00:08:19.370 --> 00:08:23.340 and then I also continued development during my Ph.D. 00:08:23.340 --> 00:08:29.120 Essentially, it does quasi-dynamic earthquake rupture, and it does 00:08:29.120 --> 00:08:32.729 the fault mechanics problem with a boundary element method. 00:08:32.729 --> 00:08:35.760 We can consider transient friction through a 00:08:35.760 --> 00:08:39.300 rate-and-state friction framework. 00:08:39.300 --> 00:08:42.500 And the really unique thing about CFRAC is that it does really tight 00:08:42.500 --> 00:08:47.130 coupling between fluid flow along the faults and the geomechanical problem. 00:08:47.130 --> 00:08:50.880 So the hydraulic properties of the faults can be functions of 00:08:50.880 --> 00:08:55.010 effective stress or slip and things like that. 00:08:55.680 --> 00:08:58.880 So, like I said, we set up a model of this fault system, 00:08:58.880 --> 00:09:06.080 and we kind of approximated it as a series of two vertical faults. 00:09:06.640 --> 00:09:11.320 And I want to kind of point out that I’m not actually considering 00:09:11.330 --> 00:09:15.690 fluid injection in this study. So I’m only considering the 00:09:15.690 --> 00:09:21.110 hydromechanical interaction between these two faults during the sequence. 00:09:21.110 --> 00:09:22.980 So there’s no flow in the surrounding rock, 00:09:22.980 --> 00:09:26.620 but there is fluid flow along these fault structures. 00:09:29.180 --> 00:09:34.540 We simulated the foreshock as a magnitude 4.8 earthquake 00:09:34.540 --> 00:09:39.420 that resulted from a uniform stress drop of about 1.7 megapascals 00:09:39.420 --> 00:09:44.360 along that upper fault in the upper right. 00:09:44.360 --> 00:09:50.339 And the state of stress resulting from that earthquake 00:09:50.339 --> 00:09:52.860 is shown here in Cartesian coordinates. 00:09:52.860 --> 00:09:56.431 If you rotate and resolve those stresses into their shear and 00:09:56.431 --> 00:09:59.800 normal components, acting on the fault that 00:09:59.800 --> 00:10:04.370 hosted the main earthquake, this is kind of what it looks like. 00:10:04.370 --> 00:10:10.560 So what you’ll notice is that there’s kind of a stress drop region here 00:10:10.560 --> 00:10:14.500 where you would expect it, but also, you got a zone of increased shear here 00:10:14.500 --> 00:10:19.360 ahead of that Fault A crack tip, and also a reduction in normal stress. 00:10:19.360 --> 00:10:24.230 And so both of those processes actually encourage earthquake rupture. 00:10:24.920 --> 00:10:28.760 So taking a little bit closer look as to what the state of stress 00:10:28.760 --> 00:10:33.880 along that fault that hosted the main earthquake actually was in the model, 00:10:33.880 --> 00:10:39.100 the changes in shear normal stress were on the order of about 2 megapascals. 00:10:39.100 --> 00:10:44.360 And see that the stress perturbations fell off over distances of about a kilometer. 00:10:44.370 --> 00:10:47.140 And also just point out that there’s this unclamping, 00:10:47.140 --> 00:10:52.100 or a reduction in normal stress, over parts of the fault. 00:10:52.460 --> 00:10:57.540 Again, both of those would encourage earthquake nucleation. 00:10:57.550 --> 00:10:59.560 Take a minute to go through some of the hydromechanical 00:10:59.560 --> 00:11:02.380 coupling that we’re considering in the model. 00:11:02.380 --> 00:11:05.680 So what I’ve just written here is the mass balance equation 00:11:05.680 --> 00:11:10.260 for single-phase flow – one-dimensional flow along a fault. 00:11:10.260 --> 00:11:15.290 And the flux term is related to the fault transmissivity, or permeability. 00:11:15.290 --> 00:11:20.600 And the storage term is related to the fault void aperture, or fault porosity. 00:11:21.520 --> 00:11:29.120 And we’re not actually solving the full poroelastic solution to this – 00:11:29.120 --> 00:11:31.750 to this problem, but we do have the flavor of hydromechanical 00:11:31.750 --> 00:11:37.930 or poroelastic effects because we’ve introduced a constitutive relationship 00:11:37.930 --> 00:11:44.470 that relates changes in effective stress to changes in the porosity of the fault. 00:11:44.470 --> 00:11:48.700 So the idea is that this fault is somehow compliant. 00:11:49.649 --> 00:11:54.160 And so, at early times, right as that fault is feeling the static stress 00:11:54.160 --> 00:12:00.680 change by the foreshock, it’ll respond in a kind of undrained behavior. 00:12:00.680 --> 00:12:06.700 And so, in the undrained limit, the flux term vanishes, and you can 00:12:06.700 --> 00:12:10.779 introduce fault compressibility and water compressibility to arrive at 00:12:10.779 --> 00:12:15.550 kind of an approximation for something analogous to the Skempton’s coefficient, 00:12:15.550 --> 00:12:21.220 which just gives you a sense of what the – for a given change in 00:12:21.220 --> 00:12:25.940 normal stress along this fault, what will the fluid pressure change be? 00:12:25.940 --> 00:12:27.700 And it’s related to, like I said, 00:12:27.700 --> 00:12:31.060 the fault compressibility and water compressibility. 00:12:32.780 --> 00:12:36.470 So now we want to think about how that will affect the state of 00:12:36.470 --> 00:12:40.490 stress along the fault and how it evolves over time. 00:12:40.490 --> 00:12:44.980 And so want to look at how the Coulomb stress evolves 00:12:44.980 --> 00:12:47.870 along this fault. And so Coulomb stress obviously 00:12:47.870 --> 00:12:52.740 depends on changes in shear stress, normal stress, and fluid pressure. 00:12:52.740 --> 00:12:57.660 And I think there is a question as to how compliant these faults are in nature. 00:12:57.660 --> 00:13:00.709 And so that’s one of the things we’re trying to address with this study. 00:13:00.709 --> 00:13:04.760 So for a very compliant fault, like I’m showing here on the left, 00:13:04.760 --> 00:13:10.690 essentially, as that unclamping occurs, and the volume of – the pore volume of 00:13:10.690 --> 00:13:17.660 the fault increases, the pore pressure in the fault drops instantaneously. 00:13:17.660 --> 00:13:22.580 And what this has the effect of doing is essentially cancelling out 00:13:22.589 --> 00:13:26.880 that decreased normal stress, at least temporarily. 00:13:26.880 --> 00:13:31.339 And so then the Coulomb stress is just dominated by the shear stress effects. 00:13:31.339 --> 00:13:37.040 But if the fault is able to behave slightly stiffer, then that 00:13:37.040 --> 00:13:41.420 fluid pressure change only partially mitigates the normal stress response. 00:13:41.420 --> 00:13:45.490 And then you get slightly larger total Coulomb stress change 00:13:45.490 --> 00:13:49.529 in the undrained limit. But of course, once pressure gradients 00:13:49.529 --> 00:13:54.630 in the fault set up, that will encourage fluid flow to occur along the fault. 00:13:54.630 --> 00:13:57.760 Because the fault is permeable, that fluid pressure wants to 00:13:57.760 --> 00:14:00.950 equilibrate back to the initial condition over time. 00:14:00.960 --> 00:14:04.820 And so that actually looks like a transient loading mechanism. 00:14:07.340 --> 00:14:10.640 So I mentioned that we’re using rate-and-state friction in this model. 00:14:10.640 --> 00:14:13.790 So that allows us to kind of keep track of how the 00:14:13.790 --> 00:14:17.220 sliding velocity along the fault evolves over time. 00:14:17.220 --> 00:14:22.240 So what I’m showing here is just time since the beginning – 00:14:22.240 --> 00:14:26.550 since the end of the foreshock. And on the Y axis is the 00:14:26.550 --> 00:14:32.600 sliding velocity at the nucleation site for the main earthquake. 00:14:32.600 --> 00:14:35.160 And essentially, it starts out at very low 00:14:35.160 --> 00:14:39.580 sliding velocities – more or less a locked fault. 00:14:39.580 --> 00:14:45.740 And then, as time progresses, this kind of behavior can kind of 00:14:45.740 --> 00:14:48.530 be characterized by three – more or less three phases. 00:14:48.530 --> 00:14:53.030 So the first is just the initial undrained response to the Coulomb stress change. 00:14:53.030 --> 00:14:55.769 And then there’s kind of a transient period where both 00:14:55.769 --> 00:14:59.430 fluid flow and friction evolution are influencing the problem. 00:14:59.430 --> 00:15:02.620 And then finally, you have this rapid acceleration 00:15:02.620 --> 00:15:07.260 towards kind of sustained rupture propagation. 00:15:07.260 --> 00:15:10.980 So when the fault reaches kind of sliding velocities on the order of 00:15:10.980 --> 00:15:16.180 about a meter per second, then the rupture is sustained into an earthquake. 00:15:19.220 --> 00:15:23.760 So we were able to just vary properties of the fault to see 00:15:23.760 --> 00:15:30.560 how this might influence the problem. And again, we’re trying to shoot for 00:15:30.560 --> 00:15:33.700 a timing of the main earthquake of about 20 hours. 00:15:33.709 --> 00:15:38.529 So that’s shown by this vertical dashed line – the black line. 00:15:38.529 --> 00:15:43.940 And, in this case, what I’m showing is I’ve kept the fault transmissivity 00:15:43.940 --> 00:15:47.380 constant, and I’ve varied the compliance of the fault. 00:15:47.390 --> 00:15:53.531 And what you’ll see is that, for a very stiff fault, the Coulomb stress 00:15:53.531 --> 00:15:57.290 change in the undrained limit is enough to kind of kick that 00:15:57.290 --> 00:16:00.981 rupture into the nucleation phase. And we tend to underestimate 00:16:00.981 --> 00:16:04.820 the timing to instability by a couple of orders of magnitude. 00:16:04.820 --> 00:16:08.640 But if we allow that fault to behave a little bit more compliantly, 00:16:08.640 --> 00:16:12.840 we can get closer to the actual timing of the rupture. 00:16:16.040 --> 00:16:18.960 So here’s another case where now I’ve kept the fault compliance, 00:16:18.970 --> 00:16:21.360 which is, in this case, relatively high, constant. 00:16:21.360 --> 00:16:25.120 And now I’m varying the fault permeability or transmissivity. 00:16:25.120 --> 00:16:31.410 And so, for a very high transmissivity fault, even if the fault is behaving 00:16:31.410 --> 00:16:37.899 relatively compliantly, that fluid flow is occurring quickly enough to kick that 00:16:37.899 --> 00:16:41.110 rupture into nucleation early. And, again, we’re underestimating 00:16:41.110 --> 00:16:46.800 the timing, but for relatively lower transmissivity faults, 00:16:46.800 --> 00:16:50.340 we were able to kind of get the right timing. 00:16:50.340 --> 00:16:55.120 So kind of our takeaways for this study is, our goal was to 00:16:55.120 --> 00:17:01.320 try and constrain some of the in-situ properties for these large basement faults 00:17:01.320 --> 00:17:04.919 using a numerical model. And what the results are showing 00:17:04.919 --> 00:17:10.110 is that it is possible to kind of explain this behavior with these 00:17:10.110 --> 00:17:13.920 hydromechanical coupling effects, but ultimately, there’s a lot of 00:17:13.920 --> 00:17:16.600 non-unique coupling going on here between 00:17:16.600 --> 00:17:20.010 both the frictional properties and the hydromechanical properties. 00:17:20.010 --> 00:17:24.390 And so we weren’t able to hone in on one unique solution. 00:17:24.390 --> 00:17:29.570 But nonetheless, if we take the hydromechanical effects to – 00:17:29.570 --> 00:17:33.500 you know, to actually be representative of what’s happening in situ, 00:17:33.500 --> 00:17:37.610 then it did give us some constraints on the hydraulic properties of the fault. 00:17:37.610 --> 00:17:40.400 We think they’re behaving relatively compliantly. 00:17:40.400 --> 00:17:43.240 So Skempton’s coefficient kind of near 1. 00:17:44.360 --> 00:17:49.600 All right, so like I said, that was kind of a quick run-through of that study. 00:17:49.610 --> 00:17:57.600 And, you know, I think it gives a sense of more – a very detailed study 00:17:57.600 --> 00:18:04.549 of a specific earthquake sequence. It required a lot of data – you know, 00:18:04.549 --> 00:18:10.340 a lot of knowledge about the state of stress in that particular location. 00:18:10.340 --> 00:18:14.640 It required quite a few modeling assumptions and also some 00:18:14.640 --> 00:18:18.610 kind of detailed interpretation. And from a modeler’s perspective, 00:18:18.610 --> 00:18:22.540 I would kind of argue that it’s really intractable to try and take that 00:18:22.540 --> 00:18:26.960 type of modeling approach at that level of detail when we’re thinking about 00:18:26.960 --> 00:18:33.790 trying to evaluate the hazard at a broader, maybe statewide-type scale. 00:18:33.790 --> 00:18:37.200 So for the second half of the talk, I’m going to be discussing some of the 00:18:37.200 --> 00:18:39.720 work that I’ve done over the last six months to, you know, 00:18:39.730 --> 00:18:46.210 still have the flavor of physics in my modeling approach, but I think 00:18:46.210 --> 00:18:53.480 you’ll see that the model I’ll talk about here is going to be a lot more practical, 00:18:53.480 --> 00:18:58.650 especially with application to understanding seismic hazard 00:18:58.650 --> 00:19:00.980 over a broader scale. 00:19:01.690 --> 00:19:05.860 So for those of you who aren’t totally keyed in on what’s happening 00:19:05.860 --> 00:19:08.770 in Oklahoma, this is kind of oil-and-gas country. 00:19:08.770 --> 00:19:11.360 So what I’m showing on the right is just a typical well pad. 00:19:11.360 --> 00:19:16.530 You might have a couple of production wells on each pad. 00:19:16.530 --> 00:19:19.240 And then you can see that this is – you know, there’s quite a few 00:19:19.240 --> 00:19:23.100 of these well pads. This is the town of Cherokee, Oklahoma. 00:19:23.100 --> 00:19:26.179 And if you zoom out to a broader scale, 00:19:26.179 --> 00:19:32.220 you’ll see that these wells just are littering this entire region of the state. 00:19:32.220 --> 00:19:37.330 And it turns out that the formations that they’re producing oil and gas 00:19:37.330 --> 00:19:41.470 from in this part of the – of the country come with a high water cut. 00:19:41.470 --> 00:19:47.040 And so there’s a lot of wastewater – saltwater to be disposed of. 00:19:48.080 --> 00:19:54.600 And in Oklahoma and Kansas, the primary target for 00:19:54.600 --> 00:19:58.460 wastewater disposal is the Arbuckle Aquifer. 00:19:58.460 --> 00:20:05.600 And it’s a basal aquifer that exists really extensively throughout the area. 00:20:05.600 --> 00:20:10.740 Like I said, it exists right above the Precambrian basement rock. 00:20:12.560 --> 00:20:16.320 And it’s actually a really good aquifer, it turns out, for wastewater disposal. 00:20:16.320 --> 00:20:22.250 It’s relatively distant from any of the shallower freshwater aquifers. 00:20:22.250 --> 00:20:27.080 It’s highly permeable. And it’s already very salty. 00:20:27.080 --> 00:20:30.130 And so, just operationally, it’s a good target. 00:20:30.130 --> 00:20:34.919 But like I said, it exists right above the basement, which suggests that there is 00:20:34.919 --> 00:20:43.100 some possibility to connect with larger faults that exist in the basement rock. 00:20:43.100 --> 00:20:47.140 So I’ve worked really hard with Andy Barbour and Justin Rubinstein 00:20:47.140 --> 00:20:50.380 to develop this kind of saltwater disposal database. 00:20:51.140 --> 00:20:54.880 What you’ll see here are 893 injection wells 00:20:54.880 --> 00:20:57.780 that are completed in the Arbuckle Aquifer. 00:20:59.060 --> 00:21:04.560 We’ve obtained the data from various sources, so the 00:21:04.570 --> 00:21:08.380 oil and gas is kind of regulated at the state level in most parts of the country. 00:21:08.380 --> 00:21:12.340 And so the injection data is reported by the Oklahoma 00:21:12.340 --> 00:21:15.170 Corporation Commission and the Kansas Corporation Commission, 00:21:15.170 --> 00:21:19.380 and then there’s also a small subset of wells that are regulated by the EPA. 00:21:19.380 --> 00:21:24.049 So all this data – we have injection rate data 00:21:24.049 --> 00:21:27.880 for all these wells going back to 1995. 00:21:27.880 --> 00:21:35.320 And primarily, for most of the time period that we’re considering, 00:21:35.320 --> 00:21:38.730 we have kind of monthly resolution of the injection data, 00:21:38.730 --> 00:21:42.680 especially during the period of highest earthquake activity. 00:21:42.680 --> 00:21:47.280 And, again, these wells are kind of spaced on the order of about 00:21:47.280 --> 00:21:54.200 2 to 5 kilometers. So really kind of a dense network of these injection wells. 00:21:57.280 --> 00:22:00.600 There’s not a lot of good pressure data, so we’re 00:22:00.600 --> 00:22:05.120 mostly working with injection rates. And I just want to point out that, 00:22:05.120 --> 00:22:09.190 you know, although there’s a lot of practical interest 00:22:09.190 --> 00:22:12.470 in understanding the link between injection and the earthquakes 00:22:12.470 --> 00:22:16.350 that have happened in Oklahoma, it’s also a great opportunity from a 00:22:16.350 --> 00:22:20.740 scientific perspective to learn more about how the sequence has evolved. 00:22:20.740 --> 00:22:28.050 You know, imagine we now have this injection data on a 2- to 5-kilometer 00:22:28.050 --> 00:22:33.950 resolution, monthly injection rates for a period of 20 years, over a huge region of 00:22:33.950 --> 00:22:39.690 about 300 kilometers in diameter there. So really compared to studies 00:22:39.690 --> 00:22:43.200 of natural earthquakes, we have a really remarkable 00:22:43.200 --> 00:22:45.500 data set here to work with in terms of 00:22:45.510 --> 00:22:48.640 understanding the connection between stressing conditions 00:22:48.640 --> 00:22:52.460 along faults and how seismicity responds to that. 00:22:53.260 --> 00:22:56.080 So just taking a look at how the earthquake sequence 00:22:56.080 --> 00:22:58.700 has evolved in Oklahoma and Kansas. 00:22:58.710 --> 00:23:03.150 This is some really fascinating work that’s been done by Martin Schoenball 00:23:03.150 --> 00:23:08.780 and Bill Ellsworth that just appeared in SRL last month, actually. 00:23:08.780 --> 00:23:12.460 And what they did was perform high-precision earthquake relocations 00:23:12.460 --> 00:23:16.910 of the seismicity throughout the area. And what you can see is all these 00:23:16.910 --> 00:23:21.280 linear features that can be interpreted as vertical strike-slip faults. 00:23:21.280 --> 00:23:26.000 They’re more or less optimally oriented for failure in the – 00:23:26.000 --> 00:23:30.480 in the tectonic stress environment in this region of the country. 00:23:30.480 --> 00:23:33.300 And just to point out some of the areas, like, this is the 00:23:33.300 --> 00:23:40.860 magnitude 5.8 Pawnee earthquake. These are the Cushing sequences. 00:23:42.160 --> 00:23:47.380 Here’s another part of the state. This is the Fairview earthquake. 00:23:47.380 --> 00:23:55.260 And, again, what I think this really points out is that these faults 00:23:55.260 --> 00:23:58.159 that are existing in the basement, predominantly between 00:23:58.159 --> 00:24:01.850 1 and 5 kilometers below the basement-Arbuckle contact 00:24:01.850 --> 00:24:03.980 are more or less ubiquitous throughout the area. 00:24:03.980 --> 00:24:08.670 And so that’s kind of the working mentality in terms of 00:24:08.670 --> 00:24:11.799 understanding what the faults actually look like in this area. 00:24:11.799 --> 00:24:14.980 So this study is really invaluable for understanding, 00:24:14.980 --> 00:24:20.200 like I said, what these faults really look like in our study area. 00:24:20.200 --> 00:24:23.200 And so taking a look at how the earthquake sequence, again, 00:24:23.200 --> 00:24:29.620 has evolved in concert with disposal throughout the region. 00:24:30.660 --> 00:24:33.300 Like I said, we have this data going back to 1995. 00:24:33.309 --> 00:24:36.169 You can see that there is significant injection occurring, 00:24:36.169 --> 00:24:42.480 even as early as the early 2000s, but relatively low levels. 00:24:42.480 --> 00:24:48.220 There was a really large ramp-up of injection in 2012, and injection 00:24:48.220 --> 00:24:53.930 actually peaked in late 2015 at about 15 million cubic meters per month. 00:24:53.930 --> 00:24:55.850 And so to kind of put that in perspective – 00:24:55.850 --> 00:25:00.320 I calculated this the other day – it was about 6,500 Olympic-sized 00:25:00.320 --> 00:25:05.160 swimming pools were month being injected underground. 00:25:05.160 --> 00:25:09.770 And then, since then, the rates of injection have decreased. 00:25:09.770 --> 00:25:11.380 And then, looking at the seismicity, 00:25:11.380 --> 00:25:17.600 really, the seismicity kind of turned on in 2009. 00:25:17.600 --> 00:25:25.350 There was a really big spike of seismicity in 2014. It started to ramp up quite a bit. 00:25:25.350 --> 00:25:29.880 And then, following the reduced injection, 00:25:29.880 --> 00:25:34.000 the seismicity rates have also decreased throughout the area. 00:25:37.040 --> 00:25:42.100 And if you look at different regions of the state, they’ve kind of – 00:25:42.110 --> 00:25:48.360 you can see kind of differences in how injection and seismicity has occurred, 00:25:48.360 --> 00:25:50.700 looking back over time. And so I’ve split it here into kind of the 00:25:50.700 --> 00:25:55.230 western Oklahoma, central Oklahoma, and southern Kansas areas. 00:25:55.230 --> 00:26:03.419 And, again, you can see just this really marked contrast in the – 00:26:03.419 --> 00:26:07.920 in both how injection has occurred and how the sequences have progressed. 00:26:07.920 --> 00:26:10.380 And if we’re thinking about applying physics-based models 00:26:10.380 --> 00:26:14.960 to understand how the sequence has evolved through time, 00:26:14.960 --> 00:26:18.100 this is something we’re going to want to be able to capture – 00:26:18.100 --> 00:26:23.760 this variety of behaviors – as a good test of how well the model is working. 00:26:26.440 --> 00:26:28.020 So getting on to kind of the modeling work 00:26:28.020 --> 00:26:30.559 that I actually performed for this study. 00:26:30.559 --> 00:26:36.330 We’re going to make use of the earthquake nucleation model that was 00:26:36.330 --> 00:26:41.880 really introduced by Jim Dieterich in his paper in JGR back in 1994. 00:26:41.880 --> 00:26:46.710 It’s, like I said, based on rate-and-state friction. 00:26:46.710 --> 00:26:49.539 And the idea is that you have a kind of set of faults 00:26:49.539 --> 00:26:53.520 distributed throughout a volume of rock that is being stressed 00:26:53.520 --> 00:26:58.600 at a tectonic rate – in Oklahoma, likely to be very small. 00:26:58.600 --> 00:27:02.370 And that yields a kind of steady background seismicity rate. 00:27:02.370 --> 00:27:07.580 And in this model, we’re not actually – unlike the first half of the talk, 00:27:07.580 --> 00:27:13.470 we’re not actually modeling individual earthquake rupture along faults. 00:27:13.470 --> 00:27:20.040 And this – kind of this constraint about the set of faults being stressed 00:27:20.040 --> 00:27:22.630 at the tectonic rate to give a background seismicity rate 00:27:22.630 --> 00:27:25.669 actually accounts for a lot of the epistemic uncertainty that would be 00:27:25.669 --> 00:27:30.260 very difficult to capture if we were doing the full physics modeling. 00:27:30.260 --> 00:27:36.779 So, again, we don’t – this kind of assumption accounts for variability 00:27:36.779 --> 00:27:41.710 in the location of the faults and variability in the local state of stress. 00:27:41.710 --> 00:27:45.920 And then, you know, the assumption that there is this kind of set of 00:27:45.920 --> 00:27:48.370 potentially active faults that exists throughout the rock is really 00:27:48.370 --> 00:27:52.539 well-founded based on that work that Martin Schoenball did. 00:27:52.539 --> 00:27:54.740 So then the question is, how do these faults respond 00:27:54.740 --> 00:27:57.980 to arbitrary loading conditions? 00:27:57.980 --> 00:28:04.240 And Paul Segall actually demonstrated, based on Dieterich’s work, that you can 00:28:04.250 --> 00:28:10.860 derive this differential equation for how seismicity rate evolves through time. 00:28:10.860 --> 00:28:15.870 It’s a first-order, nonlinear ordinary differential equation. 00:28:15.870 --> 00:28:19.250 It’s really governed by this parameter, s-dot, 00:28:19.250 --> 00:28:23.429 which is the Coulomb stressing rate affecting the set of faults. 00:28:23.429 --> 00:28:27.910 And so really, the goal of this work is to develop some approximation 00:28:27.910 --> 00:28:31.720 for what these faults are feeling in terms of the state of stress 00:28:31.720 --> 00:28:35.980 and link that to the fluid injection operations. 00:28:40.180 --> 00:28:47.300 So just as a quick aside, for those of you who aren’t super familiar with how the 00:28:47.300 --> 00:28:52.420 Dieterich model works, here’s just a couple example simulations. 00:28:52.429 --> 00:28:56.130 So this is kind of simulating a main shock-aftershock sequence. 00:28:56.130 --> 00:29:00.450 So you subject your set of faults to a static stress change and 00:29:00.450 --> 00:29:04.270 look at how those faults respond. And you see there’s initial jump 00:29:04.270 --> 00:29:09.299 in seismicity, and then a fall-off back to background rates. 00:29:09.299 --> 00:29:14.230 And one of the things I want to point out is that this kind of 00:29:14.230 --> 00:29:18.409 t-to-the-minus-1 fall-off behavior that’s characterized empirically through 00:29:18.409 --> 00:29:25.480 Omori’s law is actually an emergent property of this model. 00:29:25.480 --> 00:29:27.200 So that comes out naturally out of – 00:29:27.200 --> 00:29:33.120 out of that equation, which is really interesting. 00:29:34.080 --> 00:29:35.860 So that was first-step change in stress. 00:29:35.860 --> 00:29:38.120 What about a step change in stressing rate? 00:29:38.120 --> 00:29:42.080 So what I’ve shown here is, now it’s subjecting the set of faults 00:29:42.080 --> 00:29:45.460 to a new stressing rate that’s a factor of 1,000 higher. 00:29:45.460 --> 00:29:51.520 And you’ll see that there is some time where there’s a seismicity rate transient, 00:29:51.520 --> 00:29:56.240 but if you stress that system at a constant rate for long enough, 00:29:56.240 --> 00:29:59.580 the system evolves towards kind of a steady-state seismicity rate 00:29:59.580 --> 00:30:03.940 that depends on the level of the current stressing rate. 00:30:03.940 --> 00:30:05.920 And then the other really interesting thing that I want to 00:30:05.920 --> 00:30:09.730 draw your attention to is it turns out that the time scales 00:30:09.730 --> 00:30:13.860 over which these seismicity rate transients occur is 00:30:13.860 --> 00:30:17.260 inversely proportional to the stressing rate. 00:30:17.260 --> 00:30:21.130 So if you’re stressing it at a very high rate, seismicity rates change quickly. 00:30:21.130 --> 00:30:26.620 And if you’re stressing it at a low rate, it takes longer for the system to respond. 00:30:26.620 --> 00:30:30.240 I think you can actually see this in the data in Oklahoma. 00:30:30.240 --> 00:30:34.500 So wanted to draw your attention to it, and we’ll come back to it later. 00:30:34.500 --> 00:30:38.210 So now getting on to the main task at hand, which is to develop 00:30:38.210 --> 00:30:42.000 a relationship between the fluid – the fluid injection data 00:30:42.000 --> 00:30:45.910 and the stressing conditions along these faults. 00:30:45.910 --> 00:30:49.570 We developed a kind of relatively simple model. 00:30:49.570 --> 00:30:52.970 You know, you might think, at first, to try and develop 00:30:52.970 --> 00:30:58.010 some sort of full 3D fluid-flow numerical model of the system. 00:30:58.010 --> 00:31:01.299 And that might be a good way of going about to do it, but really, 00:31:01.299 --> 00:31:04.850 kind of like I was saying when I started this section of the talk, 00:31:04.850 --> 00:31:08.750 I think that might be intractable in terms of incorporating that level 00:31:08.750 --> 00:31:13.070 of sophistication into a hazard model type of framework. 00:31:13.070 --> 00:31:17.529 So I prefer to keep it simple. Simple is the best option. 00:31:17.529 --> 00:31:22.870 And so kind of what I imagined was, let’s consider the Arbuckle Aquifer 00:31:22.870 --> 00:31:27.200 as a closed tank. Just for – just to start with. 00:31:27.200 --> 00:31:33.160 And the idea is that we know how much mass is going into the system. 00:31:33.169 --> 00:31:36.930 And then the pressure change is just going to be essentially 00:31:36.930 --> 00:31:40.940 governed by the compressibility of the reservoir. 00:31:40.940 --> 00:31:45.200 So that’s really explained well by this equation here. 00:31:45.200 --> 00:31:50.730 It’s just showing the pressurization rate is related linearly to the injection rate. 00:31:50.730 --> 00:31:53.760 And then it’s moderated by the total compressibility of the system, 00:31:53.760 --> 00:31:57.000 which is the product of the bulk volume, the porosity, and the compressibility 00:31:57.000 --> 00:32:00.280 of the core of the rocks and the water. 00:32:01.700 --> 00:32:08.080 So, again, the really cool thing about this is that it’s – is that it’s really simple. 00:32:08.080 --> 00:32:12.360 And it at least captures kind of first-order effects. 00:32:12.360 --> 00:32:17.840 So, but it is a big assumption, and so we’ll see how well it works 00:32:17.840 --> 00:32:20.820 and see if we have to improve the model. 00:32:22.260 --> 00:32:27.799 So I usually don’t show this type of figure in these talks. 00:32:27.799 --> 00:32:34.049 This is just a list of the model parameters and the values that I used in the study. 00:32:34.049 --> 00:32:39.880 But just point out that, you know, each of the model parameters – there’s about 00:32:39.880 --> 00:32:47.540 seven here – are physical. They can all be measured in the field or in the lab. 00:32:47.540 --> 00:32:51.320 And, in fact, we have pretty good estimates for what most of these things 00:32:51.330 --> 00:32:56.530 are in the Arbuckle, for example. 00:32:56.530 --> 00:33:00.830 So hydraulic properties like porosity, compressibility, and reservoir thickness 00:33:00.830 --> 00:33:04.260 are fairly well characterized throughout the area. 00:33:04.260 --> 00:33:07.770 Some of the rate-and-state properties, like the direct effect, can be measured 00:33:07.770 --> 00:33:11.090 in the lab for granite rock, and they really don’t vary that much – 00:33:11.090 --> 00:33:17.799 a factor of 2 or 3 for most types of rock. So first-order estimate is good. 00:33:17.799 --> 00:33:20.270 And probably the least well-constrained parameter 00:33:20.270 --> 00:33:24.590 in the model is the background stressing rate. 00:33:24.590 --> 00:33:29.020 So in Oklahoma, it’s likely to be very small. 00:33:29.020 --> 00:33:33.180 But there have been – you know, at least theoretically – we can’t measure it, 00:33:33.190 --> 00:33:37.950 and there have been some attempts to do so using different types of methods. 00:33:37.950 --> 00:33:44.600 So one last bit of bookkeeping before I show some of the results 00:33:44.600 --> 00:33:47.600 is that we’re going to be comparing these kind of 00:33:47.600 --> 00:33:52.240 physics-based seismicity rate forecasts to the earthquake catalog. 00:33:52.240 --> 00:33:57.480 So we’re using the ComCat catalog. We’ve converted local duration 00:33:57.480 --> 00:34:04.120 magnitudes to a consistent set of moment magnitudes based on the 00:34:04.120 --> 00:34:08.900 CEUS Seismic Source Characterization for Nuclear Facilities report. 00:34:08.900 --> 00:34:13.260 And we’re only considering magnitudes greater than 3. 00:34:14.440 --> 00:34:19.520 And finally, to just mention that this nucleation model 00:34:19.529 --> 00:34:22.799 is not telling us anything about earthquake magnitudes, 00:34:22.799 --> 00:34:28.889 and it’s also not, like I said, considering earthquake rupture propagation along 00:34:28.889 --> 00:34:35.440 any individual fault or the stress transfer and processes associated with that. 00:34:35.440 --> 00:34:41.200 So we’re not going to expect to be able to match aftershock behavior like – 00:34:41.200 --> 00:34:46.119 kind of like you see here in late 2011 for the Prague earthquake. 00:34:46.120 --> 00:34:51.720 So we compared our results against the declustered catalog. 00:34:51.720 --> 00:34:54.620 Tested a different – couple of declustering algorithms. 00:34:54.629 --> 00:34:58.920 So the Gardner and Knopoff algorithm is shown by these red bars here. 00:34:58.920 --> 00:35:03.029 It takes away about 88% of the entire catalog. 00:35:03.029 --> 00:35:06.730 That’s actually the method that is used in the current iterations 00:35:06.730 --> 00:35:10.869 of the USGS one-year hazard model. 00:35:10.869 --> 00:35:16.160 But we are going to opt to use this Reasenberg method – you know, 00:35:16.160 --> 00:35:21.190 a slightly less aggressive declustering method, so to speak. 00:35:21.190 --> 00:35:24.729 And it captures kind of the strong changes that we’ve really seen 00:35:24.729 --> 00:35:28.680 in Oklahoma that have really driven the hazard. 00:35:28.680 --> 00:35:32.080 So, again, we’re going to use this kind of blue earthquake catalog 00:35:32.080 --> 00:35:35.480 that’s shown by the Reasenberg method. 00:35:36.520 --> 00:35:39.959 Okay, like I mentioned, there were a lot of assumptions in our 00:35:39.959 --> 00:35:43.130 simple fluid pressure model. We were going to want to test it 00:35:43.130 --> 00:35:47.600 across a variety of scales and in different parts of the state. 00:35:47.600 --> 00:35:52.860 And so I started just by considering kind of a large – the largest scale. 00:35:52.869 --> 00:35:55.900 So this kind of Oklahoma area of interest, 00:35:55.900 --> 00:35:58.970 where most of the seismicity has occurred. 00:35:58.970 --> 00:36:02.049 So there is a total of 781 wells in that area. 00:36:02.049 --> 00:36:05.239 I’m going to combine all their injection rates 00:36:05.239 --> 00:36:12.450 and plug it into that fluid pressure equation and kind of estimate 00:36:12.450 --> 00:36:16.440 an average pressurization rate over that entire volume. 00:36:17.920 --> 00:36:21.720 So when I do that, these are the values that come out of the calculation. 00:36:21.729 --> 00:36:24.849 So the black line is showing the fluid pressurization, 00:36:24.849 --> 00:36:28.289 or the stressing rate, that’s going to go into the model. 00:36:28.289 --> 00:36:32.279 The black dashed line is the background stressing rate we’re using in this study. 00:36:32.279 --> 00:36:35.039 So you can see there’s kind of two to three orders of magnitude increase 00:36:35.039 --> 00:36:40.950 over time in the stressing rate, and that’s really what’s driving the hazard. 00:36:40.950 --> 00:36:43.450 In terms of the actual values of the pressure change that 00:36:43.450 --> 00:36:46.890 we’re calculating in the Arbuckle – so that’s shown by this blue line here. 00:36:46.890 --> 00:36:54.130 Again, it’s just a simple tank model. And so, with continued injection, 00:36:54.130 --> 00:36:56.849 no flow can go out, and pressure can only rise. 00:36:56.849 --> 00:36:58.710 We’re kind of calculating pressure changes 00:36:58.710 --> 00:37:03.079 on the order of about 3 megapascals over that 20 years. 00:37:03.079 --> 00:37:09.579 So there’s not a lot of pressure data to constrain that and test whether that’s 00:37:09.579 --> 00:37:14.349 truly accurate or not, but at least it’s a reasonable value. It’s not too big. 00:37:14.349 --> 00:37:18.600 So that’s a good gut check. 00:37:21.280 --> 00:37:24.540 Okay, so we’re going to use that stressing function and 00:37:24.549 --> 00:37:28.239 plug it into this differential equation. We’re going to solve it numerically 00:37:28.239 --> 00:37:35.279 using a Runge-Kutta method to find the solution for the seismicity rate over time. 00:37:35.279 --> 00:37:39.709 And when I do that, the red line that I’m showing there 00:37:39.709 --> 00:37:46.160 is the seismicity rate forecast that comes out of that solution. 00:37:46.160 --> 00:37:51.140 And what I was really happy about seeing what that this captures the 00:37:51.150 --> 00:37:57.380 onset of seismicity there in late 2009. It captures the kind of upswing in 00:37:57.380 --> 00:38:02.869 seismicity and the timing of the peak seismicity rate in late 2015. 00:38:02.869 --> 00:38:07.890 And then it also captures a reduction in seismicity with – 00:38:07.890 --> 00:38:11.940 you know, in response to the reduced saltwater disposal. 00:38:11.940 --> 00:38:16.260 And, again, just making the point that this doesn’t depend on the earthquake 00:38:16.260 --> 00:38:19.299 catalog at all. There hasn’t been any calibration against the catalog. 00:38:19.299 --> 00:38:23.339 It’s just simply putting in reasonable physical properties 00:38:23.339 --> 00:38:27.039 for the Arbuckle rock and the rate-and-state friction properties 00:38:27.040 --> 00:38:32.380 and plugging it into this differential equation and solving that through time. 00:38:33.360 --> 00:38:36.800 And, again, you’ll see that it matches well with the declustered catalog 00:38:36.800 --> 00:38:41.560 and not the – not the full catalog that has the aftershocks included. 00:38:43.809 --> 00:38:48.559 So let’s see how it works at a slightly smaller scale in these different parts of 00:38:48.559 --> 00:38:52.979 the state where I showed that there was that marked difference in behavior. 00:38:52.979 --> 00:38:56.819 So now I’m, again, separating it into the central Oklahoma, 00:38:56.820 --> 00:39:01.880 western Oklahoma, and southern Kansas study areas. 00:39:01.880 --> 00:39:06.760 And now only considering wells within each region, scaling the bulk volume. 00:39:06.760 --> 00:39:10.630 Otherwise, I haven’t introduced any heterogeneity into the problem. 00:39:10.630 --> 00:39:13.920 Just keeping all the model properties the same. 00:39:13.920 --> 00:39:16.680 And, again, kind of calculating 00:39:16.680 --> 00:39:20.940 average pressurization rate in each of these regions. 00:39:23.060 --> 00:39:28.140 So then, when I solve that differential equation separately for the three areas, 00:39:28.140 --> 00:39:31.940 again, the red line is showing the kind of seismicity rate forecast. 00:39:31.940 --> 00:39:34.520 And what you’ll see is, again, it’s capturing that 00:39:34.520 --> 00:39:37.319 onset of seismicity really quite well. 00:39:37.319 --> 00:39:40.200 In the central Oklahoma area, you can see we’re kind of 00:39:40.200 --> 00:39:45.150 overestimating a bit in that period between 2010 and 2014. 00:39:45.150 --> 00:39:47.980 Perhaps underestimating the peak that – you know, 00:39:47.980 --> 00:39:52.910 the magnitude of the seismicity rate at the peak a little bit. 00:39:52.910 --> 00:39:56.319 But overall, it’s capturing the right order of magnitude and the timings 00:39:56.319 --> 00:40:00.900 of these sequences quite well. And you’ll see that it’s capturing the 00:40:00.900 --> 00:40:07.160 different regions just kind of as an emergent property of the model. 00:40:07.160 --> 00:40:10.400 You see it captures it really well in western Oklahoma. 00:40:10.400 --> 00:40:14.730 In southern Kansas, the scale on the injection rate and the 00:40:14.730 --> 00:40:18.199 seismicity rate is quite a bit smaller. So there’s just not as much injection 00:40:18.199 --> 00:40:21.969 going on, and also just not as much seismicity going on in that region. 00:40:21.969 --> 00:40:27.000 But, again, overall, capturing the right trends that we would like to see. 00:40:30.660 --> 00:40:36.220 So now we’ll look at a slightly smaller scale. 00:40:36.220 --> 00:40:41.839 I’m going to just look at these four regions – Oklahoma City, Cherokee, 00:40:41.839 --> 00:40:46.960 Harper, and Cushing, signified by these red circles here. 00:40:46.960 --> 00:40:50.119 And, again, only consider wells within those volumes – 00:40:50.120 --> 00:40:55.840 just scale the bulk volume and keep everything else the same. 00:40:55.840 --> 00:40:59.479 And what’s really cool about looking at these different areas is that, 00:40:59.479 --> 00:41:02.430 again, you can see that really strong contrast in behavior. 00:41:02.430 --> 00:41:07.059 So for the Cherokee area, for example, it’s kind of more typical of what we saw 00:41:07.059 --> 00:41:12.039 in western Oklahoma, in general. But contrast that to what’s going on 00:41:12.039 --> 00:41:16.269 in Oklahoma City, where the injection rate has been quite a bit 00:41:16.269 --> 00:41:20.750 less than in Cherokee area, but it’s been relatively constant 00:41:20.750 --> 00:41:26.961 since back about 2006, 2007. And what – you know, 00:41:26.961 --> 00:41:30.069 in terms of the observed seismicity, you kind of saw a trickle of seismicity 00:41:30.069 --> 00:41:36.059 that whole period between zero and five earthquakes per month. 00:41:36.059 --> 00:41:40.699 And the red line is kind of showing a similar trend in terms of that behavior. 00:41:40.699 --> 00:41:46.190 I think one area where the model is not working so well is the Cushing region. 00:41:46.190 --> 00:41:50.880 And, you know, what that tells me is that maybe we want to go in and 00:41:50.880 --> 00:41:56.130 look at that area a little bit more closely. Perhaps try and model the hydrologic 00:41:56.130 --> 00:42:01.229 system there more appropriately and see if we can do a better job at fixing that. 00:42:01.229 --> 00:42:06.500 But nonetheless, this gives us a good sense of where to go look for where the 00:42:06.500 --> 00:42:13.080 model is not behaving as close to what we’ve seen in the observed sequence. 00:42:14.190 --> 00:42:16.200 So it looked good, and then we wanted to test 00:42:16.219 --> 00:42:19.960 the accuracy of the system – or, of the model. 00:42:19.960 --> 00:42:23.749 And so we wanted to test it against the next-best model, which are kind of 00:42:23.749 --> 00:42:28.140 these observational models that are used in the one-year hazard forecasts. 00:42:28.140 --> 00:42:31.779 So these types of seismicity rate forecasts 00:42:31.780 --> 00:42:37.180 just use last year’s earthquake rate to predict this year’s rate. 00:42:37.180 --> 00:42:41.400 And we wanted to test, or evaluate, the accuracy over the whole period 00:42:41.410 --> 00:42:46.699 of the main sequence between 2008 and 2016, so we developed 00:42:46.699 --> 00:42:51.900 a set of kind of nine pseudo-perspective one-year forecasts. 00:42:51.900 --> 00:42:53.579 And then we compared our physics-based, 00:42:53.579 --> 00:42:57.380 observational-based models against the observed rates 00:42:57.380 --> 00:43:01.670 and quantified the accuracy using two different approaches. 00:43:01.670 --> 00:43:04.499 One is the – kind of a likelihood test that’s been 00:43:04.499 --> 00:43:07.059 recommended by SCEC and the CSEP community. 00:43:07.060 --> 00:43:12.280 And we also just looked at root mean square error kind of analysis. 00:43:13.240 --> 00:43:18.580 And so these two plots are showing the results of the accuracy quantification. 00:43:21.829 --> 00:43:24.080 What I’m showing on the X axis is actually 00:43:24.089 --> 00:43:26.640 kind of a forecast update frequency. 00:43:26.640 --> 00:43:32.249 And so the data point on the far left, in terms of the squares, 00:43:32.249 --> 00:43:39.360 is essentially the one-year hazard model. So that’s just using last year’s 00:43:39.360 --> 00:43:42.820 earthquakes to forecast this year’s seismicity. 00:43:42.820 --> 00:43:47.240 And what you’ll see is that this blue line is our physics-based model, 00:43:47.240 --> 00:43:52.829 and it’s doing quite a bit better than these observational approaches. 00:43:52.829 --> 00:43:57.210 Just for reference, the red line is a set of random models that were 00:43:57.210 --> 00:44:02.519 just generated by drawing seismicity rates from the set of observed rates. 00:44:02.519 --> 00:44:05.140 And so both the models are doing better than the random model, 00:44:05.140 --> 00:44:08.079 so that’s a good thing. 00:44:08.079 --> 00:44:12.539 But kind of following up on some of the work that Andrea has 00:44:12.539 --> 00:44:15.930 presented at AGU in the past and also some of the work that Morgan Moschetti 00:44:15.930 --> 00:44:21.990 is doing over in Golden, they found that, if you update your hazard model 00:44:21.990 --> 00:44:25.440 more frequently, you can actually achieve higher accuracy. 00:44:25.440 --> 00:44:28.729 And so that’s what I’m showing on the X axis. 00:44:28.729 --> 00:44:32.799 So if you use last month’s seismicity to forecast this month’s seismicity, 00:44:32.799 --> 00:44:34.920 you can obviously do a much more accurate job. 00:44:34.920 --> 00:44:37.449 And so you see that there is some kind of tradeoff 00:44:37.449 --> 00:44:43.690 where the observational models really do offer a lot of good information, 00:44:43.690 --> 00:44:48.119 and there might be some kind of optimal tradeoff where you can kind of combine 00:44:48.119 --> 00:44:50.969 both the physics-based approach and the observational approach. 00:44:50.969 --> 00:44:54.849 That’s something to look into a bit in the future. 00:44:54.849 --> 00:44:59.339 So I would be remiss if I didn’t acknowledge the only other model 00:44:59.339 --> 00:45:04.819 out there to – that’s using injection data to forecast seismicity rates. 00:45:04.819 --> 00:45:07.869 This is a model that was introduced by Cornelius Langenbruch, 00:45:07.869 --> 00:45:10.700 who is sitting here in the audience today, 00:45:10.700 --> 00:45:15.880 and Mark Zoback in a paper that appeared last year in Science Advances. 00:45:15.880 --> 00:45:18.489 It’s an empirical model. It kind of operates on the 00:45:18.489 --> 00:45:21.579 assumption that there is a Gutenberg-Richter relationship 00:45:21.579 --> 00:45:24.359 for the sequence of earthquakes. And essentially, you have a 00:45:24.359 --> 00:45:28.999 time-dependent “a” value that’s some function of the injection data. 00:45:28.999 --> 00:45:33.140 So the way the theory goes is that the cumulative earthquake count is related to 00:45:33.140 --> 00:45:37.410 the cumulative injection volume. And then they introduced this kind of 00:45:37.410 --> 00:45:42.140 seismogenic index parameter that’s calculated empirically through a 00:45:42.140 --> 00:45:46.660 regression with the injection data and the earthquake data. 00:45:47.560 --> 00:45:52.760 So, like I said – so if the cumulative earthquake count 00:45:52.769 --> 00:45:55.069 scales with the cumulative volume injected, 00:45:55.069 --> 00:45:57.999 then the seismicity rate scales with the injection rate. 00:45:57.999 --> 00:46:02.479 So actually, that’s somewhat similar to what I’m actually doing in my model. 00:46:02.480 --> 00:46:07.160 So that’s a – like I said, an interesting similarity there. 00:46:09.349 --> 00:46:15.100 If you actually look at the data, you’ll see that it doesn’t actually 00:46:15.109 --> 00:46:18.990 follow this kind of one-to-one trend. But what Cornelius found was that, 00:46:18.990 --> 00:46:23.569 if you do introduce some additional flexibility in the model, 00:46:23.569 --> 00:46:28.269 and allow for a time lag that accounts for changes – kind of a time lag 00:46:28.269 --> 00:46:32.509 between when changes in injection and changes in seismicity occur, 00:46:32.509 --> 00:46:35.499 and you also account for an injection threshold below which 00:46:35.499 --> 00:46:38.470 no seismicity is expected to occur, then you can actually do 00:46:38.470 --> 00:46:42.940 a really good job at matching the sequence. 00:46:42.940 --> 00:46:47.460 Again, it requires a sufficient amount of earthquake observations 00:46:47.460 --> 00:46:54.910 to calibrate this model against. And Cornelius found that it was – 00:46:54.910 --> 00:46:57.620 it was necessary to kind of perform that calibration separately in 00:46:57.620 --> 00:47:02.040 each of the zones, which increases 00:47:02.049 --> 00:47:06.450 the model complexity a bit. And then they also introduced, finally, 00:47:06.450 --> 00:47:11.489 an empirical Omori-decay law to capture the reduction in seismicity 00:47:11.489 --> 00:47:15.120 once it fell below the injection threshold. 00:47:15.980 --> 00:47:23.920 So, like I said, Cornelius found these various injection thresholds in 00:47:23.920 --> 00:47:27.499 different parts of the state and then calibrated that model through 00:47:27.499 --> 00:47:34.009 different time periods and found that, if you calibrated the model through 2014, 00:47:34.009 --> 00:47:37.289 and that the model was able to predict the rest of the sequence well. 00:47:37.289 --> 00:47:42.099 So they found that the model could have been used as early as 2014, 00:47:42.100 --> 00:47:45.020 which is – which is definitely very useful. 00:47:46.229 --> 00:47:49.600 So I also want to compare against this type of a model. 00:47:49.609 --> 00:47:52.619 So I’ve done some preliminary work looking at that. 00:47:52.619 --> 00:47:57.729 I used a slightly different injection and earthquake catalog than was 00:47:57.729 --> 00:48:00.859 presented in that study, so I couldn’t use the exact values that they cited 00:48:00.859 --> 00:48:07.619 for the seismogenic index parameters. But I calibrated it using the – 00:48:07.619 --> 00:48:13.299 kind of the procedure laid out in that study and calibrated it against 00:48:13.299 --> 00:48:16.700 my own data set and was able to kind of come up with similar, 00:48:16.700 --> 00:48:19.719 more or less, seismogenic index parameters, which is good – 00:48:19.719 --> 00:48:25.420 again, a good gut check that I’m following the right approach. 00:48:25.420 --> 00:48:29.660 And you can see this is kind of the comparison between my model, 00:48:29.660 --> 00:48:32.539 which is the red line, and the black dashed line, 00:48:32.540 --> 00:48:36.760 which is the seismogenic index forecast that I’ve created. 00:48:36.760 --> 00:48:40.380 And they both have their areas where they’re doing better than 00:48:40.380 --> 00:48:45.840 the other or, you know, both have their pros and cons. 00:48:47.840 --> 00:48:53.459 When you actually look at the accuracy, they’re performing relatively similarly. 00:48:53.459 --> 00:48:56.229 So at least for these parameters, the physics-based is slightly – 00:48:56.229 --> 00:49:00.800 doing slightly better than that seismogenic index model. 00:49:01.860 --> 00:49:06.459 It might require more work to calibrate that model slightly and, 00:49:06.459 --> 00:49:08.390 you know, formulate the calibration a little differently. 00:49:08.390 --> 00:49:10.839 And we’ll see how that will affect things. 00:49:10.839 --> 00:49:14.369 But just point out that both of the models are doing better 00:49:14.369 --> 00:49:16.910 than the observational approach. And what that suggests to me 00:49:16.910 --> 00:49:21.359 is that you’re really able to capture a lot of useful information 00:49:21.359 --> 00:49:24.940 using the injection data. 00:49:24.940 --> 00:49:27.920 So it’s worthwhile to think about how to actually incorporate 00:49:27.920 --> 00:49:32.700 these operational data into kind of a seismic hazard forecast. 00:49:34.800 --> 00:49:40.650 So just a few implications for the – from the study for – in terms of 00:49:40.650 --> 00:49:43.529 managing the hazard associated with induced seismicity. 00:49:43.529 --> 00:49:51.170 Really, the benefit of using a physics- based model is it provides that kind of 00:49:51.170 --> 00:49:59.319 physics-based intuition, or insight, into the processes that are actually happening. 00:49:59.319 --> 00:50:02.210 First is that, like I said, using the injection rates and the well locations 00:50:02.210 --> 00:50:10.150 did seem to improve the forecasts. So that’s a useful outcome of this study. 00:50:10.150 --> 00:50:15.160 We found that the seismicity rate is governed by the stressing rate. 00:50:15.160 --> 00:50:20.039 You know, a lot of people ask, what’s – the question of, what kind of 00:50:20.039 --> 00:50:23.989 pressure change in terms of MPa is required to cause slip on a fault. 00:50:23.989 --> 00:50:27.619 And that’s probably the right way of looking at it for any individual 00:50:27.619 --> 00:50:31.310 sequence or along any individual fault. But when you’re looking at a 00:50:31.310 --> 00:50:35.249 broader scale, we found that the pressurization rate was really 00:50:35.249 --> 00:50:38.769 a lot more important and the thing that’s driving behavior. 00:50:38.769 --> 00:50:44.670 And it seems to be characterized appropriately through this kind of 00:50:44.670 --> 00:50:47.559 compressibility model, at least to first order. 00:50:47.559 --> 00:50:50.039 It also – you know, we’re only considering pressure within 00:50:50.039 --> 00:50:54.749 the Arbuckle, and that really highlights the fact that – the fact that 00:50:54.749 --> 00:50:59.920 this basal aquifer overlying basement rock is really a key driver 00:50:59.920 --> 00:51:02.229 in terms of the hazard. So it might be worth thinking about 00:51:02.229 --> 00:51:05.520 injecting in shallower formations, for example. 00:51:06.580 --> 00:51:11.890 The – so the differential equation suggests that there is a 00:51:11.890 --> 00:51:17.519 steady-state seismicity rate that’s equal to – or essentially 00:51:17.520 --> 00:51:22.480 related to ratio of the stressing rate to the background stressing rate. 00:51:22.480 --> 00:51:26.299 And this is one of the big reasons why this model actually works. 00:51:26.299 --> 00:51:31.079 So, for example, in 2015, we kind of calculated a stressing rate 00:51:31.079 --> 00:51:36.049 increase by a factor of 900, and we observed about 865 earthquakes 00:51:36.049 --> 00:51:41.249 in that year, which is about a 900-fold increase. 00:51:41.249 --> 00:51:47.450 What’s happening this year so far, in 2017, the model is kind of calculating 00:51:47.450 --> 00:51:53.430 a stressing rate of about 300 – an increase of about a factor of 300. 00:51:53.430 --> 00:51:56.779 And in the six months from January to June, we observed 00:51:56.779 --> 00:52:02.009 about 150 earthquakes – a 300-times increase. 00:52:02.009 --> 00:52:09.969 And so what that really suggests is that there might be some seismicity – 00:52:09.969 --> 00:52:14.761 some injection rate that will result in a seismicity rate 00:52:14.761 --> 00:52:19.670 that remains below some threshold that might be tolerable. 00:52:19.670 --> 00:52:23.109 So that’s interesting from an operational standpoint. 00:52:23.109 --> 00:52:29.410 And then, finally, just looking at how – I guess another property of the model 00:52:29.410 --> 00:52:33.939 is that the seismicity rate transients – the time scale for changes in seismicity 00:52:33.939 --> 00:52:38.160 to occur – scales like 1 over the stressing rate. 00:52:38.160 --> 00:52:44.749 And so I think you can really see this in the data. 00:52:44.749 --> 00:52:50.430 So if you look at the period between, you know, 1995 and the early 2000s, 00:52:50.430 --> 00:52:54.130 you’re stressing the system at a relatively low rate. 00:52:54.130 --> 00:52:58.839 And if you plug in the numbers, it turns out to be, you know, 00:52:58.839 --> 00:53:03.529 this kind of characteristic time scale is on the order of about five to 10 years. 00:53:03.529 --> 00:53:06.559 So it really comes to no surprise that it took so long for the sequence 00:53:06.559 --> 00:53:11.819 to start up in late 2009. But if you look again at the highest 00:53:11.819 --> 00:53:17.539 stressing rate period in late 2015, again, if you plug in the numbers you kind of 00:53:17.539 --> 00:53:22.869 get a characteristic time scale of about six months or so, 00:53:22.869 --> 00:53:27.859 which again, is similar to the lag that we saw between when 00:53:27.860 --> 00:53:33.300 injection dropped and then seismicity began to drop as well. 00:53:33.300 --> 00:53:39.559 So, again, really interesting insight that the – you know, the underlying 00:53:39.559 --> 00:53:44.199 equations can give to understanding how this sequence evolved. 00:53:44.199 --> 00:53:48.130 All right, so this is the last slide. Just want to highlight 00:53:48.130 --> 00:53:51.180 some of the current work that we’re doing here at the USGS. 00:53:51.180 --> 00:53:55.319 This is being led by Andy Barbour. 00:53:55.319 --> 00:53:59.619 So we’re essentially trying to start monitoring some of the 00:53:59.619 --> 00:54:04.229 injection wells throughout the Arbuckle. So currently, we have one well 00:54:04.229 --> 00:54:06.269 logged up with a pressure logging device. 00:54:06.269 --> 00:54:08.440 So we have a co-located seismometer there. 00:54:08.440 --> 00:54:11.049 We’re preparing to try and expand that network 00:54:11.049 --> 00:54:14.630 to three or four more wells in the near future. 00:54:14.630 --> 00:54:17.769 And we’re also discussing the possibility of perhaps performing 00:54:17.769 --> 00:54:22.329 some actual injection tests that – so that we can try and back out 00:54:22.329 --> 00:54:29.890 some of the properties of the Arbuckle that would inform my model. 00:54:29.890 --> 00:54:33.309 So thanks a lot for your attention. I really appreciate it. 00:54:33.309 --> 00:54:37.520 [ Applause ] 00:54:37.520 --> 00:54:41.880 - Thank you, Jack, for a really interesting talk. Do we have questions for Jack? 00:54:46.840 --> 00:54:49.180 - Thanks, Jack. - Yeah. 00:54:49.180 --> 00:54:54.680 - That was really impressive for such a short time here at the Survey so far. 00:54:54.680 --> 00:54:59.960 So one of the sort of conceptual – you used a model of a – 00:54:59.960 --> 00:55:05.799 of a tank that has no leaks in it. So kind of the – one question 00:55:05.799 --> 00:55:08.859 that comes up, and maybe it comes up in comparison with Cornelius’ model, 00:55:08.859 --> 00:55:15.249 is are there leaks in the tank? You know, and how does that – 00:55:15.249 --> 00:55:19.339 you know, that may be something that affects the long-term predictions. 00:55:19.339 --> 00:55:23.200 And also, you know, something that may vary – or, certainly does 00:55:23.200 --> 00:55:25.520 vary over space and time. - Yeah. 00:55:25.520 --> 00:55:30.390 - Do you have any ideas on that idea? - Definitely. Yeah, that’s – I mean, 00:55:30.390 --> 00:55:32.740 that is really a major assumption of the model. 00:55:32.740 --> 00:55:36.279 I mean, the Arbuckle itself is laterally extensive. 00:55:36.279 --> 00:55:40.920 So, in theory, flow should be able to, you know, flow out laterally, at least, 00:55:40.920 --> 00:55:44.619 and then also, potentially, down into the basement, which would 00:55:44.619 --> 00:55:48.710 all appear as kind of a leaky tank. - Yeah. 00:55:48.710 --> 00:55:55.959 - One thing we’ve looked into – so actually that closed-system model 00:55:55.959 --> 00:55:59.789 is an end member case, in fact. And so we wanted to try and just 00:55:59.789 --> 00:56:05.739 test how accurate that type of model is. And so, actually Paul Hsieh helped me 00:56:05.739 --> 00:56:11.140 come up with this little analysis. So the idea was, let’s consider constant 00:56:11.140 --> 00:56:16.880 injection over, you know, a large region – a radius 100 kilometers here. 00:56:16.880 --> 00:56:20.640 And, essentially, since those wells are all distributed, we can maybe approximate 00:56:20.640 --> 00:56:27.549 that as kind of constant rate injection. We looked at how – you know, 00:56:27.549 --> 00:56:33.619 we considered a 10-year time span and looked at how pressure in the system 00:56:33.619 --> 00:56:37.019 evolved and compared it to the closed system versus the infinite domain. 00:56:37.019 --> 00:56:44.249 And so, you know, that then becomes affected by permeability or diffusivity. 00:56:44.249 --> 00:56:49.519 So what it turns out – so, you know, Arbuckle permeability is probably 00:56:49.519 --> 00:56:53.479 thought to range between 100 millidarcys and 1,000 millidarcys – 00:56:53.479 --> 00:56:55.690 so kind of these two frames on the right. 00:56:55.690 --> 00:56:59.420 And what I’m plotting here is just the radial distance away, 00:56:59.420 --> 00:57:01.779 you know, from that domain. 00:57:01.779 --> 00:57:06.989 And on the Y axis is the relative pressure, so the pressure in the 00:57:06.989 --> 00:57:09.630 open system compared to the pressure in the closed system. 00:57:09.630 --> 00:57:13.260 And what you’ll see is that, after 10 years – I’m showing the kind of 00:57:13.260 --> 00:57:16.460 average pressure as this dashed black line – and even for 00:57:16.469 --> 00:57:20.700 the highest permeability case, you’ve got a – you know, 00:57:20.700 --> 00:57:24.779 only a factor of – well, I guess about a 35% reduction. 00:57:24.779 --> 00:57:28.880 So, you know, that’s also kind of on the variability of the hydraulic 00:57:28.880 --> 00:57:33.560 properties in terms of porosity, compressibility, and things like that. 00:57:34.700 --> 00:57:37.479 You know, our model doesn’t capture anything that’s happening 00:57:37.479 --> 00:57:39.569 on the right-hand side of these plots where there’s actually 00:57:39.569 --> 00:57:42.009 some diffusion outside of the system. 00:57:42.009 --> 00:57:46.989 And so that might help if we want to consider that in a later study. 00:57:46.989 --> 00:57:50.279 But at least to first order, we’re not doing too bad with this model. 00:57:50.279 --> 00:57:55.900 And it’s really because injection is occurring over such a broad area. 00:57:55.900 --> 00:58:01.539 And then other point just to think about is, I’m only considering 00:58:01.539 --> 00:58:03.999 fluid pressure in the Arbuckle. But what we’re – you know, 00:58:03.999 --> 00:58:07.400 all the seismicity is occurring 1 to 4 kilometers below that. 00:58:07.400 --> 00:58:13.119 And so what that tells me is that these faults are behaving relatively – 00:58:13.119 --> 00:58:15.640 you know, they have relatively high permeability so that the 00:58:15.640 --> 00:58:19.150 time scales for pressure diffusion down along those faults 00:58:19.150 --> 00:58:24.880 is less than the time scales related to the seismicity rate changes. 00:58:24.880 --> 00:58:27.809 So I think that’s a useful insight that comes out of this study. 00:58:27.809 --> 00:58:31.280 - Yeah. It seems like the pressure data that you guys are planning to collect 00:58:31.280 --> 00:58:33.820 would be pretty valuable in helping to – answering … 00:58:33.820 --> 00:58:36.480 - I hope so. Yeah. - Thanks. 00:58:39.460 --> 00:58:44.189 - I was just wondering, when you were doing your forecast, what do 00:58:44.189 --> 00:58:48.940 you assume about the injection rate? Or does it not matter so much 00:58:48.940 --> 00:58:52.120 because everything – the forecast is primarily 00:58:52.130 --> 00:58:56.080 based on what has – what has – the injection rate in the past? 00:58:56.500 --> 00:59:00.459 Like, if you’re doing a forecast one year – for one year, 00:59:00.459 --> 00:59:03.469 what do you assume for the injection rate for that year? 00:59:03.469 --> 00:59:08.549 - I see. So you’re – especially that comes into the accuracy of the model because, 00:59:08.549 --> 00:59:12.660 in theory, you would never have the injection data for the forecast period. 00:59:12.660 --> 00:59:16.499 Yeah, we did some testing with that. And what we found was that, 00:59:16.499 --> 00:59:21.130 if we just assume that, for the following year, 00:59:21.130 --> 00:59:24.971 the injection rate was equal to the average of the last six months or so – 00:59:24.971 --> 00:59:29.949 or, the injection data, that it really didn’t affect that forecast calculation 00:59:29.949 --> 00:59:34.670 very significantly at all. So that suggests that, yeah, it could – 00:59:34.670 --> 00:59:39.219 this technique could be used in kind of a practical situation by just assuming that 00:59:39.220 --> 00:59:43.719 injection rate is some average of the most recent rates. 00:59:45.760 --> 00:59:49.339 - Thanks, Jack. Really clear, interesting talk. 00:59:49.339 --> 00:59:53.559 A couple questions. One – well, part A and B are pretty naive. 00:59:53.559 --> 00:59:58.059 Sorry I walked in a little late. What was the reason for the injection 00:59:58.060 --> 01:00:02.680 fall-off around 2014? I have guesses, but I’m not sure if you said that. 01:00:02.680 --> 01:00:06.600 - Yeah. So, you know, the reason why the injection rates have decreased 01:00:06.609 --> 01:00:11.759 so rapidly in the last couple years is primarily due to external factors. 01:00:11.759 --> 01:00:13.719 So the price of oil has fallen, 01:00:13.719 --> 01:00:16.410 and production has decreased a lot in that area. 01:00:16.410 --> 01:00:20.170 And so with decreased production, there’s just less water to dispose of. 01:00:20.170 --> 01:00:22.950 Now, there’s also – the Oklahoma Corporation Commission 01:00:22.950 --> 01:00:30.719 also implemented a mandate to keep injection below a certain threshold. 01:00:30.719 --> 01:00:34.699 And injection rates are – you know, were already below that when the 01:00:34.699 --> 01:00:38.359 mandate went into effect, mostly due to those external forces. 01:00:38.359 --> 01:00:41.430 - So a combination of sort of market and regulatory forces. 01:00:41.430 --> 01:00:44.410 - Yeah. - And do – could you go to the slide 01:00:44.410 --> 01:00:50.420 where you showed the different regions – the Oklahoma through Kansas? 01:00:51.800 --> 01:00:57.740 - Maybe like this one? - Yeah. So – or, well, no, I think the 01:00:57.740 --> 01:01:02.500 one with Kansas, in particular, because the difference in Kansas – yeah. 01:01:02.500 --> 01:01:06.239 The difference in Kansas in the model performance compared to the others, 01:01:06.239 --> 01:01:09.579 is that significant? You mentioned it a little bit, but is that – 01:01:09.579 --> 01:01:14.829 is that telling us something physical? Or is it just that the rate is so much 01:01:14.829 --> 01:01:19.989 less in Kansas? But then why is the rate so much less in Kansas? 01:01:19.989 --> 01:01:24.869 - Yeah. Well, you know, what comes out of this model is that the rate is so much 01:01:24.869 --> 01:01:30.329 less in Kansas because the injection rates in that area are much lower. 01:01:30.329 --> 01:01:35.860 So this is a scale of 8 million barrels – or, cubic meters per month versus 2 million. 01:01:35.860 --> 01:01:40.489 So it’s quite a bit less. Of course, it is a bit smaller area as well. 01:01:40.489 --> 01:01:44.529 So that’s taken into account in the seismicity rate forecasts. 01:01:44.529 --> 01:01:50.460 In terms of how well the red line compares to the blue bars, I kind of 01:01:50.460 --> 01:01:55.019 chalk that up to kind of stochastic variability of earthquake processes. 01:01:55.019 --> 01:01:57.259 So whether you’re talking about five earthquakes 01:01:57.259 --> 01:02:01.509 per month or two or three, you’re kind of similar there. 01:02:01.509 --> 01:02:05.260 But it would be worthwhile maybe quantifying that a little bit closer. 01:02:05.260 --> 01:02:13.360 So far, those accuracy quantifications were for the whole statewide scale, so … 01:02:15.780 --> 01:02:22.400 [ Silence ] 01:02:23.360 --> 01:02:25.239 - Thanks, Jack. That’s really interesting. 01:02:25.239 --> 01:02:30.829 I just was curious. You stated that you’re not 01:02:30.829 --> 01:02:34.179 considering the magnitude of – or, the distribution of earthquake magnitude, 01:02:34.179 --> 01:02:38.739 but what is the data on that as far as, if you lower injection rate, 01:02:38.739 --> 01:02:44.069 is there any hope that that limits the large-magnitude earthquakes? 01:02:44.069 --> 01:02:47.420 - I guess observationally, that would be a big no. 01:02:47.420 --> 01:02:53.600 Because injection rates fell down in late 2015 and continued through 2016. 01:02:53.609 --> 01:02:57.829 In fact, 2016 was the year where we had the most large earthquakes. 01:02:57.829 --> 01:03:02.180 And we had three magnitude 5s and larger. 01:03:03.080 --> 01:03:05.860 I guess, empirically, what that means is that the b value – 01:03:05.860 --> 01:03:10.360 the Gutenberg-Richter b value is changing slope over time. 01:03:10.360 --> 01:03:16.619 And if you had reason to – you know, that that was occurring, 01:03:16.619 --> 01:03:19.459 you could use that information within the framework of this model 01:03:19.459 --> 01:03:23.400 to kind of forecast the probabilities of larger earthquakes. 01:03:23.400 --> 01:03:29.809 But in terms of a physical explanation for why that might be occurring, 01:03:29.809 --> 01:03:33.010 this particular model wouldn’t really get you there. 01:03:33.010 --> 01:03:37.799 But I think that, if – you know, the way to get there would be to 01:03:37.800 --> 01:03:43.620 first use the work that Martin Schoenball and Bill Ellsworth did by – 01:03:43.620 --> 01:03:47.119 you know, where they were able to show you what the faults 01:03:47.119 --> 01:03:48.630 really look like in this part of the country. 01:03:48.630 --> 01:03:52.049 So you can kind of characterize that statistically, perhaps. 01:03:52.049 --> 01:03:55.990 And then, using, like, CFRAC, the model that I was showing 01:03:55.990 --> 01:03:59.099 in the first half of the talk, which is really forecasting 01:03:59.099 --> 01:04:04.419 earthquake rupture along individual faults, which can get you to magnitudes, 01:04:04.420 --> 01:04:09.380 coupling in the injection data and how that’s changed over time, 01:04:09.380 --> 01:04:15.420 and populating those models with stochastic fault distributions based on 01:04:15.429 --> 01:04:19.880 the data that Martin showed, you might be able to get there with that. 01:04:19.880 --> 01:04:23.919 But I don’t have a good handle for why that might be. 01:04:23.920 --> 01:04:29.000 It could just be that – it could be also just a sampling problem. 01:04:29.000 --> 01:04:33.420 So, you know, I don’t think it’s totally unlikely 01:04:33.429 --> 01:04:35.789 given the large number of earthquakes that we’ve seen, 01:04:35.789 --> 01:04:39.900 that we had three magnitude 5s and larger this last year. 01:04:39.900 --> 01:04:45.560 So it could be physics-based, or it could just be a statistics thing. 01:04:47.620 --> 01:04:55.320 [ Silence ] 01:04:56.360 --> 01:04:59.260 - Thank you. It was a really impressive study. 01:04:59.260 --> 01:05:08.410 It was conceptually clear and simple in that way. And I really enjoyed your talk. 01:05:08.410 --> 01:05:14.460 You mentioned that the parameters that you used in your model were reasonable. 01:05:14.460 --> 01:05:20.829 And my question is, how much did you tune those parameters 01:05:20.829 --> 01:05:24.779 to the particular area that you were studying? 01:05:24.779 --> 01:05:30.619 And did you do any sensitivity studies as to, if you varied these reasonable 01:05:30.620 --> 01:05:36.660 parameters within reasonable limits, what – how that affected your results. 01:05:36.660 --> 01:05:39.599 - Yeah. Yeah, thanks for that question. 01:05:39.599 --> 01:05:45.929 So really, the key fundamental takeaways for why the model works, 01:05:45.929 --> 01:05:48.900 there’s kind of a – I kind of made those two points at the end. 01:05:48.900 --> 01:05:53.869 So you want to get the steady-state seismicity rate correct. 01:05:53.869 --> 01:05:57.380 So that depends on the stressing rate to the background stressing rate. 01:05:57.380 --> 01:05:59.480 So you have that to play with. 01:05:59.480 --> 01:06:02.419 And the second thing is, you want to get the timing right. 01:06:02.419 --> 01:06:07.059 And what affects that is more or less the current stressing rate. 01:06:07.059 --> 01:06:11.099 So you have kind of these two groups of parameters to play around with, 01:06:11.100 --> 01:06:14.739 and they each impact the problem slightly differently. 01:06:16.820 --> 01:06:24.940 So, yeah, fundamentally, you know, there is free – there are free parameters 01:06:24.940 --> 01:06:28.940 in my model as well. And since we’ve already observed the 01:06:28.940 --> 01:06:34.580 earthquakes, you can kind of tune it to get the right match, like you’re guessing. 01:06:34.589 --> 01:06:39.059 And the one that was most – and I mentioned this – 01:06:39.059 --> 01:06:44.140 has the most variability is probably the background stressing rate. 01:06:44.140 --> 01:06:48.189 It really influences these results quite a bit. 01:06:48.189 --> 01:06:51.279 And unfortunately, we just do not have accurate measurements 01:06:51.279 --> 01:06:55.069 of what the background stressing rate is in this part of the country 01:06:55.069 --> 01:06:59.949 because it can’t really be detected with GPS or something like that. 01:06:59.949 --> 01:07:05.979 So I would like to try and get a closer handle on what that really is. 01:07:05.979 --> 01:07:11.509 And then, yeah, sensitivity to some of the other parameters is 01:07:11.509 --> 01:07:14.680 probably the next step as well. - Thank you. 01:07:14.680 --> 01:07:18.900 - Yeah. Thank you. - Any other questions for Jack? 01:07:19.900 --> 01:07:22.460 If there aren’t, well, thanks very much, Jack. 01:07:22.460 --> 01:07:26.800 I don’t think you opted to go to lunch, but you can find him in Building 3A. 01:07:26.800 --> 01:07:28.960 So let’s thank Jack one more time. - Thanks. 01:07:28.960 --> 01:07:33.920 [ Applause ] 01:07:36.420 --> 01:07:41.820 [ Silence ]