NEAR FIELD AND REGIONAL MODELING OF EXPLOSIONS AT THE DEGELEN TEST SITE J. L. Stevens, N. Rimer, J. R. Murphy, and H.

Xu Maxwell Technologies, Systems Division G. G. Kocharyan and B. A. Ivanov Institute for the Dynamics of the Geospheres, Moscow, Russia Sponsored by the Defense Threat Reduction Agency Contract No. DTRA01-00-C-0050

ABSTRACT Fundamental to the ability to monitor a Comprehensive Nuclear-Test-Ban Treaty is a good understanding of the nuclear explosion source. Considerable research has been done on the nuclear explosion source over the past 30 years, with mixed success. Although empirical and numerical models of explosion sources do a fairly good job of matching observed seismic signals, a good explanation of the physical basis for the explosion source has been elusive. In particular, numerical models of explosion sources developed using laboratory measurements of rock properties fail to predict observed near field ground motion in hardrock. The basic problem is that the strength of the rock measured in the laboratory is much larger than the apparent strength of the rock as determined from the near field ground motion. Furthermore, additional investigation shows that 1) The strength of the rock is not initially low, but rather decreases dynamically as the shock wave passes; and 2) the strength of the rock is reduced to a level well below that predicted for rubbleized rock under hydrostatic dynamic friction. Maxwell Technologies and the Institute for the Dynamics of the Geospheres (IDG) are in the early stages of a research program directed towards improving the capability to predict the seismic source characteristics of underground explosions in rock. This will be accomplished by development of improved dynamic failure models, constrained by a much better data set than has been available in the past. Near field waveforms are only available from a small number of U.S. nuclear tests, and until recently none have been available from the testing program of the former Soviet Union. IDG now has near field records from at least 10 nuclear explosions at the Degelen test site that will be made available for the project. IDG will also provide near source material properties measurements for all of these events. This unique data set will place strong constraints on the numerical modeling work, as well as providing data from a new and important area to augment our previous data sets which have come primarily from the western United States. In addition, the Russians have data from a seismic line located north of the Degelen test site that was maintained with consistent instrumentation for many years during the Soviet testing program. The 9 seismic stations are spread out at approximately even intervals from the Degelen test site out to a distance of about 150 km. IDG will provide seismic data from these stations for all of the events that have near field records. This provides a rare opportunity to observe and model the seismic wavefield of the explosions as they evolve from the near field of the explosion out to regional distances. Objective The objective of this program is to improve the capability to predict the seismic source characteristics of underground explosions in rock. This is to be accomplished by development of improved dynamic failure models which will be constrained by a large unique data set of near field waveforms and parametric data from historic Soviet explosions at the Degelen test site. In addition, we will analyze regional seismic data along a seismic line located north of the Degelen test site that recorded data at 9 stations spaced approximately evenly from the test site to a distance of about 150 km. This project is a collaborative effort between Maxwell Technologies and the Russian Institute for the Dynamics of the Geospheres (IDG).

Research Accomplished 1. Introduction This joint project between Maxwell Technologies and IDG has just started and still in its organizational stages. This paper therefore describes the background, motivation, and research plan for this project. The main objectives of the project are to obtain a better understanding of the explosion source function, and to model the evolution of the explosion source from the close in hydrodynamic region out to near regional distances. This is being accomplished through modeling of a unique data set of near source and near regional data to be provided by IDG. Empirical and numerical models of explosion sources do a fairly good job of matching observed seismic signals, however a good explanation of the physical basis for the explosion source is still not satisfactory. In particular, numerical models of explosion sources developed using laboratory measurements of rock properties fail to predict observed near field ground motion in hardrock. This is illustrated in Figure 1, which shows the near field waveforms from the explosion PILEDRIVER together with waveforms from a finite difference calculation of the near field ground motion made using rock properties measured in the laboratory. The basic problem is that the strength of the rock measured in the laboratory is much larger than the apparent strength of the rock as determined from the near field ground motion. Furthermore, additional investigation shows that 1. 2. The strength of the rock is not initially low, but rather decreases dynamically as the shock wave passes. The strength of the rock is reduced to a level well below that predicted for rubbleized rock under hydrostatic dynamic friction.

A number of solutions to these problems have been proposed over the years, including the effective stress model (Rimer, et al, 1984), and various types of damage models. These models all have the characteristic that the material strength is reduced dynamically to a very low level when it fails. The effective stress model says that the weakness comes from water within the rock matrix, and that the broken rock in effect floats on water that is squeezed out of pores or fractures when the rock fails. Although there are questions about the realism of this physical model, it works fairly well to explain the near field waveforms. Figure 2 shows a comparison of the PILEDRIVER waveforms with waveforms calculated using the effective stress model. The agreement is quite good, particularly at the closer two stations. Furthermore, when the PILEDRIVER solution was scaled to the appropriate yield and compared with other US explosions in granite (Hardhat and Shoal), agreement with the observed waveforms was also quite good (Stevens, et al, 1986). Under a recently completed DTRA contract, Maxwell Technologies worked together with the Russian Institute for the Dynamics of the Geospheres, (IDG) and the University of Southern California, to develop improved micromechanical material models. IDG provided extensive measurements of material properties close to nuclear and chemical explosions both before and after the explosions were detonated (Rimer, et al, 1998). In addition, we implemented the damage mechanics model which was developed by Prof. Charles Sammis at the University of Southern California into Maxwells nonlinear finite difference codes and used this model to simulate the observed explosion damage and a small set of near field waveforms that were also provided by IDG. The results of this work are discussed in detail in the final report (Rimer et al, 1999). The damage model referred to above actually applies to the growth and coalescence of cracks just prior to failure, and does not predict what happens to the rock after failure occurs. Calculations using the damage model followed by a rubbleized model with realistic values of friction did not provide enough strength reduction and thus did not match the near field data. We were more successful in matching the data by dropping the coefficient of friction to very low values (as low as 0.02), but this again leaves the question of what physical mechanism could be responsible for these very low values and corresponding low strength.

35 30 25

B-SL 204 m strong B-SL

6 5.5 5 4.5

16-SLa 470 m strong 16-SLa

Velocity (m/sec)

20 15 10 5 0 -5 -10 0

Velocity (m/sec)
0.1 0.2 0.3 0.4 0.5 Time (seconds) 1401 628 m strong 1401 0.6 0.7 0.8

4 3.5 3 2.5 2 1.5 1 0.5 0 -0.5 -1 0 0.1 0.2 0.3 0.4 0.5 Time (seconds) 1503 863 m strong 1503 0.6 0.7 0.8

2.8 2.6 2.4 2.2 2 1.8 1.6 1.4 1.2 1 0.8 0.6 0.4 0.2 0 -0.2 -0.4 -0.6 0

1.4 1.2 1

Velocity (m/sec)

Velocity (m/sec)
0.1 0.2 0.3 0.4 0.5 Time (seconds) 0.6 0.7 0.8

0.8 0.6 0.4 0.2 0 -0.2 -0.4 0

0.1

0.2

0.3

0.4 0.5 Time (seconds)

0.6

0.7

0.8

Figure 1.

Particle velocity measurements at working point depth from PILEDRIVER compared with data from a calculation using the laboratory strength of granite with no weakening.
B-SL 204 m 6 pile570 B-SL 5.5 5 30 25 4.5 4 16-SLa 470 m pile570 16-SLa

40 35

Velocity (m/sec)

20 15 10 5 0 -5 -10

Velocity (m/sec)
0 0.1 0.2 0.3 0.4 0.5 Time (seconds) 1401 628 m pile570 1401 0.6 0.7 0.8

3.5 3 2.5 2 1.5 1 0.5 0 -0.5 -1 -1.5 0 0.1 0.2 0.3 0.4 0.5 Time (seconds) 1503 863 m pile570 1503 2.5 2 0.6 0.7 0.8

4 3.5 3 2.5

Velocity (m/sec)

Velocity (m/sec)

2 1.5 1 0.5

1.5 1 0.5 0

0 -0.5 -1 -0.5 -1 0 0.1 0.2 0.3 0.4 0.5 Time (seconds) 0.6 0.7 0.8

0.1

0.2

0.3

0.4 0.5 Time (seconds)

0.6

0.7

0.8

Figure 2.

Particle velocity measurements at working point depth from PILEDRIVER compared to numerical simulation pile570 using the effective stress law.

A possible answer initially proposed by Melosh (1979) is acoustic fluidization. The physics behind this mechanism is that during the fracturing process there is a complex dynamic acoustic wavefield that causes high frequency vibrations in the broken rock. These vibrations cause rapidly changing regions of high and low normal stress, and remove the frictional normal stress from parts of the rock as it moves. Consequently parts of the rock are not confined by the frictional stress and in effect have much lower strength. Acoustic fluidization has been used to explain other phenomena such as landslides and craters (Melosh and Ivanov, 1999), which have been similarly difficult to explain because of anomalously low apparent friction.

The numerical modeling to be accomplished in this project will be constrained by a much better data set than has been available in the past. Near field waveforms are only available from a small number of U.S. nuclear tests, and except for those few referenced above, none have been available from the testing program of the former Soviet Union. IDG has near field records from a number of nuclear explosions at the Degelen test site that will be made available for this project. IDG will also provide near source material properties measurements for all of these events. This unique data set will place strong constraints on the numerical modeling work, as well as providing data from a new and important area to augment our previous data sets which have come primarily from the western United States. In addition, IDG has data from a seismic line located north of the Degelen test site that was maintained with consistent instrumentation for many years during the Soviet testing program. Figure 3 shows a map illustrating the location of the test sites and the seismic stations. The 9 seismic stations are spread out at approximately even intervals from the Degelen test site out to a distance of about 150 km. IDG will provide seismic data from these stations for all of the events that have near field records. This provides a rare opportunity to observe and model the seismic wavefield of the explosions as they evolve from the near field of the explosion out to regional distances. The data that IDG has identified to date for this project are listed in Table 1.

9 8 7 6 5 4 3 2 1

Tchinraus Fault Chara Fault BALAPAN Shagan Fault

DEGELEN Kolba-Tchingiz Fault West-Arkkalyik Fault

Figure 3. Map showing the locations of the former Soviet Degelen and Balapan test sites, faults, and seismic stations in the region. We will analyze near field waveforms from the Degelen test site, and seismic data recorded at stations 1 through 9 north of the test site for the same events. 2. Numerical Modeling of Explosions

The parameters used in numerical simulations of underground nuclear explosions are constrained by laboratory material properties tests and by direct observations of ground motion from underground explosions. Quasi-static laboratory tests on small rock samples are used to determine the following parameters: density, elastic moduli, pressure-volume (P-V) relation, material strength, porosity, and water content. Of these quantities, the density, moduli, and P-V relation seem to be quite consistent with insitu rock properties. In brittle hardrocks such as granite, porosity and water content are more variable insitu due to the presence of joints and are therefore less well constrained by laboratory data. Laboratory measurements of strength for brittle hardrocks seem to be very inconsistent with insitu strength as inferred by modeling of underground explosions. In particular, finite difference calculations of ground motion in granite, made using the laboratory measurements of shear strength, have invariably given much narrower particle velocity pulses and much smaller displacements than those measured in the field.

Table 1. Degelen events with near field and/or seismic records. Date 1988/04/22 1988/11/23 1987/07/17 1987/10/16 1989/10/04 1981/07/17 1965/02/04 1964/05/16 1966/03/20 1966/02/13 1973/12/31 1974/12/16 1978/03/26 1980/06/25 1982/12/25 1984/10/18 1968/07/12 1968/12/18 1970/06/28 1980/09/25 1971/12/15 1987/05/06 1987/12/20 1988/10/18 1989/02/17 Number of near field records 2 4 14 4 4 3 4 4 4 4 4 4 5 3 4 4 2 3 2 3 Distance range, m 69-183 280-475 170-900 55-132 45-285 115-310 147-750 150-600 300-600 350-600 110-230 100-200 76-270 155-310 90-190 35-110 80-90 200-510 200-240 100-160 Number of seismic records 3 9 7 8 8 8 9 9 9 9 Distance range, km 57-81 14-83 14.5-84 20-76 16-85 15-80 14-83 13-80 15-84 14-80

9 9 9 9 9

7.5-85 13-83 13-83 11-77 10-77

Constitutive models have been developed (Rimer and Lie, 1982, Rimer, et al., 1984) which attribute this weaker behavior of insitu granite and other rocks under explosive loading to ground-motion-induced rock damage or pore fluid pressure increases. These models have been shown by Day, et al, (1983) and Rimer, et al, (1986), (1987), using one- and two-dimensional calculations, to accurately reproduce the pulse broadening and most other aspects of the ground motion data from U.S. events in granite (PILEDRIVER, for example). The models for granite have since been used to simulate the seismic source function from explosions conducted in the Soviet Union (Stevens, et al, 1991, for example). Figure 2 compares the particle velocity measurements at working point depth from the 62 kt PILEDRIVER event with the results of the pile570 numerical simulation, made using the effective stress model discussed in Rimer, et al (1984, 1998, 1999). Numerical results from pile570 are also in good agreement with the measured cavity radius and with the estimated seismic source function. For this simulation, effective stress model parameters were calibrated to best match the velocity peak and pulse shape at the closest-in gauge station (B-SL). This required a rapid buildup of pore pressure during the loading, leading to a large strength reduction very near the propagating shock front. Note that the pulse shapes at all four stations are rather consistent, with all including a long duration shallow negative velocity pulse. The peak velocities at the two smallest ranges agree very well with the simulation, but the measured peaks at the two larger ranges are a factor of two or more lower than the calculated peaks. The two closest PILEDRIVER gauges were located along roughly a 180 degree different azimuth than the other more distant gauges. However, a connection between possible site anisotropies and the measured ground motions has never been established. It should be noted that the 5.9 kt HARDHAT event, which was detonated in similar rock, near the later PILEDRIVER event, but at a shallower depth of burial, gave particle velocity pulses which were more similar to the PILEDRIVER pulses at the larger ranges. (See Rimer, Stevens, and Day, 1987, for a comparison between the results of the pile570 simulation and the HARDHAT measurements.)

The constitutive models used in these simulations are phenomenological in nature and do not explicitly account for the dynamic response of the insitu fractures in the source environment. That is, crystalline rock massifs in the earths crust are heterogeneous over a wide range of scales from micro-cracks in grains at the micron scale to faults and joints which span many tens of kilometers. It has been observed that this fracture heterogeneity is responsible for nonlinear effects near the source of underground explosions which affect seismic coupling and the resultant seismic waveform. Very near the source, pre-existing cracks nucleate new fractures which granulate the rock, thereby significantly reducing its shear strength and elastic moduli. At greater distances, a network of radial tension cracks are formed which also degrade the mechanical properties of the rock. Even at large distances, movement on pre-existing faults and joints can produce significant deviations from numerical computer simulations which assume a simply connected elastic continuum. Over the past several years, Professor Sammis and his associates at the University of Southern California have been working to understand the physics and micro-mechanics describing the nucleation and growth of new fractures from pre-existing cracks in crystalline rock under the compressive loading states inherent in laboratory tri-axial compression tests. Similar compressive loading states characterize the shock wave loading near the source of an underground explosion. This work (see Sammis and Ashby, 1988, Ashby and Sammis, 1990, and Sammis, 1991) has resulted in a damage mechanics formulation which is suitable for incorporation into numerical simulation codes. The model introduces a damage parameter, related to the increase in flaw size from its pre-shot average value. Damage accumulates as the flaws extend during the compressive loading and reaches some maximum value at which the rock fails unstably. The model gives no guidance to rock behavior after unstable failure. Incorporation of this damage mechanics formulation into a finite difference code, and subsequent testing and validation of the revised models, were the main topics of Rimer, et al (1998, 1999). Adaptation of the Sammis damage mechanics model for use in the ground motion simulation code is summarized in Rimer, et al (1998). The key input parameter to this model, in addition to the usual elastic moduli, is the average initial flaw size at the rock site. Since unstable compressive failure of a rock element is calculated using the adapted model to occur relatively early in the dynamic motions of interest here, i.e., usually near the propagating shock front, additional modeling was incorporated to complete the description of the stress field after this failure. Limiting the magnitude of deviatoric stresses in a failed rock element through the use of a standard friction law was shown to not provide sufficient strength reduction to simulate the ground motion measurements. In particular, calculated particle velocity time histories were still much narrower than the measured pulses for all reasonable choices of initial flaw size. The additional strength reduction required to sufficiently broaden the particle velocity pulses was obtained by using a shear damage model originally developed for soft rocks such as tuff, described in Rimer and Proffer (1991). As discussed in Rimer et al (1999), this shear damage model was applied here only for rock elements that had experienced the onset of unstable compressive failure. Thus, a rock element was allowed to undergo significant damage before application of this treatment. The post-failure shear damage model performs an interpolation between the standard coefficient of friction of 0.60 and much lower effective friction values of 0.02-0.20 in the friction law used to limit deviatoric stresses. This linear interpolation is based on the maximum shear strain experienced by the failed rock element. A possible physical mechanism for this strength reduction is the acoustic fluidization concept proposed by Melosh (1979) to explain the low strengths (or very low angles of internal friction) apparent in a number of geologic processes, such as seismic faulting, impact crater slumping, and long runout landslides. The main concept of this proposed mechanism is that sufficiently strong acoustic waves in rock debris can momentarily relieve the static overburden in limited regions of the rock mass, thus allowing sliding to occur in the unloaded region. If this happens frequently enough, flow of the entire rock mass results. In terms of the explosively-induced ground motions of interest here, only the mass of failed or fractured rock, as given by the Sammis damage mechanics formulation, would have the potential for acoustic fluidization. The dynamic fracturing process, by itself, provides enough energy to generate sufficiently strong oscillations post-failure to reduce the effective normal stress/effective friction to the low values required in the weak core region near the cavity. We have not yet developed a first principles relationship between the effective friction coefficient in a (failed) rock element and calculated quantities such as perhaps total energy, kinetic energy, or accumulated damage. The effective friction coefficient for complete damage in the calculations to date is simply an input number, taken from the examples given by Melosh. Most of his examples, however, correspond to far less dynamic situations. For the

explosive induced ground motions of interest here, an effective viscosity approach, rather than the effective friction approach used so far, may prove to be a more physically reasonable treatment. A series of calculations were made using the Sammis modeling, quantifying the effect of model parameters on particle velocity pulse shapes, cavity radius, and RDP. The simulation which best fits all of the PILEDRIVER ground motion data is Run PD10. Comparisons with measured particle velocities are shown in Figure 4. In contrast with the results of pile570, made with the effective stress model, PD10 provides a much better fit to the PILEDRIVER data at the two larger ranges, while underestimating peak velocities at the closest ranges. Subsequent analyses showed that the timing of the strength reductions for the two models were somewhat different, with the effective stress model providing an earlier reduction than the present damage model. However, calculated cavity radius and RDP with the new model are in as good agreement with the measurements as were those with the more traditional model.
35 30 25 PD10 B-SL 6 5.5 5 4.5 4 PD10 16-SL

Velocity (m/sec)

Velocity (m/sec)
0 0.1 0.2 0.3 0.4 0.5 Time (seconds) 0.6 0.7 0.8

20 15 10 5 0 -5 -10

3.5 3 2.5 2 1.5 1 0.5 0 -0.5 -1 0 0.1 0.2 0.3 0.4 0.5 Time (seconds) 0.6 0.7 0.8

2.6 2.4 2.2 2 1.8 1.6 PD10 1401

1.8 1.6 1.4 1.2 PD10 1503

Velocity (m/sec)

1.4 1.2 1 0.8 0.6 0.4 0.2 0 -0.2 -0.4 -0.6 0 0.1 0.2 0.3 0.4 0.5 Time (seconds) 0.6 0.7 0.8

Velocity (m/sec)

1 0.8 0.6 0.4 0.2 0 -0.2 -0.4 0 0.1 0.2 0.3 0.4 0.5 Time (seconds) 0.6 0.7 0.8

Figure 4. Calculated particle velocity pulses at four ranges for Run PD10 and PILEDRIVER measurements.

The PD10 shear damage model begins the reduction in friction coefficient from the standard coefficient of 0.60 once failed elements experience a maximum shear strain of 0.1% or greater. The reduction to the lower effective friction coefficient of 0.02 for failed rock elements is completed for 3% shear strain or greater. In practice, calculated pulse widths were found not to be very sensitive to an order of magnitude increase in the value of shear strain used to initiate the interpolation since the onset of unstable failure usually occurred at much higher shear strains in the simulations. However, particle velocity pulse widths were sensitive to the choice of shear strain magnitude required for full strength reduction, with lower values of this parameter resulting in longer pulse duration and larger RDP. It is important to emphasize that while the shear damage model (or some other postfailure strength reduction model) is crucial to successful simulation of the ground motion data, it is the micromechanical damage mechanics model which primarily determines the size of the central core of weakened rock. Rimer, et al (1998) also contained extensive data sets on the density and morphology of the near-source fracture sets generated by several large chemical and nuclear explosions detonated in crystalline rock in the Soviet Union. These data, obtained from IDG, consist of comparisons of pre-shot and post-shot fracture densities, aperatures, lengths, and/or P-wave velocities versus scaled range (m/kt1/3) from 11 explosions (two of these were multiple events).

Included with these insitu measurements and observations are average material properties data (densities, porosities, ultrasonic velocities, and unconfined strengths) obtained from laboratory tests on pre-shot core samples from each explosion site. For several of these events, laboratory measurements on post-shot core samples quantified explosioninduced increases in porosity near the cavities. Based on these data sets, the scientists at IDG have also developed empirical relations which scale rock damage with distance and explosive yield. These data place important constraints on nonlinear source models and were used to help validate the damage mechanics formulation. More recently (Kocharyan, et al, 1998), IDG has also provided near-field dynamic ground motion data, more specifically, radial particle velocities versus time at several ranges from three of the underground explosions at Degelen Mountain for which fracture data were presented earlier. In Rimer, et al (1999), comparisons are made between the PILEDRIVER ground motion data and the Degelen data (scaled to the PILEDRIVER explosive yield). Following are some results of these comparisons: The Degelen ground motions in general do not show the long duration, shallow negative velocity pulses characteristic of PILEDRIVER and other U.S. events in granite (and other rocks). Peak velocities from Degelen events are in good agreement with the PILEDRIVER records at the larger scaled ranges, but are significantly smaller than those at the two closest ranges shown in Figures 1 and 2. Positive velocity pulses from the Degelen events, when scaled to PILEDRIVER yield, are of slightly smaller duration than those from PILEDRIVER. This appears consistent with the smaller post-shot cavities reported by IDG, 7-9 scaled meters, versus 10-11 scaled meters for PILEDRIVER.

In contrast to the differences in measured cavity radii and near-field ground motions between the two granite test sites, the seismic evidence (i.e., mb/yield relations) strongly suggests that the seismic source coupling at the Degelen site is very comparable to that of PILEDRIVER. (See Murphy, 1993.) In general, the granite at the U.S. test site has quite similar material properties to that of the Degelen site. Laboratory measurements for the U.S. rock do show a slightly higher mass density and wave speed, and much lower porosities, 0.1-0.3%, compared to 2-3% for Degelen material. Measured strengths are very similar. Traditional ground motion calculations made with these different material properties do not explain the differences between the ground motion measurements at the two sites. These ground motion differences may be the result of differences in the insitu structures of the two granite sites, such as joint, faults, and other fractures, their orientation and spacing. For the Degelen explosions, where both dynamic ground motion and fracture extent data are available, model parameters were varied to match as much of the data set as possible in order to understand the constraints of these data. The best Degelen calculation, Run DE12, on average, provided a reasonably good simulation of the IDG measurements/observations of particle velocity pulses, cavity radius, and RDP, as well as the increases in porosity and decreases in P-wave velocity as a function of post-shot range from the explosion working point. Figure 5 compares the particle velocity pulses from Run DE12 with the measurements from three Degelen explosions, labeled N1, N8, and N9, all scaled to the PILEDRIVER yield. The simulation is in good agreement with measured peak velocities at five of the six ranges. However, DE12 also shows a second peak or plateau at all six ranges, while such a second signal is readily apparent for only two of the somewhat inconsistent gauge records (although oscillations on three of the other gauge records can be construed as indications of second signals). Large negative pulses, comparable to the simulation, are measured only at some of the larger scaled ranges.

22 20 18 DE12 N1-207m

4 3.5 3 DE12 N9-415m

3.5 3 2.5
Velocity (m/sec)

DE12 N8-463m

16
Velocity (m/sec) Velocity (m/sec)

14 12 10 8 6

2.5 2 1.5 1 0.5

2 1.5 1 0.5 0 -0.5 -1

4 2 0 -2 0 0.1 0.2 0.3 0.4 0.5 Time (seconds) 0.6 0.7 0.8 0 -0.5 -1 0 0.1 0.2 0.3 0.4 0.5 Time (seconds) 0.6 0.7 0.8

0.1

0.2

0.3 0.4 0.5 Time (seconds)

0.6

0.7

0.8

3.5 3 2.5
Velocity (m/sec) Velocity (m/sec)

1.8 DE12 N9-519m 1.6 1.4 1.2 1 0.8 0.6 0.4 0.2 0 -0.2 -0.4
Velocity (m/sec)

1.4 DE12 N9-705m 1.2 1 0.8 0.6 0.4 0.2 0 -0.2 -0.4 0 0.1 0.2 0.3 0.4 0.5 Time (seconds) 0.6 0.7 0.8 -0.6 0 0.1 0.2 0.3 0.4 0.5 Time (seconds) 0.6 0.7 0.8 DE12 N8-741m

2 1.5 1 0.5 0 -0.5 -1

-0.6 0 0.1 0.2 0.3 0.4 0.5 Time (seconds) 0.6 0.7 0.8 -0.8

Figure 5. Calculated particle velocity pulses from Run DE12, which includes dilatancy, compared with the Degelen measurements.

In the calculations, the second signal is the result of a dilatant model relating decreases in bulk modulus to increases in the calculated damage parameter. This dilatant model and the non-dilatant model for decreases in shear modulus with damage accumulation are used to simulate the large observed P-wave velocity decreases near the cavity. The dilatant bulk modulus reductions also lead to calculated decreases in mass density (increases in porosity) very near the cavity. P-wave velocities from Run DE12 have been plotted versus post-shot range in Figure 6 along with the calculated increases in porosity in the rock near the cavity as a result of the explosion. (An initial porosity of 3% is assumed. Inclusion of this porosity, which was assumed to be crushable only at stress levels in the tens of kilobars, had no significant effect on the calculated pulses.) Comparisons are also shown with the data as reported by IDG. These include both the empirical fit to the post-shot P-wave velocity measurements given by IDG (scaled here to 62 kt) and reported variations in P-wave velocity and porosity data in the most damaged regions of rock surrounding the cavity.

25

Porosity (%)

20 15 10 5 0 0 100 200 300 400 500 600 700 800 900

5000

P-wave Velocity (m/s)

4000

3000

2000

1000

Range of IDG Data IDG Fit DE12

0 0 100 200 300 400 500 600 700 800 900

Post-Shot Range (m)

Near the cavity, the calculated minimum P-wave velocity of 2335 m/s was well within the variation of the measurements in this region. The slight increase in calculated P-wave velocity close to the cavity wall is not due to increases in calculated bulk or shear modulus. Rather, it results from decreased rock density (increased porosity) in this dilatant region. Further from the cavity, calculated P-wave velocity increases are also well within the data scatter which underlies the IDG empirical fit. The calculation however, shows a plateau in P-wave velocity, or more precisely, a plateau in the accumulated damage which determines the shear and bulk moduli. This plateau, which extends inward from the end of the failed region, encompasses all of the rock elements for which tensile failure was calculated before the onset of unstable compressive failure. At larger ranges, tensile failure was calculated to occur out to approximately the end of the partially damaged region, which in this computational model, also marks the end of post-shot P-wave velocity decreases. Porosity is calculated to increase from its nominal 3.0 percent initial value only near the cavity. Although the results agree fairly well with the data, Run DE12 shows on average somewhat smaller increases at a given range, with very large porosity (42.5%) calculated only in the grid zone adjacent to the cavity. (As a result of grid rezoning during the calculation, this grid zone actually encompasses the first eight of the initial grid zones.) It is very possible that the very simple model used to introduce dilatancy here is inadequate. However, it is much more likely that the insitu rock near the cavity has experienced significant asymmetric block motion which cannot be reproduced using our spherically symmetric code. The calculated RDP and cavity radius from Run DE12 are in good agreement with the observations. Dilatancy in the calculations acts to increase RDP and decrease cavity radius. Although we found it difficult, it may also be possible to simulate the observed RDP and cavity radius without introducing dilatancy by changing other parameters in our constitutive model. Without some dilatancy model however, rock porosity does not increase near the cavity.

In summary, it has been shown in Rimer et al (1999) that most of the available ground motion and damage data from Degelen Mountain explosions in granite may be simulated reasonably well using spherically symmetric calculations with the modified Sammis micro-mechanical damage mechanics model. For reasonable initial flaw size estimates, however, this model predicts compressive failure to roughly a factor of two larger range than the range of intensive failure estimated by IDG from their data analysis. This is not understood as yet. It should be emphasized that simulation of the Degelen ground motion and damage data requires the use of model parameters which are considerably different from those used to simulate the ground motion data from U.S. PILEDRIVER explosion in similar granite. For example, Run DE12 used a damaged effective friction of 0.20 rather than 0.02, and larger shear strain magnitudes (6%) to reach full damage in order to simulate the narrower particle velocity pulses. It is possible that other combinations of model parameters might provide as good a simulation of these Degelen data as those used in Run DE12. Nevertheless, it has been demonstrated that the following model features are important in simulating these data: A post-failure damage model gives the central core of weakened rock required to sufficiently widen the particle velocity pulses. The strength of this weakened core is determined by parameters of this damage model. Parameters of the micro-mechanical damage mechanics model determine the radial extent of this weakened core. Dilatant reductions in bulk modulus are required to simulate the P-wave velocity and porosity measurements.

Key Words: Rock Mechanics, Degelen, explosion, numerical modeling REFERENCES Ashby, M.F. and C.G. Sammis (1990), The Damage Mechanics of Brittle Solids in Compression, PAGEOPH, Vol. 133, No. 3. Day, S. M., N. Rimer and J. T. Cherry (1983), Surface Waves from Underground Explosions with Spall: Analysis of Elastic and Nonlinear Source Models, BSSA, 73, pp. 247-264. Kocharyan, G.G., B.A. Ivanov, and V.N. Kostyuchenko (1998), Rock Damage and Estimated Seismic Source Parameters of Underground Explosions in Hardrock, Russian Acad. Sci., IDG Final Report on Subcontract 9-11439-KLW, December. Melosh, H.J. (1979), Acoustic Fluidization: A New Geologic Process?, J. Geophys. Res., 84, 7513, December. Melosh, H.J., and B.A. Ivanov (1999), Impact Crater Collapse, Annual Rev. Earth Planet Science, 27, pp385-415. Murphy, J.R. (1993),"Comment on 'Q for short-period P waves: is it frequency dependent?' by Alan Douglas", Geophys. J. Int., 113, p535-540. Rimer, N. and K. Lie (1982), Numerical Simulation of the Velocity Records from the SRI Grout Spheres Experiments, S-CUBED Topical Report DNA-TR-82-54, September. Rimer, N., S. M. Day, G. A. Hegemier, H. E. Read, S. K. Garg, and S. Peyton (1984), Effect of Pore Fluid Pressure on Explosive Ground Motions in Low Porosity Brittle Rocks, S-CUBED Final Report DNA-TR-85-245, July. Rimer, N., J. L. Stevens, and K. H. Lie (1986), Effects of Pore Pressure and Fractures on Ground Motion in Granite, S-CUBED Final Report DNA-TR-86-227, June. Rimer, N., J. L. Stevens, and S.M. Day (1987), Effects of Pore Pressure, Fractures, and Dilatancy on Ground Motion in Granite, S-CUBED Final Report AFGL-TR-87-0136, April. Rimer, N. and W. Proffer (1991), Containment Phenomenology Using a New Shear-Strain-Based Computational Damage Model for Tuff, S-CUBED Technical Report SSS-DTR-91-12612, September, also SSS-DTR-9011961, October. Rimer, N., K. Lie, J.L. Stevens, J.R. Murphy, and G.G. Kocharyan (1998), A Micro-Mechanical Damage Model for Estimating Seismic Source Characteristics of Underground Explosions in Hardrock, Maxwell Technologies Technical Report MFD-DTR-98-15985, January. Rimer, N., J.L. Stevens, J.R. Murphy, and G.G. Kocharyan (1999), Estimating Seismic Source Characteristics of Explosions in Hardrock Using a Micro-Mechanical Damage Model, Maxwell Technologies Final Report MTSD-DTR-99-16423, July.

Sammis, C.G. and M.F. Ashby (1988), The Damage Mechanics of Brittle Solids in Compression, AFGL-TR-880160, July. Sammis, C.G. (1991), A Damage Mechanics Source Model for Underground Nuclear Explosions, PL-TR-912103, August. Stevens, J. L., N. Rimer, and S. M. Day (1986), "Constraints on modeling of underground explosions in granite," SCUBED annual report to Air Force Geophysics Laboratory, AFGL-TR-86-0264, SSS-R-87-8312, October. Stevens, J.L., T.G. Barker, E.J. Halda, J.R. Murphy, and N. Rimer (1991), Simulation of Teleseismic Body Waves from Underground Explosions Beneath Sloping and Flat Surfaces Using Nonlinear Axisymmetric Finite Difference Calculations, S-CUBED Technical Report DNA-TR-91-5, July.