EFFECTS OF TILTED AND FAULTED STRATA ON 

SEISMIC GROUND MOTION 

Prepared for 

U.S. Nuclear Regulatory Commission 

Contract NRC–02–07–006 

Prepared by 

Amitava Ghosh 

Simon Hsiung

Center for Nuclear Waste Regulatory Analyses 

San Antonio, Texas 

September 2011 
ii 

ABSTRACT 

This report documents the results of estimating seismic motion at the ground surface from an 

earthquake. The earthquake motion amplifies while traveling upward through the geologic 

medium. This amplification is dependent on the properties of the geologic medium. Reliable 

estimation of the ground motion is necessary for designing surface structures. Although 

traditionally a one-dimensional analysis is conducted to estimate the ground motion, in this 

study two-dimensional analyses are conducted to estimate the ground motion from an 

earthquake. The primary objective of this study is to assess the potential spatial variation of 

ground response from seismic waves in fractured and/or layered deposits, especially with a 

thick soil cover having varying thickness and extent. The purpose of the study is to develop an 

understanding of seismic response in a complex medium with varying soil thickness and extent. 

The UDECTM code is used to simulate a fractured medium without material degradation, and the 

FLAC code is used to simulate a layered medium with material degradation. Results using 

UDEC show that the amplification of seismic motion is small, generally limited to approximately 

1.5 times the input motion. However, with material degradation incorporated in the simulation 

using FLAC, the amplification could be large—several times the input motion. In addition, the 

estimated ground motion at the surface shows resonant behavior with long duration oscillations. 

These phenomena have been observed in different earthquake records, especially over thick 

soil deposits similar to those used in this study. The results are preliminary and require 

further investigation. 
iii 

CONTENTS 

Section Page 

ABSTRACT ................................................................................................................................... ii 

FIGURES ..................................................................................................................................... iv 

TABLES ........................................................................................................................................ v 

ACKNOWLEDGEMENTS ............................................................................................................ vi 

1 INTRODUCTION ............................................................................................................ 1-1

1.1 Background and Purpose ................................................................................... 1-1 

1.2 Objectives and Scope ........................................................................................ 1-1 

2 DYNAMIC SOIL/ROCK BEHAVIOR .............................................................................. 2-1 

3 SITE RESPONSE ANALYSIS ........................................................................................ 3-1 

4 ANALYSIS OF SURFACE RESPONSE ........................................................................ 4-1 

4.1 Fractured Rock Mass Model: UDEC Analysis ................................................... 4-1 

4.1.1 Rock Mass Model in UDEC .................................................................... 4-1 

4.1.1.1 Results and Discussion .......................................................................... 4-7 

4.2 Horizontal and Inclined Strata Models: FLAC Analysis ................................... 4-11 

4.2.1 Strata Model ......................................................................................... 4-11 

4.2.2 Model With Horizontal Strata Layers .................................................... 4-14 

4.2.3 Model With Inclined Strata Layers ........................................................ 4-18 

4.2.4 Basin Model: FLAC Analysis ............................................................... 4-18 

5 SUMMARY AND CONCLUSIONS ................................................................................. 5-1 

6 REFERENCES ............................................................................................................... 6-1 

APPENDIX—HYSTERETIC DAMPING 
iv 

FIGURES 

Figure Page 

2-1 Typical Shear Stress and Shear Strain Relationship ..................................................... 2-2 

2-2 Variation of Normalized Shear Modulus With Shear Strain Using Data for Upper 

Mean Tuff and Upper Mean Alluvium ............................................................................ 2-2 

2-3 Variation of Material Damping Ratio With Shear Strain Using Data for Upper Mean 

Tuff and Upper Mean Alluvium ...................................................................................... 2-3 

4-1 Particle Acceleration Versus Time of the Coalinga1 Earthquake Used In 

FLAC Simulations .......................................................................................................... 4-2

4-2 Input Motion in Terms of Particle Velocity Versus Time Used in UDEC Simulations ..... 4-2 

4-3 Cross Section of Strata Used in Seismic Response Modeling Using 

the UDEC Code ............................................................................................................. 4-3

4-4 The UDEC Model of the Strata Cross Section With Monitoring Points 1–9 ................... 4-5 

4-5 Comparison of Ground Motion at Point 1 With Input Motion (Point 6) 

(Case 1) ......................................................................................................................... 4-8 

4-6 Comparison of Ground Motion at Points 1–5 With Input (Point 6) 

(Case 2) ......................................................................................................................... 4-8 

4-7 Comparison of Ground Motion at Points 1 and 5 With Input (Point 6) 

(Case 2) ....................................................................................................................... 4-10 

4-8 Comparison of Ground Motion at Points 1 and 5 With Input (Point 6) 

(Case 3) ....................................................................................................................... 4-10 

4-9 UDEC Model of the Fractured and Jointed Rock Mass, One Fracture Set Dips 

 Toward Left and Another Toward Right ....................................................................... 4-12 

4-10 Comparison of Ground Motion at Point 5 (on Soil) With Input (Point 6) 

(Case 4) ....................................................................................................................... 4-13 

4-11 Horizontal Strata Model Used in FLAC Analysis .......................................................... 4-13 

4-12 Inclined Strata Model Used in FLAC Analysis .............................................................. 4-14 

4-13 Estimated Velocity History at the Top of the Model for First 18 Seconds With the 

Default Material Degradation Model [y-axis: Velocity (m/s), x-axis: Time (s)] ............ 4-16 

4-14 Comparison of Acceleration at the Top of the Model With Input Motion at the Base 

With the Default Material Degradation Model [y-axis: Acceleration (m/s2

), 

x-axis: Time (s)] ........................................................................................................... 4-16 

4-15 Spectral Acceleration at the Top Surface of the Horizontal Layer Model as a 

Function of Frequency ................................................................................................. 4-17 

4-16 Comparison of Acceleration at the Top Surface With Input Motion at the Base with 

Sig3 Material Degradation Model [y-axis: Acceleration (m/s2

), x-axis: Time (s)] ........ 4-17 

4-17 Comparison of Acceleration at the Top Surface With Input Motion at the Base with 

Sig4 Material Degradation Model [y-axis: Acceleration (m/s2

), x-axis: Time (s)] ........ 4-19 

4-18 Comparison of Acceleration Time Histories at Left and Right Corners of the Inclined 

Strata Model Shown in Figure 4-12 [y-axis: Acceleration (m/s2

), x-axis: Time (s)] ..... 4-19 

4-19 Basin Model Used With Monitoring Points (Notes All Three Rock Blocks Have 

Same Properties .......................................................................................................... 4-20

4-20 Comparison of Seismic Wave Amplification Between Rock and Soil at the Center 

of the Basin (Impedance Ratio 5.0) [y-axis: Acceleration (m/s2

), x-axis: Time (s)] ..... 4-21 

4-21 Comparison of Seismic Wave Amplification Between Rock and Soil at the Edge of 

the Basin (Impedance Ratio 5.0) [y-axis: Acceleration (m/s2), x-axis: Time (s)] ........ 4-21 

4-22 Comparison of Seismic Wave Amplification Across the Soil Basin 

(Seismic Impedance Ratio 2.79) [y-axis: Acceleration (m/s2

), x-axis: Time (s)] ......... 4-22 
v 

FIGURES (continued) 

Figure Page 

4-23 Comparison of Seismic Wave Amplification Across the Soil Basin 

(Seismic Impedance Ratio 3.05) [y-axis: Acceleration (m/s2

), x-axis: Time (s)] ......... 4-22 

4-24 Comparison of Seismic Wave Amplification Across the Soil Basin 

(Seismic Impedance Ratio 5.0) [y-axis: Acceleration (m/s2

), x-axis: Time (s)] ........... 4-23 

4-25 Comparison of Spectral Acceleration Across the Soil Basin (Seismic Impedance 

Ratio 3.05) .................................................................................................................... 4-24 

4-26 Comparison of Spectral Acceleration Across the Soil Basin (Seismic Impedance 

Ratio 5.0) ...................................................................................................................... 4-24 
vi 

TABLES 

Table Page 

4-1 Material Properties Used in UDEC Simulations ............................................................. 4-5 

4-2 Properties of Faults and Strata Contacts Used in UDEC Simulations ........................... 4-6 

4-3 Properties of Joints Used in UDEC Simulations (Case 4) .............................................. 4-6 

4-4 Results From Four Simulations Using UDEC................................................................. 4-9 

4-5 Fitted Parameters for Default Model Representing Material Degradation With 

Shear Strain ................................................................................................................. 4-15 

4-6 Fitted Parameters for Sig3 Model Representing Material Degradation With 

 Shear Strain ................................................................................................................. 4-15 

4-7 Fitted Parameters for Sig4 Model Representing Material Degradation With 

Shear Strain ................................................................................................................ 4-15 

4-8 Material Properties Used in Simulations ...................................................................... 4-20 

4-9 Impedance Ratio Between Rock and Soil Used in Simulations ................................... 4-20 
vii 

ACKNOWLEDGMENTS 

This report was prepared to document work performed by the Center for Nuclear Waste 

Regulatory Analyses (CNWRA®) for the U.S. Nuclear Regulatory Commission (NRC) under 

Contract No. NRC–02–07–006. The activities reported here were performed on behalf of the 

NRC Office of Nuclear Material Safety and Safeguards, Division of High-Level Waste 

Repository Safety. The report is an independent product of CNWRA and does not necessarily 

reflect the views or regulatory position of NRC. 

The authors gratefully acknowledge B. Dasgupta for his technical review, B. Sagar for his 

programmatic review, L. Mulverhill for her editorial review, and A. Ramos for his 

secretarial support. 

QUALITY OF DATA, ANALYSES, AND CODE DEVELOPMENT 

DATA: All CNWRA-generated data contained in this report meet quality assurance 

requirements described in the Geosciences and Engineering Division Quality Assurance 

Manual. Sources of other data should be consulted for determining the level of quality of 

those data. 

ANALYSES AND CODES: UDEC Version 5.0 and FLAC Version 5.0 (Itasca Consulting Group, 

Inc., 2011; 2005) have been used in the analyses presented in this report. These codes are 

controlled under the CNWRA quality assurance program. Analyses described in this report are 

documented in CNWRA Scientific Notebook 889 (Ghosh, 2007). 

REFERENCES 

Ghosh, A. “Characterization of Effects of Faulted and Inclined Strata on Ground Response at 

the Proposed Site for Surface Facilities at the Potential Repository at Yucca Mountain.” 

Scientific Notebook 889. San Antonio, Texas: Center for Nuclear Waste Regulatory Analyses. 

33 p. 2007. 

Itasca Consulting Group, Inc. “UDEC Version 5.0 User’s Guide.” Minneapolis, Minnesota: 

Itasca Consulting Group, Inc. 2011. 

_____. “FLAC Fast Lagrangian Analysis of Continua User’s Guide (Version 5.0).” 

Minneapolis, Minnesota: Itasca Consulting Group, Inc. 2005. 
1-1 

1 INTRODUCTION 

1.1 Background and Purpose 

An estimate of the seismic ground motion at the surface of the Earth is necessary to (i) develop 

design response spectra, (ii) evaluate dynamic stresses and strains for assessing liquefaction 

potential, and (iii) determine earthquake-induced force for assessing stability of Earth-retaining 

structures (Kramer, 1996). Generally, empirical ground motion relationships (attenuation 

relationships), based on characteristics of previously recorded earthquakes, are used to 

estimate the seismic motion at the bedrock of a given site. The ground response analysis is 

used to assess the response of the soil or rock medium subjected to the seismic motion at 

the bedrock. 

Over the years, a number of different techniques have been developed to estimate the surface 

motion from an earthquake. Both linear and nonlinear models have been used to approximate 

the material behavior when subjected to an earthquake motion. Analyses also have been 

conducted in one, two, and three dimensions. The Center for Nuclear Waste Regulatory 

Analyses (CNWRA®) conducted analyses to simulate the ground response using a 

one-dimensional (1D) model with linearized material behavior (Hsiung, et al., 2007). CNWRA 

staff also conducted additional two-dimensional (2D) analyses using both UDEC Version 5.0 

(Itasca Consulting Group, 2011) and FLAC Version 5.0 (Itasca Consulting Group, 2005) codes 

to assess the spatial variation of the ground response at the surface. 

This report documents the analysis conducted so far in estimating the ground motion at different 

locations at the surface when the seismic waves from a nearby strong motion earthquake, 

Coalinga1 (Hsiung, et al., 2007), propagate through the rock and soil layers. The rock and 

layers form a rock mass fractured by a number of faults and joint sets. Additionally, materials in 

all rock and soil layers sustain significant deterioration of the shear modulus (modulus of rigidity) 

when the waves propagate through them. 

1.2 Objectives and Scope 

The objectives of the analyses reported in this document are to assess the potential spatial 

variation of ground response at the surface when the seismic waves from the underlying 

bedrock propagate through a fractured and/or layered deposit with material variability and a 

thick soil cover. In addition, the shear modulus of these materials deteriorates when subjected 

to high shear strain from the propagating seismic waves. This report documents the results 

obtained from the analyses for understanding the phenomena that take place when these 

materials are subjected to large-amplitude seismic waves from a strong-motion earthquake. 

The analyses were conducted using various readily available data associated with the proposed 

geologic repository site at Yucca Mountain, Nevada. However, note that the results reported 

here are preliminary approximations and do not pertain to any licensing action; the geometry of 

the rock mass and, in some cases, material properties were significantly simplified for ease 

in analysis. 

Behavior of soil and rock under dynamic load is discussed in Chapter 2. Traditional site 

response analysis is discussed in Chapter 3. Discussion also includes spatial variation and 

resonant behavior observed in the ground motion in 2D analyses and in actual earthquake 

records. Results obtained in this study using UDEC and FLAC codes are given in Chapter 4. 

Chapter 5 summarizes this study and discusses proposed future work. 
2-1 

2 DYNAMIC SOIL/ROCK BEHAVIOR 

Generally, soils including soft rocks exhibit nonlinear response under shear loads especially at 

larger strain amplitudes. The shear modulus of the material decreases as the shear strain 

increases with an accompanying increase in material damping. These changes in material 

properties affect the amplitude and frequency content of the resulting ground motion. The 

constitutive model for the material response must account for these changes in the material 

properties to appropriately model the ground response. 

It is well known that soil and rock exhibit nonlinear hysteresis behavior, as shown in Figure 2-1 

(Kramer, 1996; Bechtel SAIC Company, LLC, 2004; EPRI, 1993). The inclination of the 

hysteresis loop is dependent on the soil stiffness and expressed by the tangent shear modulus 

Gtan. Gtan varies throughout the loading cycle, its average value over the loop is approximated 

by the secant shear modulus Gsec 

ܩ௦௘௖ ൌ ߬௖

௖ߛ

(2-1) 

where τc and γc are the magnitudes of shear stress and shear strain, respectively. 

The breadth of the hysteresis loop represents the energy dissipated in a loading cycle and is 

described by the damping ratio ξ

 ൌ 1

ߨ2

௣௢௢௟ܣ

௖ߛ௖௘௦ܩ

ଶ

(2-2) 

where Aloop is the area of the hysteresis loop. Gsec and ξ are called the equivalent nonlinear 

material parameters of the soil. 

Secant modulus Gsec of a soil varies with the cyclic shear strain amplitudes; at low shear strain, 

Gsec is large but it decreases as the strain amplitude increases. The secant modulus at zero 

cyclic shear strain amplitude has the highest value Gmax. The secant modulus normalized by 

this highest value Gmax as a function of the shear strain for tuff rock (Upper Mean Tuff) and 

alluvium (Upper Mean Alluvium) at the Yucca Mountain site is shown in Figure 2-2 

(Bechtel SAIC Company, LLC, 2004) and is consistent with the typical curves observed for soils 

(EPRI, 1993). 

Although there should not be any hysteretic dissipation of energy at strains below the linear 

cyclic threshold shear strain, experimental results do show some energy dissipation even at 

very low shear strain levels. The breadth of the hysteresis loops increases with the increasing 

cyclic shear strain amplitude indicating that the damping ratio increases with increasing strain 

amplitudes. Figure 2-3 (Bechtel SAIC Company, LLC, 2004) shows the variation of the 

damping ratio with strain amplitude for tuff and alluvium at the Yucca Mountain site. 
2-2 

Figure 2-1. Typical Shear Stress and Shear Strain Relationship 

(Modified From Kramer, 1996) 

Figure 2-2. Variation of Normalized Shear Modulus With Shear Strain Using Data for 

Upper Mean Tuff and Upper Mean Alluvium From Bechtel SAIC Company, LLC (2004)
2-3 

Figure 2-3. Variation of Material Damping Ratio With Shear Strain Using Data for Upper 

Mean Tuff and Upper Mean Alluvium From Bechtel SAIC Company, LLC (2004) 
3-1 

3 SITE RESPONSE ANALYSIS 

1D ground response analyses assume that the medium above the bedrock is composed of 

horizontal strata and the ground response is predominantly due to the shear waves propagating 

vertically (Kramer, 1996). It is inherently assumed that the bedrock and the overlain strata are 

infinite in the horizontal direction. 

In a linear approach to ground response analysis, transfer functions are used to estimate the 

response parameters, such as displacement, velocity, acceleration, shear stress, and shear 

strain, to an input motion, such as acceleration at the bedrock (Kramer, 1996), using the 

superposition principle. The transfer function is defined as the ratio of the seismic motion 

amplitude at the soil surface to the amplitude at the bedrock. SHAKE (Schnabel, et al., 1972) is 

the well-known computer program that implement the super position approach. The known 

ground motion at the bedrock is decomposed into a Fourier series (i.e., a series of sine waves 

with different amplitudes, frequencies, and phase angles), generally using the Fast Fourier 

Transform (FFT). The contribution of each Fourier component is estimated using the transfer 

function. Using the inverse FFT, the ground motion at the top surface is estimated (Kramer, 

1996). Nonlinear behavior of the medium is generally approximated using an iterative 

procedure with equivalent linear properties (Kramer, 1996). As an alternative approach, the 

actual nonlinear response of a soil deposit can be analyzed using direct numerical integration in 

the time domain (Kramer, 1996). Currently available computer programs for nonlinear 1D 

ground response analysis represent the stress-strain behavior of the medium by different 

models (e.g., hyperbolic, modified hyperbolic, Ramberg-Osgood) (Kramer, 1996). 

The equivalent linear approach is computationally simple and provides reasonable estimate of 

soil response subjected to seismic excitation in many situations; however, it is still an 

approximation to the nonlinear seismic ground response. The shear modulus and damping ratio 

would remain constant throughout the seismic excitation, although, as discussed later in this 

report, they are strongly dependent on the amplitude of shear strain (Arslan and Siyahi, 2006). 

The difference between the results from equivalent linear and nonlinear analyses is dependent 

on the degree of nonlinearity, as shown in Figure 2-2, of the medium response, in terms of 

degradation of rigidity modulus to dynamic load. For cases with high strain levels (i.e., shear 

stress approaches the shear strength of the medium), a nonlinear analysis is likely to produce 

reasonable results. On the other hand, both types of analyses are expected to produce 

reasonable results if the level of strain is low (i.e., stiff medium and/or weak input seismic 

motion) (Kramer, 1996). The differences in estimated ground response using linear and 

nonlinear approaches follow (Kramer, 1996; Arslan and Siyahi, 2006): 

• An inherent linearity assumption may generate spurious resonance. As the stiffness of 

the medium decreases over the duration of large earthquakes, such high amplification 

levels will not develop in the field. 

• Use of an effective shear strain in an equivalent linear analysis may lead to an 

oversoftened and overdamped system in cases where the peak shear strain is 

significantly larger than the remainder of the shear strain pulse. Alternatively, in cases 

with nearby uniform shear strain amplitude, it will develop an undersoftened, 

underdamped system. 
3-2 

• Nonlinear analysis requires an appropriate constitutive model; however, the parameters 

describing such models are generally not well established compared to the equivalent 

linear model. A substantial laboratory and field testing program is needed to 

appropriately evaluate the parameters characterizing the nonlinear stress-strain model. 

1D analysis is appropriate for cases with horizontal strata or gently sloping sites with 

parallel material boundaries. If the strata are not horizontal, stiff materials are present, or 

materials with nonparallel boundaries are present, 2D and in some cases, 3D analysis will be 

required. In this study, a time domain (i.e., direct integration) method has been used to address 

the problem. 

Site response assessment is typically conducted using 1D analysis assuming horizontally 

bedded strata with linear elastic material behavior. Many recent studies, however, have 

shown that 2D and 3D analyses may be required to fully capture the complexities in the 

subsurface, such as tilted layers and sedimentary basin effects (e.g., Bardet, et al., 1997, 1992; 

Frankel, 1993; Graves, 1993; Hartzell, et al., 1996; Zhang and Papageorgiou, 1996; 

Bielak, et al., 1999, 2003; Xu, et al., 2003; Yoshimura, et al., 2003). For example, 

Bardet, et al. (1992) analyzed the site response of the Marina District of San Francisco during 

the Loma Prieta earthquake. A nonlinear 1D analysis predicted a peak acceleration of 0.12 g at 

the ground surface in the Marina District with a bedrock peak acceleration of 0.067 g. The 

predicted ground surface motion is approximately 80 percent larger than the motion at the 

bedrock. SHAKE predicted a larger amplification because it uses an equivalent linear material 

model that does not consider complete failure. The 2D analyses predicted a peak ground 

acceleration of 0.20 g, approximately three times the amplification of the bedrock motion. 

Predicted larger peak accelerations with the 2D analyses indicate that the 1D site-response 

analyses used in traditional engineering practice may not necessarily be conservative. 

Bielak, et al. (1999) compared results of 1D and 2D analyses for ground responses in an alluvial 

valley for the magnitude 6.8 earthquake in northern Armenia on December 7, 1988. 

Bielak, et al. (1999) made the following conclusions. 

• The 1D amplification ratio exhibited the resonant behavior typical of a 

flat-layered system. 

The 2D amplification ratio exhibited resonant behavior in the vicinity of the 1D resonant 

frequencies; however, the corresponding 2D frequencies were slightly higher because of 

the lateral confinement the valley provided. The values of the respective peaks were 

considerably larger in the 2D simulations. 

• The 2D amplification ratios greatly exceeded the 1D value. Interestingly, at certain 

frequencies the 2D amplification ratio was smaller than unity, indicating deamplification, 

or destructive interference of seismic waves. 
4-1 

4 ANALYSIS OF SURFACE RESPONSE 

The software FLAC Version 5.0 (Itasca Consulting Group, 2005) and UDEC Version 5.0 

(Itasca Consulting Group, 2011) have been used to simulate the waves propagating upward 

through the layered strata. FLAC is a 2D explicit finite difference program for engineering 

mechanics computation. Materials are represented by elements, or zones, which form a grid 

that the user adjusts to fit the shape of the object to be modeled. Each element behaves 

according to a prescribed linear or nonlinear stress-strain law in response to the applied 

forces or boundary restraints. The explicit, Lagrangian calculation scheme and the 

mixed-discretization zoning technique used in FLAC ensure that plastic collapse and flow are 

modeled accurately. UDEC is a 2D distinct element code that can simulate the response of a 

fractured rock mass subjected to dynamic loads. The medium is represented as an assemblage 

of discrete blocks. Large displacements are allowed along discontinuities between blocks. 

Rotation of individual blocks is allowed. 

The recorded motion of the Coalinga1 earthquake with a 5 × 10−4

 annual exceedance 

frequency, as used in Hsiung, et al. (2007), has been used in this study. This earthquake had a 

moment magnitude of 5.2 measured at the Transmitter Hill Station 10.4 km [6.5 mi] from the 

epicenter. The earthquake time history in terms of particle acceleration is shown in Figures 4-1 

and is used in FLAC simulations. Figure 4-2 shows the time history in terms of particle velocity 

used in UDEC simulations. It should be noted that the amplitude of the particle velocity time 

history has been amplified to induce fracture opening and/or failure in the UDEC simulations. 

As the UDEC simulations assumed elastic behavior of the blocks, an increase of input motion 

amplitude will not cause any rock failure and, thereby, affect the resulting ground motion at 

the surface. 

4.1 Fractured Rock Mass Model: UDEC Analysis 

This analysis assesses the effects of the fractured rock mass on the propagating seismic 

waves. The rock mass is modeled to comprise many smaller sized rock blocks. Different types 

of fractures—namely, faults, layer contact boundaries, and joints—are assumed to exist in the 

rock mass. Large deformation of the rock mass as the seismic waves pass through would 

mostly be realized by the fractures. Currently available constitutive models for intact blocks and 

joints cannot simulate the degradation of shear modulus as a function of the shear strain, as 

shown in Figures 2-1 and 2-2.1

 Therefore, relative deformation between the blocks and the 

fractures would influence the propagation of the seismic waves and an inherent assumption in 

this analysis is that the shear strain in both rock and alluvium is small. 

4.1.1 Rock Mass Model in UDEC 

Figure 4-3 shows the fractured rock mass that has been modeled in this analysis using the 

UDEC code. An area where the proposed surface facilities of the geologic repository would be 

located was selected in this analysis because significant information has been collected for this 

 

1

Development of a material constitutive model in UDEC to represent degradation of shear modulus and damping 

ratio as a function of shear strain similar to that in the FLAC code is a major endeavor and is outside the scope of 

this activity. 
4-2 

Figure 4-1. Particle Acceleration Versus Time of the Coalinga1 Earthquake Used in 

FLAC Simulations

Figure 4-2. Input Motion in Terms of Particle Velocity Versus Time Used in 

UDEC Simulations
4-3 

Figure 4-3. Cross Section of Strata Used in Seismic Response Modeling Using the UDEC Code (Desc

given in Table 4-1)
4-4 

area. Figure 4-3 shows the surface topography is not uniform and is reproduced in the UDEC 

model, as shown in Figure 4-4. The UDEC model is 600 m [1,969 ft] in length. It starts 50 m 

[164 ft] from the left side of Figure 4-3 and extends up to 650 m [2,134 ft] in the horizontal 

direction. The UDEC model extends to a depth of 200 m [656 ft] at one side and approximately 

260 m [853 ft] at the other side due to presence of a small hill. 

Additionally, in a few UDEC simulations, two joint sets are incorporated in the UDEC model in 

addition to the faults and strata layers shown in Figure 4-4. One set of joint dips toward the left 

at an angle of 30° with a spacing of 8 m [26 ft] and trace length of 30 m [98 ft]. There is a 2-m 

[6-6 ft] variation of spacing about the mean spacing value. Trace length has a variation of 10 m 

[33 ft] about the mean trace length value. Properties of faults and strata contacts are given in 

Table 4-1. Fracture properties used in different simulations are given in Table 4-2. 

Property values used in both Tables 4-2 and 4-3 are within the range of rock fracture properties 

Hsiung, et al. (1994) observed for Apache Leap tuff joints. 

Note that measured in-situ stiffness and strength properties of faults, in some cases with 

substantial thickness of gauge materials, are generally not available in the literature and are 

assumed to be within the measured properties of natural joints, both fresh (previously not 

sheared) and previously sheared with none to a small amount of gauge materials. The UDEC 

model has a zero horizontal velocity boundary condition at both lateral boundaries and a zero 

vertical velocity boundary condition at the bottom boundary. During simulation of wave 

propagation, the “free field” boundary condition is applied at both lateral boundaries. As UDEC 

requires a “quiet” boundary for application of the input seismic motion as particle velocity, the 

bottom boundary of the model is made of a “quiet” boundary. A Rayleigh damping of 5 percent 

of the critical damping, applied at the approximate natural frequency, was used in all 

simulations. The natural frequency of the model was estimated using the same model by 

suddenly applying the gravity load. The vertical velocity was measured at different points in the 

model. As the model is highly nonuniform with different types of fractures oriented at different 

directions and different types of materials with significantly different properties, the estimated 

frequency varies somewhat at different monitoring locations. An average frequency was used in 

these simulations. 

UDEC solves a dynamic problem using an explicit time-marching scheme. In UDEC simulation, 

the timestep is dependent on the stiffness properties. The timestep UDEC estimated to solve 

these complex simulations was in the neighborhood of a fraction of a microsecond. 

Consequently, runtime required for simulating the effects of the first 25 seconds of the seismic 

signal would be at least a week for the available computer resources. Consequently, it was 

decided to run up to 5 seconds for each problem to complete several simulations in a 

reasonable time. Additionally, UDEC has an option to partially scale up the mass of the blocks 

to increase the speed of calculation. A partial mass scaling was used to achieve a timestep in 

the neighborhood of 50 microseconds. 

Both frequency content and velocity of propagation of the input wave will affect the numerical 

accuracy of wave transmission. For accurate representation of wave transmission through a 

numerical model, the element (zone) size S should be 

 (4-1) 

10

λ S ≤
4-5 

Figure 4-4. The UDEC Model of the Strata Cross Section With Monitoring Points 1–9. 

Monitoring Points 1–5 are on the Surface, Points 6 and 8 at the Bottom Boundary 

(on Bedrock), and Points 7 and 9 Inside the Model.

Table 4-1. Material Properties Used in UDEC Simulations

Unit 

Description 

Unit Weight 

Shear Wave 

Velocity 

Compression 

Wave Velocity 

Shear 

Modulus 

(Pa) 

Bulk 

Modulus 

(Pa) 

Ttpul Upper 

Lithophysal 

Zone of 

Crystal-Poor 

Member of 

Topopah 

Spring Tuff 

2,200 kg/m3

[137.3 pcf] 

1,759 m/s 

[5,771 ft/s] 

3,404 m/s 

[11,165 ft/s] 

6.81 × 109 1.65 × 1010

Ttr Crystal-Rich 

Member of 

Topopah 

Spring Tuff 

2,127 kg/m3

[132.8 pcf] 

1,978 m/s 

[6,490 ft/s] 

3,228 m/s 

[10,591 ft/s] 

8.32 × 109 1.11 × 1010

Tbt2 Pre-Pah 

Canyon Tuff 

1,000 kg/m3

[62.4 pcf] 

700 m/s 

[2,300 ft/s] 

1,247 m/s 

[4,091 ft/s] 

4.90 × 108 9.02 × 108

Tpp Pah Canyon 

Tuff 

1,200 kg/m3

[74.9 pcf] 

1,102 m/s 

[3,615 ft/s] 

1,809 m/s 

[5,935 ft/s] 

1.46 × 109

 1.98 × 109

Tpy Yucca 

Mountain Tuff 

1,500 kg/m3

[93.6 pcf] 

1,230 m/s 

[4,035 ft/s] 

2,011 m/s 

[6,598 ft/s] 

2.27 × 109 3.04 × 109
4-6 

Table 4-1. Material Properties Used in UDEC Simulations (continued) 

Unit 

Description 

Unit 

Weight 

Shear Wave 

Velocity 

Compression 

Wave Velocity 

Shear 

Modulus 

(Pa) 

Bulk 

Modulus 

(Pa) 

Tcpln Lower 

Nonlithophysial 

Zone of Tiva 

Canyon Tuff 

2,127 kg/m3

[132.8 pcf] 

1,894 m/s 

[6,213 ft/s] 

3,221 m/s 

[10,568 ft/s] 

7.63 × 109 1.19 × 1010

Tcpll Lower 

Lithophysial 

Zone of Tiva 

Canyon Tuff 

2,207 kg/m3

[137.8 pcf] 

1,639 m/s 

[5,378 ft/s] 

3,107 m/s 

[10,194 ft/s] 

5.93 × 109 1.34 × 1010

Tcpmn Middle 

Nonlithophysial 

Zone of Tiva 

Canyon Tuff 

2,321 kg/m3

[144.9 pcf] 

1,583 m/s 

[5,193 ft/s] 

3,062 m/s 

[10,046 ft/s] 

5.82 × 109 1.40 × 1010

Tcpul Upper 

Lithophysial 

Zone of Tiva 

Canyon Tuff 

2,070 kg/m3

[129.2 pcf] 

1,132 m/s 

[3,713 ft/s] 

2,139 m/s 

[7,018 ft/s] 

2.65 × 109 5.94 × 109

Tcr Crystal-Rich 

Member of Tiva 

Canyon Tuff 

1,158 kg/m3

[72.3 pcf] 

1,277 m/s 

[4,190 ft/s] 

2,535 m/s 

[8,317 ft/s] 

1.89 × 109 4.92 × 109

Tpkl Tuff Unit “x” 1,572 kg/m3

[113.4 pcf] 

1,118 m/s 

[3,318 ft/s] 

2,029 m/s 

[6,657 ft/s] 

1.96 × 109 3.85 × 109

Qtac Undifferentiated 

Quaternary/Terti

ary Alluvial and 

Colluvial 

Deposit 

1,852 kg/m3

[115.6 pcf] 

712 m/s 

[2,336 ft/s] 

1,485 m/s 

[4,872 ft/s] 

9.40 × 108 2.83 × 109

Table 4-2. Properties of Faults and Strata Contacts Used in UDEC Simulations 

 Faults Strata Contacts 

Case Number 1 2 3 4 1 2 3 4 

Normal Stiffness 

kn [GPa/m] 

20 40 10 50 20 40 10 50 

Shear Stiffness ks

[GPa/m] 

10 20 5 30 10 20 5 30 

Friction Angle φ

[degree] 

30 30 10 15 30 30 20 30 

Cohesion C [MPa] 1 1 0.001 1 1 1 0.001 1 

Table 4-3. Properties of Joints Used in UDEC Simulations (Case 4) 

Left Dipping Right Dipping 

Normal Stiffness kn [GPa/m] 50 50 

Shear Stiffness ks [GPa/m] 30 30 

Friction Angle φ [degree] 30 30 

Cohesion C [MPa] 1 1 
4-7 

where λ is the wave length associated with the highest frequency component of the input wave. 

Hsiung, et al. (2007, Table 2-3) show that the maximum shear velocity occurs for the case of 

“mean + one standard deviation” of Lower Nonlithophysal Tiva Canyon rock. It is 2,260 m/s 

[7,415 ft/s]. Therefore, the wavelength associated with the wave having frequency f is 

2260/ f . For a maximum frequency of interest ݂ ൌ 50 Hz, the associated wavelength is 45 m 

[148 ft]. Consequently, the finite difference grid size in accordance with Eq. (4-1) should be less 

than or equal to 4.5 [15.6 ft] for a wave with frequency of 50 Hz to propagate adequately. The 

zone size in this model is smaller than 4.5 m [15.6 ft] for appropriate transmission of the input 

seismic waves applied at the base. 

4.1.1.1 Results and Discussion 

Locations of the monitoring points in the UDEC model are shown in Figure 4-4. Point 1 is on 

the top of the hill and is located on the rock layer. Points 2 through 5 are on the alluvium. The 

thickness of the alluvium layer is smallest at Point 2. It increases to the maximum at Point 5. 

Points 6 and 8 are at the bottom model boundary, which is assumed to be the top surface of the 

bedrock, and thereby show the seismic motion applied at the bottom boundary of the model. 

Points 7 and 9 are a small distance away: 20 m [66 ft] and 10 m [33 ft], respectively. 

Models Fractured Rock Mass Only

Case 1: Fault Properties (kn = 20 GPa/m, ks = 10 GPa/m, φ = 30°, C = 1 MPa) 

 Contact Properties (kn = 20 GPa/m, ks = 10 GPa/m, φ = 30°, C = 1 MPa)

This case represents a typical scenario expected in nature. Both types of fractures have 

average stiffness and frictional strength properties and are well within the range observed 

(e.g., Hsiung, et al., 1994). It has been assumed that these fractures do not have any tensile 

strength, a generally conservative assumption. Figure 4-5 compares the particle velocity 

estimated at Point 1 with respect to the input motion for an analysis conducted for the first 

5 seconds of the earthquake motion signal. The seismic motion after traveling through different 

layers is amplified as it passes through Point 1; the amplification is approximately 1.4, and peak 

amplitude is 0.43 m/s [1.41 ft/s] compared to 0.31 m/s [1.02 ft/s] of the input motion. 

Additionally, the shift observed in the time-varying signal is due to the finite propagation velocity 

of the seismic motion through different layers. Amplification of the seismic motion increases as 

the soil thickness increases. The soil thickness is maximum at Point 5, and the peak amplitude 

is 0.45 m/s [1.48 ft/s], an amplification approximately 1.5 times the input motion and 

approximately 7 percent higher than at Point 1. Figure 4-6 compares the amplification of 

seismic motion at Points 1 through 5 relative to the input motion at Point 6. Figures 4-5 and 4-6 

show that individual pulses of the input seismic motion did not appreciably change shape as the 

waves propagate upward, when an elastic model without any damage parameter to represent 

the shear modulus and damping ratio degradation as a function of the shear strain was used in 

this analysis. Consequently, the seismic amplification estimated at the surface is smaller than in 

reality, where rigidity modulus degradation of the material may enhance the amplification of the 

seismic motion at the surface. UDEC results also show minor opening {approximately 0.6 mm 

[0.24 in]} of a strata contact near the surface. These results are compared with those from other 

cases in Table 4-4. 
4-8 

Figure 4-5. Comparison of Ground Motion at Point 1 With Input Motion (Point 6) 

(Case 1) 

Figure 4-6. Comparison of Ground Motion at Points 1–5 With Input (Point 6) 

(Case 2)
4-9 

Table 4-4. Results From Four Simulations Using UDEC 

Case Amplification Comments 

 Peak Value Percentage of Input 

1 0.43 m/s [1.41 ft/s] 140 Amplification is approximately 7 percent 

higher at Point 5, which is over the 

deepest part of the soil layer. 

2 0.35 m/s [1.15 ft/s] 113 Significantly stiffer fracture planes 

resulting in an overall stiffener rock mass. 

3 0.49 m/s [1.61 ft/s] 160 Significantly softer fracture planes 

resulting in an overall softer rock mass 

4 0.49 m/s [1.61 ft/s] 160 Rock mass has two sets of joints; 

however, the fracture planes and joints 

are significantly stiffer. 

Case 2: Fault Properties (kn = 40 GPa/m, ks = 20 GPa/m, φ = 30°, C = 1 MPa) 

Contact Properties (kn = 40 GPa/m, ks = 20 GPa/m, φ = 30°, C = 1 MPa)

This case also may be a typical one in nature. Both normal and shear stiffness values are well 

within the range observed (e.g., Hsiung, et al., 1994). As the stiffness properties of both faults 

and strata contact planes increase, the timesteps required to complete 5 seconds of the 

Coalinga1 earthquake motion increases. Figure 4-7 compares the estimated seismic motion at 

the surface (Points 1–5) when the seismic motion is applied at the model base. 

The peak amplitude of the seismic motion at Point 1 is approximately 0.35 m/s [1.15 ft/s], an 

amplification of 1.13 times of the input motion. The amplitude at Point 5 is 0.39 m/s [1.28 ft/s], 

an amplification of approximately 1.30 times. Reduction of the amplification magnitude 

compared to Case 1 can be attributed to the stiffer fracture planes, which make the rock mass 

stiffer overall. The time shift of the pulse at the top surface compared to the pulse applied at the 

base is smaller compared to that observed in Case 1, which has a less stiff rock mass. Again, 

the seismic pulses do not undergo any significant distortions, because an elastic model with no 

material degradation parameter has been assumed in the constitutive model for the rock blocks 

in this case. Results are provided in Table 4-4 for comparing with other cases. 

Case 3: Fault Properties (kn = 10 GPa/m; ks = 5 GPa/m, φ = 10°, C = 0.001 MPa) 

Contact Properties (kn = 5 GPa/m, ks = 5 GPa/m, φ = 20°, C = 0.001 MPa 

This case may represent a scenario with low-end stiffness and strength values for both faults 

and strata contact planes; both faults and strata contact planes are assumed to be weak. 

Contact planes are assumed somewhat stronger than the fault planes. The amplification of the 

seismic motion at all points increases to approximately 1.6 times the input motion and the 

estimated amplitude is 0.49 m/s [1.61 ft/s], as shown in Figure 4-8. However, the difference 

between amplifications at Points 1 and 5 is small, less than 5 percent. Low stiffness of the faults 

and contact planes made the medium relatively soft, and, consequently, it can undergo larger 

deformation. Figure 4-8 also shows that transmission of seismic waves did not introduce 

distortion to the pulse shape; the materials did not degrade because elastic behavior of the 

materials was assumed. Information in Table 4-4 compares the results with those from other 

cases analyzed. 

 
4-10 

Figure 4-7. Comparison of Ground Motion at Points 1 and 5 with Input (Point 6) 

(Case 2)

Figure 4-8. Comparison of Ground Motion at Points 1 and 5 With Input (Point 6) 

(Case 3) 
4-11 

Case 4: Fault Properties (kn = 50 GPa/m; ks = 50 GPa/m, φ = 15°, C = 1 MPa) 

Contact Properties (kn = 50 GPa/m, ks = 30 GPa/m, φ = 15°, C = 1 MPa) 

Joint Properties (kn = 50 GPa/m, ks = 30 GPA/m, φ = 30°, C = 1 MPa) 

In this case, two sets of joints were incorporated into the rock mass model of faults and strata 

contact planes, as shown in Figure 4-9. One set of joints dips toward the left, and the other 

set dips toward the right. These additional fractures split each block into several small-sized 

blocks and make the model complex. The number of intersecting fractures and resulting blocks 

increases several fold. Tracking movement of these blocks in each calculation timestep adds 

significant bookkeeping tasks in UDEC runs. Consequently, this simulation model in UDEC will 

run significantly slower. The runtime for the input motion (Coalinga1 earthquake) up to 

5 seconds took several days. Consequently, only one case could be successfully analyzed. 

Adding joints into the rock mass makes the overall material less stiff, which decreased the 

estimated average natural frequency to approximately 15 Hz. Figure 4-10 shows the estimated 

wave motion at Point 5 (on soil) compared to the input motion. The peak amplitude at Point 5 is 

0.49 m/s [1.61 ft/s] compared to an input of 0.31 m/s [1.02 ft/s], approximately 1.6 times 

amplification. All measuring points on the surface (Points 1–5) recorded higher magnification of 

the seismic waves than the cases with no joints. As discussed, joints made the overall stiffness 

of the medium softer resulting in more deformation (as well as larger particle velocity) at the 

surface. Results are presented in Table 4-4. 

4.2 Horizontal and Inclined Strata Models: FLAC Analysis 

Two models were developed to assess the effects of inclined strata on the ground response to 

be realized on the surface as a result of seismic excitation at the base. The first model is 

composed of horizontal strata, and the second model has the same material layers dipping to 

the right. The difference in ground response between these two models would be due to 

inclined strata: the longer path to be traversed by the vertically propagating waves. The models 

are shown in Figures 4-11 and 4-12. 

4.2.1 Strata Model 

In this analysis, the material properties of the strata are assumed to be those of rock and 

alluvium at the Yucca Mountain site because the material properties data, especially the 

reduction of shear modulus with increasing shear strain, are readily available for this region 

(Bechtel SAIC Company, LLC, 2004). However, note that the analysis assumes a somewhat 

arbitrary thickness for each layer, which does not necessarily conform to strata thicknesses 

at the Yucca Mountain site. The material properties assumed for each stratum are given in 

Table 4-1. The grid size in the FLAC models was smaller than necessary to transmit waves 

with a frequency of 20 Hz. 

The free-field boundary condition, available in FLAC, was used to simulate the infinite extent of 

the rock mass in the horizontal direction. This boundary condition accounts for the free-field 

motion that would exist in the subsurface strata in the absence of the “artificial” model 

boundaries as the input seismic waves propagate upward. This boundary condition not only 

allows placing the model boundaries sufficiently close resulting in a small-sized model, but also 

makes the boundaries nonreflecting. Plane waves, while propagating upward, suffer no 

distortion at the boundary as the free-field boundary; they are modeled as a 1D grid attached to 

the main grid representing the medium, supplying identical conditions of a model having infinite 

extent. As this boundary condition requires the model to be in static equilibrium, the model is 
4-12 

Figure 4-9. UDEC Model of the Fractured and Jointed Rock Mass; One Fracture Set Dips 

Toward the Left and Another Toward the Right 
4-13 

Figure 4-10. Comparison of Ground Motion at Point 5 (on Soil) With Input (Point 6) 

(Case 4) 

Figure 4-11. Horizontal Strata Model Used in FLAC Analysis 
4-14 

Figure 4-12. Inclined Strata Model Used in FLAC Analysis 

first cycled 20,000 calculation steps with only gravity load to keep the unbalanced force in the 

model negligible. This state of model was used as the initial condition for propagating the 

seismic motion. The seismic motion of the Coalinga1 earthquake is applied at the bottom 

boundary as a table of discrete values of particle acceleration. 

4.2.2 Model With Horizontal Strata Layers 

In these simulations, the strata used are composed of horizontal layers, as shown in 

Figure 4-11. A grid of 50 horizontal and 200 vertical units was used to represent a block of 

materials 250 × 250 m [820 × 820 ft]. The materials deform as the seismic waves pass through 

and sustain damage in terms of degradation of rigidity modulus, as shown in Figures 2-2 

and 2-3. 

Reduction of normalized shear modulus and material damping ratio with shear strain has been 

represented by three different models available in FLAC and used in this analysis: (i) default, 

(ii) Sig3, and (iii) Sig4. Information for both tuff (Upper Mean Tuff) and alluvium (Upper Mean 

Alluvium) from Bechtel SAIC Company, LLC (2004) was fitted to these models. Bechtel SAIC 

Company, LLC (2004) combined results of laboratory tests of all tuff rocks as one material. In 

this study, the same degradation model parameters were used for all tuff layers. Results of 

fitting are discussed in the appendix. The fitted parameters are given in Tables 4-5 through 4-7. 

Default Material Degradation Model

Figure 4-13 shows the estimated velocity history at the top surface. As this is a horizontal strata 

model, there is no spatial variation of the velocity (as well as acceleration) history at the base 

and the top surface. Figure 4-14 compares the input acceleration at the base and the top 

surface. The acceleration input at the base with peak amplitude of 0.37 m/s2

, as shown in 
4-15 

Table 4-5. Fitted Parameters for Default Model Representing Material Degradation With 

Shear Strain 

Tuff Alluvium 

L1 −3.425 −3.425 

L2 3.23 0.823 

Table 4-6. Fitted Parameters for Sig3 Model Representing Material Degradation With 

Shear Strain

 Tuff Alluvium 

a 0.997 1.012 

b −0.3034 −0.4753 

x0 −0.2334 −1.1418 

Table 4-7. Fitted Parameters for Sig4 Model Representing Material Degradation With 

Shear Strain 

 Tuff Alluvium 

a 341.99 1.0217 

b −0.5125 −0.4837 

x0 3.227 −1.1322 

y0 −340.9961 −8.08 × 10−3

Figure 4-1, amplifies by approximately 3.4 m/s2

 peak acceleration at the top of material column, 

almost 9 times amplification. Note in Figure 4-14 that the first major pulse of the input motion 

corresponds with a relatively large amplitude response followed by several large amplitude 

responses with similar amplitudes. It seems the first major pulse of the input motion 

degrades the rigidity modulus of the material significantly. This degradation is permanent. 

Input pulses that follow this first major pulse generate a significantly large response of the 

degraded materials. 

The vibration at the top continues for a long time even after the input motion has died down. 

This is expected as the materials will continue to vibrate for some time. Note that the simulation 

only captured the first 20 seconds of the input motion. As a result, the trailing end of the 

surfaceground motion was not calculated. Although pseudo-acceleration response spectra, 

shown in Figure 4-15, have a peak around 1.2 Hz, other peaks occur at 3, 5, and 8 Hz. 

Sig3 Material Degradation Model

Figure 4-16 compares the estimated horizontal acceleration at the surface with the input motion 

at the base. Input acceleration at the base is amplified to 3.4 m/s2

 peak acceleration, 

approximately nine times amplification, as the waves propagate upward. Like the previous case 

with the default model simulating degradation of material rigidity this model also shows peak 

pulse at the surface occurring at much later than the peak pulse of the input motion. The top 

surface seems to vibrate for a long time after cessation of a ground acceleration pulse 

containing significant energy. 
4-16 

 

Figure 4-13. Estimated Velocity History at the Top of the Model for First 18 Seconds With 

the Default Material Degradation Model 

[y-axis: Velocity (m/s), x-axis: Time (s)]

Figure 4-14. Comparison of Acceleration at the Top of the Model With Input Motion at the 

Base With the Default Material Degradation Model 

[y-axis: Acceleration (m/s2

), x-axis: Time (s)]
4-17 

Figure 4-15. Spectral Acceleration at the Top Surface of the Horizontal Layer Model as a 

Function of Frequency 

Figure 4-16. Comparison of Acceleration at the Top Surface with Input Motion at the 

Base With Sig3 Material Degradation Model 

[y-axis: Acceleration (m/s2

), x-axis: Time (s)]
4-18 

Sig4 Material Degradation Model

Figure 4-17 compares the estimated horizontal acceleration at the surface of the material 

column with the input ground acceleration at the base. The input acceleration is significantly 

amplified at the top the estimated peak acceleration is almost nine times the input motion. The 

amplification is extremely high even at 18 seconds after the start of the seismic motion. 

4.2.3 Model With Inclined Strata Layers 

In this simulation, the strata layers dip toward the right (e.g., a limb of a fold structure), as shown 

in Figure 4-12. Material degradation increases with increase in shear strain, as shown in 

Figures 2-2 and 2-3. It is expected that the spatial variation of acceleration measured at the 

surface would be due to different materials encountered at different locations. Reduction of 

normalized shear modulus and material damping ratio with shear strain has been represented 

by the default model only in this analysis. 

Figure 4-18 shows the variation of estimated horizontal acceleration at different locations on the 

surface. The input acceleration is amplified to an estimated peak acceleration of 3.5 m/s2

, and 

amplification of approximately 9.5 times the input motion. The top surface vibrates for a 

considerably longer time. Although not markedly different in these simulation result shown in 

Figure 4-18, the amplification increases to the right side of the model, where the thickness of the 

alluvium top layer increases continuously. 

4.2.4 Basin Model: FLAC Analysis 

As Bielak, et al. (1999) observed, the impedance contrast between the soil and rock layers, 

especially in a valley or basin-type deposit, controls the amplification of the input ground motion 

as the waves propagate to the surface. A higher contrast ratio produces larger amplification. In 

this analysis, a simple basin model has been developed to analyze the phenomena. The model 

is shown in Figure 4-19: a basin of relatively softer material (e.g., a soil deposit) is surrounded 

by a stiffer material (e.g., a rock deposit). All three rock types have the same material 

properties and, therefore, belong to the same rock unit. 

The model is 500 m [1,640 ft] long and has a depth of 250 m [820 ft]. The soil basin is in the 

middle and has a length of 300 m [984 ft] with a depth of 150 m [492 ft]. The material properties 

used are given in Table 4-8 (QTac for soil and Tcpln for rock). Hysteretic damping of both soil 

and rock is represented by the default model available in FLAC with parameters, as given in 

Table 4-4. The impedance ratios between rock and soil used in these simulations are given in 

Table 4-9. 

Figures 4-20 and 4-21 show the amplification of the input seismic motion in the large soil basin 

with an impedance ratio (ZRock/ZSoil) of 5.0. Amplification is not uniform across the soil basin. 

Amplification increases toward the center of the basin and reaches maximum at the center. 

This phenomenon is present in all three cases analyzed here. Amplification also increases as 

the impedance ratio increases between the rock and soil, as can be seen in Figures 4-22 

through 4-24 for three cases with increasing impedance ratio. 
4-19 

Figure 4-17. Comparison of Acceleration at the Top Surface With Input Motion at the 

Base With Sig4 Material Degradation Model 

[y-axis: Acceleration (m/s2

), x-axis: Time (s)] 

Figure 4-18. Comparison of Acceleration Time Histories at Left and Right Corners of the 

Inclined Strata Model Shown in Figure 4-12 

[y-axis: Acceleration (m/s2

), x-axis: Time (s)] 
4-20 

Figure 4-19. Basin Model Used With Monitoring Points 

(Note All Three Rock Blocks Have The Same Properties) 

Table 4-8. Material Properties Used in Simulations 

Case Soil Rock 

Density 

(kg/m3

) 

Bulk 

Modulus 

(GPa) 

Shear 

Modulus 

(GPa) 

Density 

(kg/m3) 

Bulk 

Modulus 

(GPa) 

Shear 

Modulus 

(GPa) 

1 1,852 2.83 0.94 2,127 11.90 7.63 

2 1,852 2.83 0.94 2,321 14.00 5.82 

3 1,852 2.83 0.94 2,200 26.40 19.78 

Table 4-9. Impedance Ratio Between Rock and Soil Used in Simulations 

Case Impedance Ratio (ZRock/ZSoil) 

1 2.79 

2 3.05 

3 5.00 
4-21 

Figure 4-20. Comparison of Seismic Wave Amplification Between Rock and Soil at the 

Center of the Basin (Impedance Ratio 5.0) 

[y-axis: Acceleration (m/s2

), x-axis: Time (s)] 

Figure 4-21. Comparison of Seismic Wave Amplification Between Rock and Soil at the 

Edge of the Basin (Impedance Ratio 5.0) 

[y-axis: Acceleration (m/s2

), x-axis: Time (s)] 
4-22 

Figure 4-22. Comparison of Seismic Wave Amplification Across the Soil Basin 

(Seismic Impedance Ratio 2.79) 

[y-axis: Acceleration (m/s2

), x-axis: Time (s)] 

Figure 4-23. Comparison of Seismic Wave Amplification Across the Soil Basin 

(Seismic Impedance Ratio 3.05) 

[y-axis: Acceleration (m/s2

), x-axis: Time (s)] 
4-23 

Impedance contrast between the soil basin and the surrounding rock mass allows a small 

portion of the input seismic energy to dissipate in the rock mass. A large part of the input 

waves, once these propagate inside the soil basin, gets trapped by repeated reflections from the 

soil–rock interface. Only a small portion of the energy will be able to be transmitted back into 

the rock. This amplification in the soil basin will be more severe as the impedance contrast 

between the rock and soil increases. As the soil is significantly weak, it undergoes severe 

deterioration of the shear modulus. The soil mass may transform practically into an almost 

liquefied state with small rigidity to resist the seismic motion, especially with Case 3 having the 

impedance ratio 5.0. As a result, the soil mass vibrates significantly under the load from the 

trapped seismic waves. It is expected that the oscillation will continue for a significantly 

longer duration. 

Figures 4-25 and 4-26 show the spectral acceleration plots across the soil basin with impedance 

ratios 3.05 and 5.0, respectively. These figures show that the middle of the basin vibrates more 

in resonance than the regions near the soil–rock interface. 

Figure 4-24. Comparison of Seismic Wave Amplification Across the Soil Basin 

(Seismic Impedance Ratio 5.0) 

[y-axis: Acceleration (m/s2

), x-axis: Time (s)] 
4-24 

Figure 4-25. Comparison of Spectral Acceleration Across the Soil Basin 

(Seismic Impedance Ratio 3.05) 

Figure 4-26. Comparison of Spectral Acceleration Across the Soil Basin 

(Seismic Impedance Ratio 5.0) 
5-1 

5 SUMMARY AND CONCLUSIONS 

In this study, 2D analyses have been conducted to assess the ground motion to be expected at 

the top surface of a layered, faulted, and fractured deposit subjected to an earthquake. The 

rock mass in this model is overlain by a thick soil layer. The rock mass is composed of 

several different layers of rock with widely varying properties, including shear wave velocity. 

The analyses used the reference geology and materials properties data, particularly the 

measured values of degraded shear modulus with increasing shear strain, from the proposed 

geologic repository site at Yucca Mountain, Nevada (Bechtel SAIC Company, LLC, 2004). . The 

recorded motion from the Coalinga1 earthquake (Hsiung, et al., 2007) was applied at the bottom 

of the modeled region, generally at a depth of 200 to 260 m [656 to 853 ft] from the ground 

surface. Simulations were conducted using two widely used computer programs in soil and rock 

engineering UDEC (Itasca Consulting Group, 2011) and FLAC (Itasca Consulting Group, 2005). 

This report documents the results obtained from the analyses. The analyses results are 

preliminary and should not be used in a licensing proceeding until further investigation is 

conducted on the estimated ground response. 

UDEC has the capability to model numerous intersecting fractures easily. Consequently, UDEC 

was used to assess the effects of a faulted and jointed rock mass on wave propagation. UDEC 

is a discontinuum code: a major portion of the rock mass deformation is realized by fracture 

deformation and slippage, which moves the intact rock blocks. Currently, the UDEC code does 

not have a constitutive model for the block material that can simulate the nonlinear shear 

modulus degradation as a function of shear strain, as shown in Figures 2-2 and 2-3 In this 

analysis, only linear elastic behavior of the intact rock blocks was used. Consequently, these 

analyses do not capture the rock shear strain effects on the seismic waves as they pass through 

a rock or soil layer. 

Results of UDEC analysis show the following: 

(1) Amplification of the input seismic waves increases with increase in soil thickness above 

the rock mass as the soil has a lower modulus of rigidity and, hence, slower shear wave 

velocity than the rock. 

(2) Stiffness properties of the fractures (both joints and faults) in the rock mass also 

influence the seismic motion at the ground surface; lower stiffness properties make the 

overall material less stiff, and as a result, the material amplifies the seismic motion 

more. The amplification observed in this study was generally less than 1.6 times the 

original peak motion. 

(3) Because UDEC can only simulate linear elastic deformation of the blocks, the estimated 

amplification of the seismic waves at the ground surface is low, compared to the wave 

amplification obtained from FLAC analysis. 

To assess the effects of shear modulus reduction of the rock and soil layers as a function of 

induced shear strain on the surface ground motion, the computer code FLAC was used. Three 

built-in models of FLAC were used to represent the modulus degradation as the waves pass 

through. Note that FLAC is a continuum code, unlike UDEC. Consequently, it cannot simulate 

the presence of multiple intersecting fractures. Effects due to both horizontal and inclined strata 
5-2 

layers were assessed using the FLAC code. In addition, a soil basin model was analyzed to 

assess the effects of impedance mismatch between the soil and the surrounding rock medium. 

Results of FLAC analysis show the following: 

(1) Response at the ground surface is affected by material layers with different thickness 

and material properties, similar to that in UDEC analysis. Response at the ground 

surface is also affected by the shear modulus reduction with shear strain as a result 

of the hysteretic damping of each layer at the site. An amplification of approximately 

9.5 times was observed for the geometry of soil and rock layers including their 

material properties used in this analysis. This amplification factor is expected to 

be different for other material geometries with different properties because these are 

site-specific parameters. 

(2) In the inclined strata model, the length of the path to be traversed by the vertically 

propagating wave through each material layer increases. This increased path length 

renders additional effects on the surface ground motion, which is expected to be 

proportional to the extra path length, other parameters remaining the same. 

(3) In the basin model, amplification of the ground motion is largely controlled by 

the impedance contrast between the soil and the surrounding rock. As the impedance 

contrast increases (higher ZRock/ZSoil ratio), the seismic wave amplitude inside the soil 

basin increases significantly. Moreover, the wave amplification varies spatially. The 

amplification is lower near the edges and increases toward the middle of the basin. 

(4) A soil basin with high impedance contrast with the surrounding medium appears to trap 

seismic wave energy within the basin. If the seismic energy is sufficient (large shear 

strain) to reduce the modulus significantly, the response motion could have a large 

magnitude and significantly long duration. 

(5) Strong ground motion may significantly affect the response at the ground surface. It 

seems that the shear modulus of rock and soil reduces significantly and the damping 

ratio increases, following the curves shown in Figures 2-2 and 2-3, when the first major 

pulse of the input ground motion arrives. Alluvium (soil) undergoes more shear modulus 

reduction than rock (Figures 2-2 and 2-3). Both materials thereafter respond to the 

incoming waves by developing oscillatory pulses with large amplitude; amplification 

approximately 9.5 times has been observed in this study. Bielak, et al. (1999) observed 

eight times amplification of the seismic motion over thick soil deposits in Armenia from 

the 1988 Armenia earthquake. Kramer (1996) indicated alluvial valleys can generate a 

very complex, even chaotic, motion. 

(6) Estimated acceleration histories using the layered media and the basin model 

(Figures 4-14, 4-16, 4-17, 4-18, 4-20–4-24) indicate resonant behavior. Plots of spectral 

acceleration (Figures 4-25 and 4-26) also show concentration of energy at different 

frequencies, indicating resonant behavior. Moreover, the frequencies at which resonant 

behavior was observed vary with the material properties. The resonant frequencies are 

higher with larger seismic impedance ratio, as can be seen in Figures 4-25 and 4-26. 

Similar resonant response has been observed in measured and analyzed responses 

over thick soil deposits in other studies [e.g., Bielak, et al. (1999, Figure 6); Kramer 

(1996, Figure 8.15)]. Bielak, et al. (1999) observed vibration of soil at four resonant 

frequencies near the deepest part of the soil deposit. The stiffer region vibrated at only 

one resonant frequency. 
5-3 

This study did not conduct 1D simulations. As a result, quantitative differences between 1D 

and 2D analyses could not be assessed; however, the trend observed in the estimated 

ground motion in the 2D analyses is similar to that observed in other studies in which 1D 

and 2D simulations were directly compared (e.g., Bielak, et al., 1999). Bielak, et al. (1999) 

observed higher amplification across the valley with deep soil deposits. Similar 

observations can be made for the basin model in this analysis (Section 4.2.4). 

Bielak, et al. (1999) noted that the spatial variation of site response cannot be reproduced 

using 1D analysis. In addition, Bielak, et al. (1999) observed resonant behavior in 2D 

simulations at slightly higher frequencies than that observed in 1D simulation. They 

attributed this behavior to lateral confinement provided by the valley. Additionally, they 

observed resonant behavior at frequencies not observed in 1D simulation. 
6-1 

6 REFERENCES 

Arslan, H. and B. Siyahi. “A Comparative Study of Linear and Nonlinear Site Response 

Analysis.” Environmental Geology. Vol. 50. pp. 1,193–1,200. 2006. 

Bardet, J.P., I.M. Idriss, T.D. O’Rourke, N. Adachi, M. Hamada, and K. Ishihara. “North America 

Japan Workshop on the Geotechnical Aspects of the Kobe, Loma Prieta, and Northridge 

Earthquakes.” Report to National Science Foundation, Air Force Office of Scientific Research, 

and Japanese Geotechnical Society. Kobe, Japan: National Science Foundation. 1997. 

Bardet, J.P., M. Kapuskar, G.R. Martin, and J. Proubet. “Site Response of the Marina District of 

San Francisco During the Loma Prieta Earthquake.” The Loma Prieta, California, Earthquake of 

October 17, 1989; Marina District; Strong Ground Motion and Ground Failure. T.D. O’Rourke 

and T.L. Holzer, eds. U.S. Geological Survey Professional Paper 1551-F. Denver, Colorado: 

U.S. Geological Survey. pp. 85-140. 1992. 

Bechtel SAlC Company, LLC. “Development of Earthquake Ground Motion Input for Preclosure 

Seismic Design and Postclosure Performance Assessment of Geologic Repository at 

Yucca Mountain, Nevada.” MDL–MGR–GS–000003. Rev. 01. Las Vegas, Nevada: 

Bechtel SAlC Company, LLC. 2004. 

Bielak, J., J. Xu, and O. Ghattas. “Earthquake Ground Motion and Structural Response in 

Alluvial Valleys.” Journal of Geotechnical and Geoenvironmental Engineering. Vol. 125. 

pp. 413-423. 1999. 

Bielak, J., K. Loukakis, Y. Hisada, and C. Yoshimura. “Domain Reduction Method for 

Three-Dimensional Earthquake Modeling in Localized Regions, Part I: Theory.” Bulletin of the 

Seismological Society of America. Vol. 93. pp. 817-824. 2003. 

EPRI. “Guidelines for Determining Design Basis Ground Motion.” Palo Alto, California: 

Electric Power Research Institute. Vol. 1, Appendix 7.A. 1993. 

Frankel, A. “Three-Dimensional Simulations of Ground Motions in the San Bernardino Valley, 

California, for Hypothetical Earthquakes on the San Andreas Fault.” Bulletin of the 

Seismological Society of America. Vol. 83. pp. 1,020-1,041. 1993. 

Graves, R.W. “Modeling Three-Dimensional Site Response Effects in the Marina District Basin, 

San Francisco, California.” Bulletin of the Seismological Society of America. Vol. 83, No. 4. 

pp. 1,042-1,063. 1993. 

Hartzell, S., A. Leeds, A. Frankel, and J. Michael. “Site Response for Urban Los Angeles Using 

Aftershocks of the Northridge Earthquake.” Bulletin of the Seismological Society of America. 

Vol. 86. pp. 168-192. 1996. 

Hsiung, S.M., A. Ghosh, and A.H. Chowdhury. “Site Response and Soil–Structure Interaction 

Analyses—A Progress Report.” San Antonio, Texas: Center for Nuclear Waste Regulatory 

Analyses. 2007. 
6-2 

Hsiung, S.M., D.D. Kana, M.P. Ahola, A.H. Chowdhury, and A. Ghosh. NUREG/CR–6178, 

“Laboratory Characterization of Rock Joints.” San Antonio, Texas: Center for Nuclear Waste 

Regulatory Analyses. 1994. 

Itasca Consulting Group, Inc. “UDEC Version 5.0 User’s Guide.” Minneapolis, Minnesota: 

Itasca Consulting Group, Inc. 2011. 

Itasca Consulting Group, Inc. “FLAC Fast Lagrangian Analysis of Continua User’s Guide 

(Version 5.0).” Minneapolis, Minnesota: Itasca Consulting Group, Inc. 2005. 

Kramer, S.L. Geotechnical Earthquake Engineering. Upper Saddle River, New Jersey: 

Prentice Hall. 1996. 

Schnabel, P.B., J. Lysmer, and H.B. Seed. “SHAKE: A Computer Program for Earthquake 

Response Analysis of Horizontally Layered Sites.” EERC72–12. Berkeley, California: 

University of California. 1972. 

Xu, J., J. Bielak, O. Ghattas, and J. Wang. “Three-Dimensional Nonlinear Seismic Ground 

Motion Modeling in Basins.” Physics of the Earth and Planetary Interiors. Vol. 137. pp. 81-95. 

2003. 

Yoshimura, C., J. Bielak, Y. Hisada, and A. Fernandez. “Domain Reduction Method for 

Three-Dimensional Earthquake Modeling in Localized Regions, Part II: Verification and 

Applications.” Bulletin of the Seismological Society of America. Vol. 93. pp. 825-840. 2003. 

Zhang, B. and A.S. Papageorgiou. “Simulation of the Response of the Marina District Basin, 

San Francisco, California, to the 1989 Loma Prieta Earthquake.” Bulletin of the Seismological 

Society of America. Vol. 86. pp. 1,382-1,400. 1996. 
APPENDIX 

HYSTERETIC DAMPING
A–1 

APPENDIX—ESTIMATION OF CONSTANTS FOR HYSTERETIC 

DAMPING MODELS 

The equivalent linear method (e.g., SHAKE analysis) assumes linearity in the solution process. 

Therefore, strain-dependent modulus and damping functions are taken in an average sense. To 

approximate the nonlinear effects of material softening and damping, as shown in Figures 2-1 

and 2-2, FLAC can model hysteretic damping, which does not drastically reduce the timesteps. 

In the FLAC formulation, the mean strain rate tensor is calculated before any calls to the 

constitutive model. The constitutive model then uses the normalized tangent modulus as the 

multiplier for adjusting the apparent tangent shear modulus of the full zone being processed. 

In earthquake engineering, the modulus degradation curves are usually available as tables of 

values with cyclic strains spaced logarithmically. Because a derivative of this curve is 

necessary to represent the normalized tangent modulus, the coarse spacing of the tabular 

values may lead to unacceptable errors when derivatives are calculated numerically. 

Consequently, analytical models should be fitted to the tabular data before simulations using 

FLAC can be conducted. In this analysis, three models of FLAC (Itasca Consulting Group, Inc., 

2005) have been used to represent the hysteretic damping of materials, as described next. 

The normalized shear modulus and material damping with respect to shear strain used in 

this analysis are given in Tables A–1 and A–2 for tuff (Tmbt1 through Tpcyln) and alluvium 

(Qal), respectively. 

Data presented in Tables A–1 and A–2 are plotted in Figures A–1 and A–2, respectively. 

Default Model

The modulus versus logarithm of cyclic strain curve is S shaped and can be represented by a 

cubic equation with zero slope at both ends (at both low and high strains). The secant modulus 

Ms is 

M S= s

2

(3− 2s)

where 

s= L2− L

L2− L1

and L is the logarithmic strain 

L= log10 (γ) , 

where L1 and L2 are the strain values with zero slopes. Thus, putting L1 = 3 and L2 = 1 extends 

that the S curve from a low cyclic strain of 103

 or 0.001 percent to a high cyclic strain of 101 or 

10 percent. 
A–2 

Sigmoidal Models Sig3 and Sig4

Sigmoidal curves can be used to represent modulus degradation curves. Two models with 3 

(Sig3) and 4 (Sig4) parameters are 

M s= a

1+ exp [− (L− x0)/b]

and 

M s= y0+ a

1+ exp [− (L− x0 /b)]

These three models were fitted to the data presented in Tables A–1 and A–2. The fitted 

parameters are given next: 

Default Model

Tuff Alluvium 

L1 3.425 L1 3.425 

L2 3.23 L2 0.823 

Sig3 Model Parameters

Tuff Alluvium 

a 0.997 a 1.012 

b −0.3034 b −0.4753 

x0 −0.2334 x0 −1.1418 

Sig4 Model Parameters

Tuff Alluvium 

a 341.99 a 1.0217 

b −0.5125 b −0.4837 

x0 3.2271 x0 −1.1322 

y0 −340.9961 y0 −8.07 × 10−3

Figures A–3 and A–4 show the fit of these parameters to the data obtained at the 

Yucca Mountain site for “Upper Mean Tuff” and “Upper Mean Alluvium” (Bechtel SAIC, 

Company, LLC, 2004). 

REFERENCES 

Bechtel SAlC Company, LLC. “Development of Earthquake Ground Motion Input for Preclosure 

Seismic Design and Postclosure Performance Assessment of Geologic Repository at 

Yucca Mountain, Nevada.” MDL–MGR–GS–000003. Rev. 01. Las Vegas, Nevada: 

Bechtel SAlC Company, LLC. 2004. 
A–3 

Itasca Consulting Group, Inc. “FLAC Fast Lagrangian Analysis of Continua User’s Guide 

(Version 5.0).” Minneapolis, Minnesota: Itasca Consulting Group, Inc. 2005. 

Table A–1. Measured Damping Data for Tuff (Tmbt1 Through Tpcyln)* 

Average Shear Strain 

(Percent) Normalized Shear Modulus Material Damping Ratio 

1.0000 × 10−4 1.0000 0.5000 

3.0000 × 10−4 1.0000 0.5000 

1.0000 × 10−3 1.0000 0.5000 

3.2000 × 10−3 1.0000 0.5000 

0.0100 0.9900 0.5000 

0.0316 0.9700 0.5000 

0.1000 0.8700 2.5400 

0.3162 0.7000 5.6300 

1.0000 0.4600 10.5900 

1.3162 0.0800 15.0000 

*Bechtel SAlC Company, LLC. “Development of Earthquake Ground Motion Input for Preclosure Seismic Design and 

Postclosure Performance Assessment of Geologic Repository at Yucca Mountain, Nevada.” MDL–MGR–GS–

000003. Rev. 01. Las Vegas, Nevada: Bechtel SAlC Company, LLC. 2004.

Table A–2. Measured Damping Data for Alluvium (Qal)* 

Average Shear Strain 

(Percent) Normalized Shear Modulus Material Damping Ratio 

1.0000 × 10−4 1.0000 0.7000 

3.0000 × 10−4 1.0000 0.8000 

1.0000 × 10−3 1.0000 0.8000 

3.2000 × 10−3 0.9700 1.1200 

0.0100 0.8700 1.8000 

0.0316 0.6800 3.5300 

0.1000 0.4300 7.1000 

0.3162 0.2200 12.7800 

1.0000 0.0900 15.0000 

1.3162 0.0500 15.0000 

*Bechtel SAlC Company, LLC. “Development of Earthquake Ground Motion Input for Preclosure Seismic Design and 

Postclosure Performance Assessment of Geologic Repository at Yucca Mountain, Nevada.” MDL–MGR–GS–

000003. Rev. 01. Las Vegas, Nevada: Bechtel SAlC Company, LLC. 2004.
A–4 

Figure A–1. Variation of Normalized Shear Modulus With Shear Strain Using Data for 

Upper Mean Tuff and Upper Mean Alluvium From Bechtel SAIC Company, LLC (2004) 

Figure A–2. Variation of Material Damping Ratio With Shear Strain Using Data for Upper 

Mean Tuff and Upper Mean Alluvium From Bechtel SAIC Company, LLC (2004) 
A–5 

Figure A–3. Fit of Three Hysteretic Damping Models to Damping Data of Tuff From 

Bechtel SAIC Company, LLC (2004)

Figure A–4. Fit of Three Hysteretic Damping Models to Data of Alluvium From Bechtel 

SAIC Company, LLC (2004)