Manual

User Manual:

Open the PDF directly: View PDF PDF.
Page Count: 77

DownloadManual
Open PDF In BrowserView PDF
Strata Technical Manual
Albert R. Kottke, Xiaoyue Wang, and Ellen M. Rathje
January 22, 2018

Abstract
The computer program Strata performs equivalent-linear site response
analysis in the frequency domain using input time series or random vibration theory (RVT) methods, and allows for randomization of the site
properties. The following document explains the technical details of the
program, as well as provides a user’s guide to the program.
Strata is distributed under the GNU General Public License v31 . The
Strata source code and executables can be downloaded from: https://
github.com/arkottke/strata.

1

https://www.gnu.org/licenses/gpl-3.0-standalone.html

i

Acknowledgements
This project was sponsored by the Pacific Earthquake Engineering Research Center’s Program of Applied Earthquake Engineering Research of
Lifelines Systems supported by the California Energy Commission, California Department of Transportation, and the Pacific Gas and Electric
Company.
This work made use of the Earthquake Engineering Research Centers
Shared Facilities supported by the National Science Foundation, under
award number EEC-9701568 through the Pacific Earthquake Engineering
Research Center (PEER). Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do
not necessarily reflect those of the National Science Foundation.
Additional support was provided by the Nuclear Regulatory Commission.

ii

Contents
Abstract

i

Acknowledgements

ii

Contents

iii

List of Figures

v

List of Tables

vii

1

Introduction

2

Site Response Analysis
2.1 Equivalent Linear Site Response Analysis
2.1.1 Linear Elastic Wave Propagation .
2.1.2 Equivalent-Linear Analysis . . . .
2.1.3 Dynamic Soil Properties . . . . . .
2.2 Site Response Methods . . . . . . . . . . .
2.2.1 Time Series Method . . . . . . . .
2.2.2 Random Vibration Theory Method

3

1

Variation of Site Properties
3.1 Introduction . . . . . . . . . . . . . . . . .
3.2 Random Variables . . . . . . . . . . . . . .
3.3 Statistical Models for Soil Properties . . .
3.3.1 Layering and Velocity Model . . .
3.3.2 Depth to Bedrock Model . . . . . .
3.3.3 Non-Linear Soil Properties Model

iii

.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.

2
2
2
8
10
16
16
17

.
.
.
.
.
.

31
31
32
35
35
47
47

CONTENTS

iv

4

49
49
49
51
51
55
58
58
60
63

Using Strata
4.1 Strata Particulars . . . . . . . . . . . . . .
4.1.1 Auto-Discretization of Layers . . .
4.1.2 Interaction with Tables . . . . . . .
4.1.3 Non-Linear Curves . . . . . . . . .
4.1.4 Recorded Motion Dialog . . . . . .
4.1.5 Output Widget . . . . . . . . . . .
4.2 Examples . . . . . . . . . . . . . . . . . . .
4.2.1 Example 1: Basic time domain . .
4.2.2 Example 2: RVT and Site Variation

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

Bibliography

66

Index

69

List of Figures
2.1
2.2
2.3
2.4

2.5
2.6
2.7
2.8
2.9
2.10
2.11
2.12
2.13
2.14
3.1

The notation used in the wave equation . . . . . . . . . . . .
Nomenclature for the theoretical wave propagation. . . . . .
The input to surface transfer functions site in Table 2.1 considering different types of input. . . . . . . . . . . . . . . . .
Outcrop describes upward and downward waves being equal,
within is used when the upward and downward waves are
not equal. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
An example of the strain-time history and effective strain
(γeff ). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
An example shear-wave velocity profile . . . . . . . . . . . .
Examples of shear modulus reduction and material damping curves for soil. . . . . . . . . . . . . . . . . . . . . . . . . .
The nonlinear soil properties predicted by the Darendeli
(2001) model. . . . . . . . . . . . . . . . . . . . . . . . . . . .
The mean and mean ±σ nonlinear soil properties predicted
by Darendeli (2001). . . . . . . . . . . . . . . . . . . . . . . . .
Time domain method sequence . . . . . . . . . . . . . . . . .
The comparison between the target response spectrum and
the response spectrum computed with RVT. . . . . . . . . . .
The FAS computing through the inversion process. . . . . . .
The relative error between the computed response spectra
and the target response spectrum. . . . . . . . . . . . . . . .
RVT method sequence . . . . . . . . . . . . . . . . . . . . . .
Using a random variable from a uniform distribution between 0 and 1 to generate a variable from a distribution
with a zero mean and unit standard deviation. . . . . . . . .

v

3
4
6

7
9
11
12
14
15
18
25
26
27
29

33

LIST OF FIGURES
3.2
3.3
3.4

3.5

3.6
3.7

3.8
4.1
4.2
4.3
4.4
4.5
4.6
4.7
4.8

Two variables with a correlation coefficient of: (a) 0.0, (b),
0.99, and (c) -0.7. . . . . . . . . . . . . . . . . . . . . . . . . . .
A layering profile consisting of 8 layers modeled by a homogeneous Poisson process with a rate of 1. . . . . . . . . . .
Toro (1995) layering model. (a) The occurrence rate (λ) as a
function of depth (d), and (b) the expected layer thickness
(h) as a function of depth. . . . . . . . . . . . . . . . . . . . .
Transformation between a homogeneous Poisson process
with rate 1 to the Toro (1995) non-homogeneous Poisson
process. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
A layering simulated with the non-homogeneous Poisson
process defined by Toro (1995). . . . . . . . . . . . . . . . . .
Ten generated shear-wave velocity (vs ) profiles for a USGS
C site class. (a) Using generic layering and median vs , (b)
using user defined layering and median vs . . . . . . . . . . .
Generated nonlinear properties assuming perfect negative
correlation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Strata location selection . . . . . . . . . . . . . . . . . . . . .
By clicking on the button circled in red all rows in the table
are selected. . . . . . . . . . . . . . . . . . . . . . . . . . . . .
The nonlinear curve manager. . . . . . . . . . . . . . . . . . .
The initial view of the Recorded Motion Dialog. . . . . . .
An example a completed Recorded Motion Dialog. . . . . .
Using the Output view to examine the results of a calculation.
The shear-wave velocity profile of the Sylmar County Hospital Parking Lot site (Chang, 1996). . . . . . . . . . . . . . .
Settings to enable RVT site response and variation of the
shear-wave velocity. . . . . . . . . . . . . . . . . . . . . . . .

vi

34
37

38

40
40

46
48
50
52
54
56
57
59
60
64

List of Tables
2.1
2.2
2.3

The site properties of an example site. . . . . . . . . . . . . . 5
The values of the RVT calculation for the input motion. . . . 28
The values of the RVT calculation for the input and surface
motions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30

3.1

The categories of the geotechnical subsurface conditions
(third letter) in the Geomatrix site classification (Toro, 1995). 43
The USGS site categories, where vs,30 is the time weighted
average shear-velocity of the top 30 m (Toro, 1995). . . . . . . 43
Coefficients for the Toro (1995) model, calculated by maximum likelihood. . . . . . . . . . . . . . . . . . . . . . . . . . . 44
Median shear-wave velocity (m/s) based on the generic site
classification. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45

3.2
3.3
3.4
4.1

Soil profile at the Sylmar County Hospital Parking Lot site (Chang,
1996). The mean effective stress (σ0m ) is computed assuming
a k0 of 1/2 and a water table depth of 46 meters. . . . . . . . . 61

vii

Chapter 1
Introduction
The computer program Strata performs equivalent linear site response
analysis in the frequency domain using time domain input motions or
random vibration theory (RVT) methods, and allows for randomization of
the site properties. Strata was developed with financial support provided
by the Lifelines Program of the Pacific Earthquake Engineering Research
(PEER) Center under grant SA5405- 15811 and funding from the Nuclear
Regulatory Commission. Strata is distributed under the GNU General
Public License v31 .
The following document provides explanation of the technical details
of the program. Chapter 2 provides an introduction to equivalent-linear
elastic wave propagation using both time series and random vibration theory methods. Using the time series method, a single motion is propagated
through the site to compute the strain compatible ground motion at the
surface of the site or at any depth in the soil column. Using random vibration theory, the expected maximum response is computed from a mean
Fourier amplitude spectrum (amplitude only) and duration. Chapter 3
introduces random variables and the models that Strata uses to govern
the variability of the site properties (nonlinear properties, layering thickness, shear-wave velocity, and depth to bedrock). Chapter 4 discusses in
the interaction with Strata along with several tutorials that introduce the
features found in Strata.

1

https://www.gnu.org/licenses/gpl-3.0-standalone.html

1

Chapter 2
Site Response Analysis
Strata computes the dynamic site response of a one-dimensional soil column using linear wave propagation with strain dependent dynamic soil
properties. This is commonly referred as equivalent linear analysis method,
which was first used in the computer program SHAKE (Idriss and Sun, 1992;
Schnabel et al., 1972). Similar to SHAKE, Strata only computes the response
for vertically propagating, horizontally polarized shear waves propagated
through a site with horizontal layers.
The following chapter introduces strain dependent soil properties,
linear-elastic wave propagation through a layered medium, and the equivalent linear approach to site response analysis.

2.1
2.1.1

Equivalent Linear Site Response Analysis
Linear Elastic Wave Propagation

For linear elastic, one-dimensional wave propagation, the soil is assumed
to behave as a Kelvin-Voigt solid, in which the dynamic response is described using a purely elastic spring and a purely viscous dashpot (Kramer,
1996). The solution to the one-dimensional wave equation for a single wave
frequency (ω) provides displacement (u) as a function of depth (z) and time
(t) (?):
u(z, t) = A exp [i (ωt + k∗ z)] + B exp [i (ωt − k∗ z)]
(2.1)
In equation (2.1), A and B represent the respective amplitudes of the upward (−z) and downward (+z) waves, respectively (Figure 2.1). The com2

CHAPTER 2. SITE RESPONSE ANALYSIS

3

Figure 2.1: The notation used in the wave equation
plex wave number (k∗ ) in equation (2.2) is related to the shear modulus (G),
damping ratio (D), and mass density (ρ) of the soil using:
ω
(2.2)
k∗ = ∗
vs
s
G∗
v∗s =
(2.3)
ρ


√
G∗ = G 1 − 2D2 + ı2D 1 − D2 ' G(1 + ı2D)
(2.4)
G∗ and v∗s are called the complex shear modulus and complex shear-wave
velocity, respectively. If the damping ratio (D) is small (< 10 − 20%),
then the approximation of the complex shear modulus in equation (2.4)
is appropriate. Strata uses the complete definition of the complex shearmodulus, not the approximation, in the calculations.
Equation 2.1 applies only to a single layer with uniform soil properties and the wave amplitudes (A and B) can be computed from the layer
boundary conditions. For a layered system, shown in Figure 2.2, the
wave amplitudes are calculated using recursive formulas developed by
maintaining compatibility of displacement and shear stress at the layer
boundaries. Using these assumptions, the following recursive formulas
are developed (Kramer, 1996):




∗
∗
hm + 21 Bm 1 − α∗m exp −ıkm
hm
Am+1 = 21 Am 1 + α∗m exp ıkm




(2.5)
∗
∗
Bm+1 = 12 Am 1 − α∗m exp ıkm
hm + 21 Bm 1 + α∗m exp −ıkm
hm

CHAPTER 2. SITE RESPONSE ANALYSIS
1
2

A1 6?
B1
A2 6?
B2

4

ρ1 h1 G1 D1
ρ2 h2 G2 D2

Am 6?
m
Bm
ρm hm Gm Dm
m + 1 Am+1 6?
Bm+1 ρm+1 hm+1 Gm+1 Dm+1

n

An 6?
Bn

ρn hn Gn Dn

Figure 2.2: Nomenclature for the theoretical wave propagation.
where m is the layer number, hm is the layer height and α∗m is the complex
impedance ratio. The complex impedance ratio is defined as:
α∗m

∗
ρm v∗s,m
km
G∗m
=
= ∗
km+1 G∗m+1 ρm+1 v∗s,m+1

(2.6)

and quantifies the relative amplitudes of the upward and downward
waves. At the surface of the soil column (m = 1), the shear stress must
equal zero, therefore the amplitudes of the upward and downward waves
must be equal (A1 = B1 ).
The wave amplitudes (A and B) within the soil profile are calculated at
each frequency (assuming known stiffness and damping within each layer)
and used to computed the response at the surface of a site. This calculation is performed by setting A1 = B1 = 1.0 at the surface and recursively
calculating the wave amplitudes (Am+1 ,Bm+1 ) in successive layers until the
input (base) layer is reached. The transfer function between the motion in
the layer of interest (m) and in the rock layer (n) at the base of the deposit
is defined as:
um (ω) Am + Bm
=
(2.7)
TFm,n (ω) =
un (ω)
An + Bn
where ω is the frequency of the harmonic wave. The transfer function
is the ratio of the amplitude of motion–either displacement, velocity, or
acceleration–between two layers of interest and varies with frequency.

CHAPTER 2. SITE RESPONSE ANALYSIS
Property
Mass Density (ρ)
Height (h)
Shear-wave Velocity (vs )
Damping ratio (D)

Rock
2.24 g/cm3
inf
1500 m/s
1%

5
Soil
1.93 g/cm3
50 m
350 m/s
7%

Table 2.1: The site properties of an example site.
The transfer function for the site with the properties presented in Table 2.1
is shown in Figure 2.3. The locations of the peaks in the transfer function
are controlled by the modes of vibration of the soil deposit. The peak at the
lowest frequency represents the fundamental (i.e. first) mode of vibration
and results in the largest amplification. The peaks at higher frequencies
are the higher vibrational modes of the site. The first natural frequency of
a site is inversely related to the site period, where the site period is defined
as (Kramer, 1996):
4 · hsoil
(2.8)
Ts =
vs
In equation (2.8), hsoil is the total height of the soil and vs average velocity
of the site.
For the example site (Table 2.1), the site period is calculated to be 0.57
s which corresponds to a natural frequency of 1.75 Hz. In the transfer
function (Figure 2.3), the peak with the highest amplification occurs at this
frequency. The amplitudes of the peaks are controlled by the damping
ratio of the soil. As the damping of the system increases, the amplitudes
of the peaks decrease, which results in less amplification at the surface.
The response at the layer of interested is computed by multiplying
the Fourier amplitude spectrum of the input rock motion by the transfer
function:
Ym (ω) = TFm,n (ω) · Yn (ω)
(2.9)
where Yn is the input Fourier amplitude spectrum at layer n and Ym is the
Fourier amplitude spectrum at the top of the layer of interest. The Fourier
amplitude spectrum of the input motion can be defined using a variety of
methods and is discussed further in Sections 2.2.1 and 2.2.2.
One issue that must be considered is that the input Fourier spectrum
typically represents a motion at the ground surface, where the upgoing
and downgoing wave amplitudes are equal (A1 = B1 ), not at the base of

CHAPTER 2. SITE RESPONSE ANALYSIS

6

10
Surface / Within
Surface / Outcrop

|TF|

7.5

5

2.5

0
0

5

10

15

20

25

Frequency (Hz)
Figure 2.3: The input to surface transfer functions site in Table 2.1 considering different types of input.
a soil deposit, where the wave amplitudes are not equal (Figure 2.4). The
change in boundary conditions ( An = Bn for the surface, An , Bn at the
base of a soil deposit) must be taken into account. The motions at any
free surface are referred to as outcrop motions and their amplitudes are
described by twice the amplitude of the upward wave (2A). Equation (2.7)
can be modified to transfer an outcrop motion to a surface motion. This
result is obtained by multiplying equation (2.7) by a transfer function that
takes the outcrop motions and makes it a within motion at the base of the
soil column.
to an outcrop motion a second transfer function is required to translate
from an outcrop motion to a within motion. The combined transfer function

CHAPTER 2. SITE RESPONSE ANALYSIS

7

Figure 2.4: Outcrop describes upward and downward waves being equal,
within is used when the upward and downward waves are not equal.
is defined as:
TFm,n (ω) =

An + Bn
Am + Bm
Am + Bm
·
=
2 · An
An + Bn
2 · An
| {z } | {z }
| {z }
outcrop→within within→layern

(2.10)

outcrop→layern

Motions recorded at depth (e.g. recorded in a borehole) are referred
to as within motions and for these motions the transfer function given in
equation (2.7) can be used.
Figure 2.3 shows the transfer function for the site profile presented in
Table 2.1 using equation (2.10 where the input motion is specified as outcrop.
In comparison with the within-to-outcrop transfer function, the outcrop-tooutcrop transfer function displays less amplification for all modes.

CHAPTER 2. SITE RESPONSE ANALYSIS

2.1.2

8

Equivalent-Linear Analysis

The previous section assumed that the soil was linear-elastic. However,
soil is nonlinear, such that the dynamic properties of soil (shear modulus,
G, and damping ratio, D) vary with shear strain, and thus the intensity of
shaking. In equivalent-linear site response analysis, the nonlinear response
of the soil is approximated by modifying the linear elastic properties of the
soil based on the induced strain level. Because the induced strains depend
on the soil properties, the strain compatible shear modulus and damping
ratio values are iteratively calculated based on the computed strain.
A transfer function is used to compute the shear strain in the layer based
on the outcropping input motion. In the calculation of the strain transfer
function, the shear strain is computed at the middle of the layer (z = hm /2)
and used to select the strain compatible soil properties. Unlike the previous
transfer functions that merely amplified the Fourier amplitude spectrum,
the strain transfer function amplifies the motion and converts acceleration
into strain. The strain transfer function based on an outcropping input
motion is defined by:



∗
∗
hm /2 − Bm exp −ıkm
hm /2
γ(ω, z = hm /2) ıkm Am exp ıkm
strain
=
TFmn (ω) =
ün,outcrop (ω)
−ω2 (2 · An )
(2.11)
The strain Fourier amplitude spectrum within a layer is calculated by
applying the strain transfer function to the Fourier amplitude spectrum
of the input motion. The maximum strain within the layer is derived
from this Fourier amplitude spectrum – either through conversion to the
time domain or through RVT methods, further discussed in Section 2.2).
However, it is not appropriate to use the maximum strain within the layer
to compute the strain-compatible soil properties, because the maximum
strain only occurs for an instant in time. Instead, an effective strain (γeff )
is calculated from the maximum strain. Typically, the effective strain is
65% of the maximum strain. An example of a strain time-series and the
effective strain is shown in Figure 2.5.
Equivalent-linear site response analysis requires that the strain dependent nonlinear properties (i.e. G and D) be defined. The initial (small
strain) shear modulus (Gmax ) is calculated by:
Gmax = ρv2s

(2.12)

CHAPTER 2. SITE RESPONSE ANALYSIS

9

0.02

Shear Strain (%)

0.01
0
−0.01
−0.02
−0.03

Time Series
Effective Strain
0

10

20

30

40

Time (s)
Figure 2.5: An example of the strain-time history and effective strain (γeff ).
where ρ is the mass density of the site, and vs is the measured shear-wave
velocity. Characterizing the nonlinear behavior of G and D is achieved
through modulus reduction and damping curves that describe the variation of G/Gmax and D with shear strain (discussed in the next section). Using
the initial dynamic properties of the soil, equivalent-linear site response
analysis involves the following steps:
1. The wave amplitudes (A and B) are computed for each of the layers
2. The strain transfer function is calculated for each of the layers.
3. The maximum strain within each layer is computed by applying the
strain transfer function to the input Fourier amplitude spectrum and
finding the maximum response (see Section 2.2).
4. The effective strain (γeff ) is calculated from the maximum strain
within each layer.

CHAPTER 2. SITE RESPONSE ANALYSIS

10

5. The strain compatible shear modulus and damping ratio are recalculated based on the new estimate of the effective strain within each
layer.
6. The new nonlinear properties (G and D) are compared to the previous
iteration and an error is calculated. If the error for all layers is below
a defined threshold the calculation stops.
After the iterative portion of the program finishes, the dynamic response
of the soil deposit is computed.

2.1.3

Dynamic Soil Properties

In a dynamic system, the properties that govern the response are the mass,
stiffness, and damping. In soil under seismic shear loading, the mass of
the system is characterized by the mass density (ρ) and the layer height (h),
the stiffness is characterized by the shear modulus (G), and the damping
is characterized by the viscous damping ratio (D). The dynamic behavior
of soil is challenging to model because it is nonlinear, such that, both the
stiffness and damping of the system change with shear strain. Section 2.1.2
introduced equivalent-linear site response analysis in which the nonlinear
response of the soil was simplified into a linear system that used straincompatible dynamic properties (G and D). The analysis requires that
the strain dependence of the nonlinear properties within a layer be fully
characterized.
Defining the mass density of the system is a straight forward process because the density of soil falls within a limited range for soil and a
good estimate of the mass density can be made based on soil type alone.
Characterization of the stiffness and damping properties of soil is more
complicated, the most rigorous approach requiring testing in both the field
and laboratory.
The shear modulus and material damping of the soil are characterized
using the small strain shear modulus (Gmax ), modulus reduction curves
that relate G/Gmax to shear strain, and damping ratio curves that relate D
to shear strain. The small strain shear modulus is best characterized by in
situ measurement of the shear-wave velocity as a function of depth. An
example shear-wave velocity profile is shown in Figure 2.6. The profile
tends to be separated into discrete layers with a generally increasing shearwave velocity with increasing depth. Examples of modulus reduction and

CHAPTER 2. SITE RESPONSE ANALYSIS

11

0
−10
−20

Depth (m)

−30
−40
−50
−60
−70
−80
−90
−100
0

100

200

300
400
500
600
Shear−Wave Velocity (m/s)

700

800

Figure 2.6: An example shear-wave velocity profile
damping curves for soil are shown in Figure 2.7. These curves show a
decrease in the soil stiffness and an increase in the damping ratio with an
increase in shear strain.
The modulus reduction and damping curves may be obtained from laboratory measurements on soil samples or derived from empirical models
based on soil type and other variables. One of the most comprehensive
empirical models was developed by Darendeli (2001) and is included with
Strata. The model expands on the hyperbolic model presented by Hardin
and Drnevich (1972) and accounts for the effects of confining pressure (σ00 ),
plasticity index (PI), over-consolidation ratio (OCR), frequency ( f ), and
number of cycles of loading (N) on the modulus reduction and damping
curves.
In the Darendeli (2001) model, the shear modulus reduction curve is a
hyperbola defined by:
1
G
=
(2.13)
 
Gmax 1 + γ a
γr

where a is 0.9190, γ is the shear strain and γr is the reference shear strain.

CHAPTER 2. SITE RESPONSE ANALYSIS

12

25
1
20

Damping (%)

G/Gmax

0.8

0.6

0.4

10

5

0.2

0
−4
10

15

−3

10

−2

10
Strain (%)

−1

10

0

10

0
−4
10

−3

10

−2

10
Strain (%)

−1

10

0

10

Figure 2.7: Examples of shear modulus reduction and material damping
curves for soil.
The reference shear strain (not in percent) is computed from:
γr =

σ00
pa

!0.3483



0.0352 + 0.0010 · PI · OCR0.3246

(2.14)

where σ00 is the mean effective stress and pa is the atmospheric pressure
in atm. In the model, the damping ratio is calculated from the minimum
damping ratio at small strains (Dmin ) and from the damping ratio associated with hysteretic Masing behavior (DMasing ). The minimum damping is
calculated from:



Dmin (%) = (σ00 )−0.2889 0.8005 + 0.0129 · PI · OCR−0.1069 1 + 0.2919 ln f
(2.15)
where f is the excitation frequency (Hz). The computation of the Masing
damping requires the calculation of the area within the stress-strain curve
predicted by the shear modulus reduction curve. The integration can be
approximated by:
DMasing (%) = c1 DMasing,a=1 + c2 D2Masing,a=1 + c3 D2Masing,a=1

(2.16)

CHAPTER 2. SITE RESPONSE ANALYSIS
where:

13


 
 γ+γ  
r






γ
−
γ
ln
r



100  
γr


4
−
2
Dmasing,a=1 (%) =




2




γ


π 

 
γ+γ

(2.17)

c1 = −1.1143a + 1.8618a + 0.2533
c2 = 0.0805a2 − 0.0710a − 0.0095
c3 = −0.0005a2 + 0.0002a + 0.0003

(2.18)

r

2

The minimum damping ratio in equation (2.15) and the Masing damping
in equation (2.16) are combined to compute the total damping ratio (D)
using:


G 0.1
D=b
· DMasing + Dmin
(2.19)
Gmax
where b is defined as:
b = 0.6329 − 0.0057 ln (N)

(2.20)

where N is the number of cycles of loading. In most site response applications, the number of cycles (N) and the excitation frequency ( f ) in the
model are defined as 10 and 1, respectively. Figure 2.8 shows the predicted
nonlinear curves for a sand (PI = 0, OCR = 1) at an effective confining
pressure of 1 atm.
A Bayesian approach was used in the Darendeli (2001) model to calculate the model coefficients. One of the unique aspects of this model is that
the scatter of the data about the mean estimate is quantified. In the Darendeli (2001) model, the variability about the mean value is assumed to be
normally distributed. The normal distribution is described using a mean
and standard deviation. The mean values are calculated from equations
(2.13) and (2.19). The standard deviation is a function of the amplitude of
the nonlinear property (i.e. G/Gmax and D). The standard deviation of the
normalized shear modulus (σNG ) is computed by:
s
(G/Gmax − 0.5)2
0.25
−
(2.21)
σNG = exp(−4.23) +
exp(3.62)
exp(3.62)
This model results in small σNG when G/Gmax is close to 1 or 0 and relatively
large σNG when G/Gmax is equal to 0.5. The standard deviation of the
damping ratio (σD ) is computed by:
p
σD = exp(−5.0) + exp(−0.25) D(%)
(2.22)

CHAPTER 2. SITE RESPONSE ANALYSIS

14

25
1

σ0‘ = 1.00 atm
20

Damping (%)

G/Gmax

0.8

0.6

0.4

15

10

5

0.2

0
−4
10

OCR = 1.00
PI = 0

−3

10

−2

10
Strain (%)

−1

10

0

10

0
−4
10

−3

10

−2

10
Strain (%)

−1

10

0

10

Figure 2.8: The nonlinear soil properties predicted by the Darendeli (2001)
model.
In the damping ratio model, the σD increases with increasing damping
ratio. Using these definitions for the standard deviation, the ±σ modulus
reduction and damping curve for sand at a confining pressure of 1 atm are
shown in Figure 2.9.

CHAPTER 2. SITE RESPONSE ANALYSIS

15

25
1

σ0‘ = 1.00 atm
20

Damping (%)

G/Gmax

0.8

0.6

0.4

15

10

5

0.2

0
−4
10

OCR = 1.00
PI = 0

−3

10

−2

10
Strain (%)

−1

10

0

10

0
−4
10

−3

10

−2

10
Strain (%)

−1

10

0

10

Figure 2.9: The mean and mean ±σ nonlinear soil properties predicted by
Darendeli (2001).

CHAPTER 2. SITE RESPONSE ANALYSIS

2.2

16

Site Response Methods

The previous section introduced transfer functions which transform the
input Fourier amplitude spectrum (FAS) into a FAS of strain, acceleration,
and even a single-degree-of-freedom oscillator at any depth within the site
profile. In both the time domain and random vibration theory methods,
the same transfer functions are applied to the input FAS. The difference in
the methods is in how this FAS in the frequency domain is converted into
the time domain information.

2.2.1

Time Series Method

In the time series method, an input acceleration-time history is provided
and the input FAS is computed from that time series using the Fast Fourier
transform (FFT) to compute the discrete Fourier transformation on the
provided time series. The computed FAS is complex valued, and can be
converted into amplitude and phase information. Strata uses the free and
open-source FFTW library (http://www.fftw.org). The inverse discrete
Fourier transform is used to compute a time series for a given FAS. The
details of the FFT process are not discussed here, but can be found on the
FFTW webpage.
In Strata, the time series is padded with zeros to obtain a number of
points that is a power of two. If a time series contains a power of two
values, then it is padded with zeros until the next power of two.
The frequencies associated with the FAS are computed from the time
step between points and the number of points (N) in the record. The
highest possible sampling frequency is known as the Nyquist frequency
and is defined as:
1
(2.23)
fNyquist =
2∆t
where ∆t is the time increment between data points. The increment between the frequencies is calculated by:
∆f =

fNyquist
1
=
N/2 − 1 2∆t (N/2 − 1)

(2.24)

After the FAS of the motion has been computed it is possible to perform
site response analysis with the motion. The following is a summary of the

CHAPTER 2. SITE RESPONSE ANALYSIS

17

steps to compute the surface acceleration time-series for the site described
in Table 2.1:
1. Read the acceleration-time series file (Figure 2.10a).
2. Compute the input FAS with the Fast Fourier transformation (FFT)
(Figure 2.10b, only amplitude is shown).
3. Compute the transfer function for the site properties (Figure 2.10c,
only amplitude is shown).
4. Compute the surface FAS by applying the transfer function to the
input FAS (Figure 2.10d, only amplitude is shown).
5. Compute the surface acceleration-time series through the inverse FFT
of the surface FAS (Figure 2.10e).

2.2.2

Random Vibration Theory Method

The random vibration theory (RVT) approach to site response analysis
was first proposed in the engineering seismology literature (e.g. Schneider et al. (1991)) and has been applied to site response analysis (Rathje et al.,
2005; Silva et al., 1997). RVT does not utilize time domain input motions,
but rather initiates all computtions with the input FAS (amplitude only,
no phase information). Because RVT does not have the accompanying
phase angles to the Fourier amplitudes, a time history of motion cannot
be computed. Instead, extreme value statistics are used to compute peak
time domain parameters of motion (e.g. peak ground acceleration, spectral acceleration) from the Fourier amplitude information. Due to RVT’s
stochastic nature one analysis can provide a median estimate of the site
response with a single analysis and with the need for time domain input
motions.
The Basics of RVT
Random vibration theory can be separated into two parts: (1) conversion
between time and frequency domain using Parseval’s theorem, and (2)
estimation of the peak factor using extreme value statistics.

CHAPTER 2. SITE RESPONSE ANALYSIS

18

1
Accel (g)

(a)
0
−1

5

7.5

10

12.5

15

17.5

20

Fourier Amp. (g-s)

Time (s)

75
(b)

50
25
0
0

2.5

5

7.5

10

12.5

15

Frequency (Hz)

4
|T F |

(c)
2
0
0

2.5

5

7.5

10

12.5

15

Fourier Amp. (g-s)

Frequency (Hz)
75
(d)

50
25
0
0

2.5

5

7.5

10

12.5

15

Frequency (Hz)
Accel (g)

1
(e)
0
−1

5

7.5

10

12.5

15

17.5

20

Time (s)

Figure 2.10: Time domain method sequence: (a) input acceleration-time
series, (b) input Fourier amplitude spectrum, (c) transfer function from
input to surface, (d) surface Fourier amplitude spectrum, and (e) surface
acceleration-time series.

CHAPTER 2. SITE RESPONSE ANALYSIS

19

Consider a time varying signal x(t) with its associated Fourier amplitude spectrum, X( f ). The root-mean-squared value of the signal (xrms )
is a measure of its average value over a given time period, Trms , and is
computed from the integral of the times series over that time period:
s
Z Trms
1
[x(t)]2 dt
xrms =
(2.25)
Trms 0
Parseval’s theorem related the integral of the time series to the integral of
its Fourier Transform, such that Equation 2.25 can be written in term of the
FAS of the signal:
s
r
Z ∞
2
m0
2
|X( f )| d f =
(2.26)
xrms =
Trms 0
Trms
where m0 is defined as the zero-th moment of the FAS. The n-th moment
of the FAS is defined as:
Z ∞
n
mn = 2
2π f |X( f )|2 d f
(2.27)
0

The peak factor (PF) represents the ratio of the maximum value of the
signal (xmax ) to its rms value (xrms ), such that if xrms and PF are known, then
xmax can be computed using:
xmax = PF · xrms

(2.28)

Cartwright and Longuet-Higgins (1956) studied the statistics of ocean
wave amplitudes, and considered the probability distribution of the maxima of a signal to develop expressions for the PF in terms of the characteristics of the signal. Cartwright and Longuet-Higgins (1956) derived an
integral expression for the expected values of the peak factor in terms of
the number of extrema (Ne ) and the bandwidth (ξ) of the time series (Boore,
2003):
Z ∞
h
 iNe
√
E[PF] = 2
1 − 1 − ξ exp −z2
dz
(2.29)
0

where the bandwidth is defined as:
s
ξ=

m22
m0 m4

(2.30)

CHAPTER 2. SITE RESPONSE ANALYSIS
and the number of extrema is defined as:
r
1 m4
T gm
Ne =
π m2

20

(2.31)

Boore and Joyner (1984) illustrated the need to modify the duration
used in the rms calculation when considering oscillator responses, and
they introduced the concept of an rms duration (Trms ). The rms duration
requires modification for spectral acceleration to account for the enhanced
duration due to the oscillator response. Generally, adding the oscillator
duration to the ground motion duration will suffice, except in cases where
the ground motion duration is short (Boore and Joyner, 1984). Boore and
Joyner (1984) recommend the following expressions to define Trms :
!
γn
(2.32)
Trms = T gm + To n
γ +α
T gm
Tn
Tn
To =
2πβ
γ=

(2.33)
(2.34)

where To is the oscillator duration, Tn is the oscillator natural period, and
β is the damping ratio of the oscillator. Based on numerical simulations,
Boore and Joyner (1984) proposed n = 3 and α = 31 for the coefficients in
Equation 2.32.
Defining the Input Motion
The input motion in an RVT analysis is defined by a ground motion duration (T gm ) and a Fourier amplitude spectrum (FAS). The FAS can be directly
computed using seismological source theory (e.g (e.g. Brune, 1970, 1971)),
or it can be back-calculated from an acceleration response spectrum (see
section 2.2.2 on page 21). Strata does not provide implementation for any
seismological source theory models but does allow for the Fourier amplitude spectrum to be directly provided. The frequencies provided with the
Fourier amplitude spectrum is the frequency range used by the program
so it is critical that enough points be provided.

CHAPTER 2. SITE RESPONSE ANALYSIS

21

Calculation of the duration for use in RVT analysis can be done using
seismological source theory or empirical models. Boore (2003) recommends the following description of ground motion duration (T gm ) for the
Western United States:
T gm =

1
f0
|{z}

+

0.05R
|{z}

(2.35)

Path duration, Tp

Source duration, Ts

where R is the distance in km, and the corner frequency ( fo ) in hertz is
given by:


∆σ 1/3
f0 = 4.9 · 106 βs
M0
where ∆σ is the stress drop in bar, βs is the shear-wave velocity in units
of km/s, and M0 is the seismic moment is units of dyne-cm (Brune, 1970,
1971). The seismic moment (M0 ) is related to the moment magnitude (Mw )
by:
3
M0 = 10 2 (M0 +10.7)
For the Eastern United States, Campbell (2003) proposes that the path
duration effect be distance dependent:


0
R ≤ 10km




10km < R ≤ 70km
 0.16R
Tp = 
(2.36)

−0.03R
70km < R ≤ 130km



 0.04R
R > 130km
Empirical ground motions duration models (e.g. Abrahamson and Silva,
1996) can also be used to estimate the duration of the scenario event (T gm ).
When such a model is applied, it is recommended that T gm be taken as time
between the build up from 5% to 75% of the normalized arias intensity
(D5−75 ).
Calculation of a FAS from an Acceleration Response Spectrum
The input rock FAS (Y( f )) can be derived from an acceleration response
spectrum using an inverse technique. The inversion technique follows
the basic methodology proposed by Gasparini and Vanmarcke (1976) and
further described by Rathje et al. (2005). The inversion technique makes use

CHAPTER 2. SITE RESPONSE ANALYSIS

22

of the properties of the single-degree-of-freedom (SDOF) transfer function
used to compute the response spectral values. The square of the Fourier
amplitude at the SDOF oscillator natural frequency fn (|Y( fn )|2 ) can be
written in terms of the spectral acceleration at fn (Sa, fn ), the peak factor
(PF), rms duration of the motion (Trms ), the square of the Fourier amplitudes
(|Y( f )|2 ) at frequencies less than the natural frequency, and the integral of
the SDOF transfer function (|H fn ( f )|2 ):


 Trms S2a, fn Z fn

1
2


(2.37)
−
|Y( fn )|2 ≈ R ∞
|Y(
f
)|
d
f


2
2
2
PF
|H
(
f
)|
d
f
−
f
0
f
n
n
0
Within Equation 2.37, the integral of the transfer function is constant for a
given natural frequency and damping ratio (β), allowing the equation to
be simplified to (Gasparini and Vanmarcke, 1976):


 Trms S2a, fn Z fn

1
2

 
|Y( fn )|2 ≈ 
−
|Y(
f
)|
d
f
(2.38)

2
π
2
PF
0
fn 4β − 1
The peak factors in Equation 2.38 depend on the moments of the FAS
which is currently undefined. So the peak factors for all natural frequencies
are initially assumed to be 2.5. Equation 2.38 is applied first to the spectral
acceleration of the lowest frequency (longest period) provided by the user.
At this frequency, the FAS integral term in Equation 2.38 can be assumed
to be equal to zero. The equation is then applied at successively higher
frequencies using the previously computed values of |Y( f )| to assess the
integral.
After an initial estimate of the FAS is developed using Equation 2.38
with assumed constant peak factors, it would be possible to recompute
the peak factors for each period. The variable peak factors would provide
a better estimate of the FAS from the target response spectrum. This
approach was originally implemented, but the FAS based on the variable
peak factors still resulted in a response spectrum that was more than 10%
different than the target response spectrum at short periods.
To improve the agreement between the RVT-derived FAS (and associated response spectrum) and the target response spectrum by correcting
the FAS by the ratio of the two response spectra. This correction technique
is possible because of the narrow band property of the SDOF transfer function. Using an iterative process, the FAS from iteration i is corrected by the

CHAPTER 2. SITE RESPONSE ANALYSIS

23

ratio of the spectral acceleration computed by RVT (SRVT
( f )) to the target
a
target
(Sa ( f )) at each frequency:
|Yi+1 ( f )| =

SRVT
(f)
a
target

Sa

(f)

|Yi+1 ( f )|

(2.39)

The following process is used in the generate a corrected FAS:
1. Initial FAS is computed using the Vanmarcke (1983) technique.
2. The acceleration response spectrum associated with this FAS is computed using RVT.
3. The FAS is corrected using equation (2.39).
4. Using the corrected FAS, a new acceleration response spectrum is
calculated.
This process is repeated with three conditions:
1. maximum of 30 iterations,
2. a root-mean-square-error of 0.005 is achieved, or
3. the change in the root-mean-square-error is less than 0.001.
This ratio correction works very well in producing a FAS that agrees with
the target response spectrum, but the resulting FAS may have an inappropriate shape at some frequencies, as discussed below.
To demonstrate the inversion process consider a scenario event of a
magnitude 6.5 earthquake with a strike-slip faulting mechanism at a distance of 20 km. The target response spectrum is computed using the
Abrahamson and Silva (1997) attenuation model (Figure 2.11). An initial
estimate of the FAS is computed using the Gasparini and Vanmarcke (1976)
method and then the ratio correction algorithm is applied. This methodology results in a good agreement – less than 5% relative error as shown
in Figure 2.13 – with the target response spectrum (Figure 2.11), but the
associated FAS tends to curl up at low and high frequencies (curve labeled
“Ratio Corrected” in Figure 2.12). The curling up at low frequencies can
be mitigated by extending the frequency domain.

CHAPTER 2. SITE RESPONSE ANALYSIS

24

The frequency domain is extended to half of the minimum frequency
and twice the maximum frequency specified in the target response spectrum. For example, if the target response spectrum is provided from 0.2 to
100 Hz (5 to 0.01 seconds), then the frequencies of the FAS are defined to be
2048 points equally spaced in log space from 0.1 to 200 Hz. The resulting
response spectrum shows even better agreement with a maximum error of
about 3% (curve labeled “Ratio and Extrapolated” in Figure 2.13).
Theoretically the slope of the FAS at high frequencies should be increasingly negative due to a path-independent loss of the high-frequency
motion (Boore, 2003). However, in the computed FAS the slope actually
increases at around 50 Hz (Figure 2.12). At first glance it may not appear
that this relatively small amount of high frequency energy is important.
However, in site response analysis the damping in the soil further attenuates this high frequency portion of the FAS, and the peak factor depends on
the 4th moment which is sensitive to this high frequency energy. The result
is that the strain (or acceleration) near the surface – after attenuation of the
high frequency energy has occurred – is significantly less than it computed
using traditional time domain methods. To solve for this shortcoming, the
slope of the FAS at high frequencies is forced down (Figure 2.12). However, this solution results in a slight under prediction of the peak ground
acceleration (Figure 2.11 and 2.13).

CHAPTER 2. SITE RESPONSE ANALYSIS

25

1

Spectral Accel. (g)

0.5
0.2
0.1
0.05
0.02
0.01
0.01

Ratio Corrected
Ratio & Extrapolated
Ratio, Extrap., & Slope Forced
Target
0.02

0.05

0.1

0.2

0.5

1

2

5

10

Period (s)

Figure 2.11: The comparison between the target response spectrum and
the response spectrum computed with RVT.

CHAPTER 2. SITE RESPONSE ANALYSIS

26

|FAS| (g-s)

10−1

10−2

10−3
Ratio Corrected
Ratio & Extrapolated
Ratio, Extrap., & Slope Forced
10−4
0.1

1

10

100

1000

Frequency (Hz)

Figure 2.12: The FAS computing through the inversion process.

CHAPTER 2. SITE RESPONSE ANALYSIS

27

6

Relative Error (%)

4
2
0
Ratio Corrected
Ratio & Extrapolated
Ratio, Extrap., & Slope Forced

−2
−4
0.01

0.02

0.05

0.1

0.2

0.5

1

2

5

10

Period (s)

Figure 2.13: The relative error between the computed response spectra and
the target response spectrum.

CHAPTER 2. SITE RESPONSE ANALYSIS

28

Example of the RVT Procedure
The following is an example of the random vibration theory based site
response analysis to estimate the peak acceleration at the top of the site
described in Table 2.1: The earthquake scenario is a magnitude 6.5 event
at a distance of 20 km, as described in the previous section.
1. Empirical relationships are used to specify the input rock response
spectrum (Figure 2.11) and ground motion duration (T gm = D5−75 =4.55s).
2. Using the inversion technique, the FAS corresponding to the target
response spectrum is computed (Figure 2.14a). In this example, the
maximum acceleration of the input motion is computed with RVT to
allow for a comparison in the peak response between the surface and
the input. The RVT calculation results are shown in Table 2.2.
Parameter
Moments of FAS (m0 , m2 , and m4 )

Bandwidth (ξ)
Number of extrema (Ne )
Peak factor (PF)
Root-mean-square acceleration (arms )
Expected maximum acceleration (amax )

Input Motion Value
0.0280,
93.8435, and
1.7382 ·107
0.1346
623.3158
3.1406
0.0784 g
0.2462 g

Equation
2.27

2.30
2.31
2.29
2.26
2.28

Table 2.2: The values of the RVT calculation for the input motion.
3. Compute the transfer function for the site properties (Figure 2.14b).
4. Compute the surface FAS by applying the absolute value of the transfer function to the input FAS (Figure 2.14c). Using the surface FAS
the maximum expected acceleration can be computed as presented
in Table 2.3. The calculation shows that the site response increases
the expected peak ground acceleration by approximately 34%.

CHAPTER 2. SITE RESPONSE ANALYSIS

29

Fourier Amp. (g-s)

100
(a)
10−1
10−2
10−3
10−4
0.1

0.2

0.5

1

2

5

10

20

50

100

Frequency (Hz)

10
(b)

|T F |

1

0.1

0.01
0.1

0.2

0.5

1

2

5

10

20

50

100

Frequency (Hz)
0

Fourier Amp. (g-s)

10

(c)
10−1
10−2
10−3
10−4
0.1

0.2

0.5

1

2

5

10

20

50

100

Frequency (Hz)

Figure 2.14: Random vibration theory method sequence: (a) input Fourier
amplitude spectrum, (b) transfer function from input to surface, and (c)
surface Fourier amplitude spectrum.

0.2462 g

0.3375 g

Value
Input Motion Surface Motion
0.0280,
0.0635,
93.8435, and 39.6356, and
1.7382 ·107
1.6306 ·105
0.1346
0.3895
623.3158
92.8944
3.1406
2.8568
0.0784 g
0.1181 g

2.28

2.30
2.31
2.29
2.26

Equation
2.27

Table 2.3: The values of the RVT calculation for the input and surface motions.

Expected maximum acceleration (amax )

Bandwidth (ξ)
Number of extrema (Ne )
Peak factor (PF)
Root-mean-square acceleration (arms )

Parameter
Moments of FAS (m0 , m2 , and m4 )

CHAPTER 2. SITE RESPONSE ANALYSIS
30

Chapter 3
Variation of Site Properties
3.1

Introduction

A soil profile consists of discrete layers that vary in thickness based on
the properties of the soil. The layers are typically discretized based on the
soil type, recorded from borehole samples or inferred from a shear-wave
velocity profile. In seismic site response analysis, each layer is characterized by a thickness, mass density, shear-wave velocity, and nonlinear
properties (G/Gmax , and D). One of the challenges in defining values for
these properties is the natural variability across a site and the uncertainty
in their measurement. Because the dynamic response of a site is dependent
on the soil properties, any variation in the soil properties will change both
the expected surface motion and its standard deviation.
In a simple system, the variability of the components can be analytically combined to quantify the variability of the complete system, thus
allowing for the expected value and variability of the system response to
be computed. In seismic site response analysis, the nonlinear response of
the system does not allow an exact analytic quantification of the variability
of the site response. Instead, an estimate of the expected surface response
and its standard deviation due to variations in the soil properties can be
made through Monte Carlo simulations. Monte Carlo simulations estimate
the response of a system by generating parameters of the system based on
defined statistical distributions and computing the response for each set
of input parameters. The following chapter introduces Monte Carlo simulations as applied to site response analysis and presents the models that

31

CHAPTER 3. VARIATION OF SITE PROPERTIES

32

describe the variability of the layering, shear-wave velocity, and nonlinear
properties (G/Gmax and D).

3.2

Random Variables

The goal of a Monte Carlo simulation is to estimate the statistical properties
of the response of a complex system. To achieve this goal, each of the
properties of the system is selected from defined statistical distributions
and the response of the system is computed. The response is computed
for many realizations and the calculated response from each realization is
then used to estimate statistical properties of the system’s response. While
Monte Carlo simulations can be used on a wide variety of problems, a
major disadvantage is that a large number of simulations is required to
achieve stable results.
Monte Carlo simulations require that each of the components in the
system has a complete statistical description. The description can be in the
form of a variety of statistical distributions (i.e. uniform, triangular, normal, log-normal, exponential, etc.), however the normal and log-normal
distributions typically are used because they can be easily described using
a mean (µ) and standard deviation (σ). For normally distributed variables,
a random value (x) can be generated by:
x = σx · ε + µx

(3.1)

where µx is the mean value, σx is the standard deviation, and ε is a random
variable with zero mean and unit standard deviation.
There are a variety of methods for generating random variables for a
given distribution. The most general technique involves using the cumulative density function (CDF) to convert between probability and value.
For example, the calculation of a normally distributed random variable is
achieved by first calculating a random variable from a uniform distribution
between 0 and 1. The random variable (ε) is then calculated by finding the
inverse to this random value on the standard (σ=1) normal CDF, shown
in Figure 3.1. If the distribution of the variable is not normal, a similar
technique of using a random value from a unit uniform distribution and
the CDF of the probability distribution of the variable can be used (Ang
and Tang, 1990). While this method is applicable for a set of independent

CHAPTER 3. VARIATION OF SITE PROPERTIES

33

Cumulative Probability

1
0.8
0.6
0.4
0.2
0
−3

−2

−1

0
ε

1

2

3

Figure 3.1: Using a random variable from a uniform distribution between
0 and 1 to generate a variable from a distribution with a zero mean and
unit standard deviation.
random variables, it may not be used when the variables are related (i.e.
correlated).
In the case of random variables that are not independent, a more complicated procedure is required for the generation of values. Before discussing
the technique for generating correlated random variables, the concepts of
correlation and linear functions of random variables must be introduced.
Consider random variables x1 and x2 . The covariance of the two random
variables is defined as:
cov(x1 , x2 ) = E[(x1 − µx1 )(x2 − µx2 )] = E[x1 x2 ] − E[x1 ]E[x2 ]

(3.2)

where E is the expected value and µx1 and µx2 are the mean values of x1
and x2 , respectively (Ang and Tang, 1975). The covariance quantifies the
strength of relationship between x1 and x2 . If the variables are independent
of each other (Figure 3.2a), then the covariance is zero, however a covariance of zero does not necessarily indicate the variables are independent.
Instead, it indicates that the variables do not have a linear dependence.
As the covariance becomes more positive, two variables have a greater
tendency to both differ from their respective mean values in the same
direction (Figure 3.2b). Conversely, as the covariance becomes more negative, variables have a greater tendency to differ in the opposite direction

CHAPTER 3. VARIATION OF SITE PROPERTIES

4

4

2

2

2

0
−2
−4
−4

X2

4

X2

X2

ρ = −0.7

ρ = 0.99

ρ = 0.0

34

0
−2

−2

0
X1

2

(a)

4

−4
−4

0
−2

−2

0
X1

2

4

−4
−4

(b)

−2

0
X1

2

4

(c)

Figure 3.2: Two variables with a correlation coefficient of: (a) 0.0, (b), 0.99,
and (c) -0.7.
( Figure 3.2c). The covariance matrix ([C]) of a set of random variables is
defined as:
[C]i,j = cov(Xi , X j )
(3.3)
For two variables, the covariance matrix expands to:
#
"
# "
σ2x1
ρx1 ,x2 σx1 σx2
cov(x1 , x1 ) cov(x1 , x2 )
[C] =
=
cov(x2 , x1 ) cov(x2 , x2 )
ρx1 ,x2 σx1 σx2
σ2x2

(3.4)

where σx1 and σx2 are the standard deviations of x1 and x2 , respectively, and
ρx1 ,x2 is the correlation coefficient, defined as:
ρx1 ,x2 =

cov(x1 , x2 ) E[x1 x2 ] − E[x1 ]E[x2 ]
=
σx1 σx2
σx1 σx2

(3.5)

The correlation coefficient can range -1 to 1.
Independent random variables from a normal distribution are generated using equation (3.1) independently for each random variable. By combining the multiple applications equation (3.1) into a system of equations,
the generation of two independent variables is achieved by multiplying a
vector of random variables (~
ε) by a matrix ([σ]) and adding a constant (~
µ),
defined as:
(
) "
#(
) (
)
x1
σx1 0
ε1
µ1
=
+
(3.6)
x2
0 σx2
ε2
µ2

CHAPTER 3. VARIATION OF SITE PROPERTIES

35

where ε1 and ε2 are random variables randomly selected from a standard
normal distribution (µ = 0 and σ = 1), σx1 and σx2 are the standard deviations of x1 and x2 , respectively, and µx1 and µx2 are the mean values of x1 and
x2 , respectively. Because the random variables x1 and x2 are independent
(ρx1 ,x2 = 0), the off-diagonal values in the matrix ([σ]) are zero.
Using the same framework, a linear system of equations is used to
generate a pair of correlated random variables. However, because of the
correlation between x1 and x2 the off diagonal values in the matrix will
no longer be zero. Instead, a pair of correlated random variables (~
x) are
generated by:
(
) (
)
(
)  σ
0
 ε1

x1
µ1
x1
q



+
(3.7)
= 
µ2
x2
ρx1 ,x2 σx2 σx2 1 − ρ2x1 ,x2  ε2
Here, the first random variable (x1 ) is calculated based on the value of ε1
alone, while the second random variable (x2 ) is a function of both ε1 and
ε2 . Note that ε1 and ε2 still represent random and independent variables
generated from the standard normal distribution.

3.3

Statistical Models for Soil Properties

For the properties of the soil to be randomized and incorporated into Monte
Carlo simulations, the statistical distribution and properties of the soil need
to be characterized. In this research, two separate models are used. The
first model, developed by Toro (1995), describes the statistical distribution
and correlation between layering and shear-wave velocity. The second
model by Darendeli (2001) was previously introduced in Section 2.1.3 and
is used to describe the statistical distribution of the nonlinear properties
(G/Gmax and D).

3.3.1

Layering and Velocity Model

In Strata, the randomization of the layering and the shear-wave velocity is
done through the use of the velocity profile model proposed by Toro (1995).
The Toro (1995) models provides a framework for generating layering and
then to vary the shear-wave velocity of these layers. This model improves
upon previous work by quantifying the correlation between the velocities

CHAPTER 3. VARIATION OF SITE PROPERTIES

36

in adjacent layers. In previous models, one of two assumptions were made
that simplified the problem: the velocities at all depths were perfectly
correlated and could be randomized by applying a constant random factor
to all velocities (McGuire et al., 1989; Toro et al., 1992), or the velocities
within each of the layers are independent of each other, and therefore
can be randomized by applying an independent random factor to each
layer (Costantino et al., 1991). While these two assumptions simplify the
problem, they represent two extreme conditions. The Toro (1995) model
makes neither of these assumptions, instead the model incorporates a
correlation between layers.
Layering Model
The layering is modeled as a Poisson process, which is a stochastic process
with events occuring at a given rate (λ). For a homogeneous Poisson
process this rate is constant, while for a non-homogeneous Poisson process
the rate varies. Generally, a Poisson process models the occurrence of
events over time, but for the layering problem the event is a layer interface
and its rate is defined in terms of length (i.e., number of layer interfaces
per meter).
In the Toro (1995) model, the layering thickness is modeled as nonhomogeneous Poisson process where the rate changes with depth (λ(d),
where d is depth). Before considering the non-homogeneous Poisson process, first consider the simpler homogeneous Poisson process with a constant rate. For a Poisson process with a constant occurrence rate (λ), the
distance between layer boundaries, or layer thickness (h), is an exponential
distribution with rate λ. The probability density function of an exponential
distribution is defined as:
(
λ exp(−λh) for h ≥ 0
f (h; λ) =
(3.8)
0
for h < 0
The cumulative density function for the exponential distribution is given
by:
(
1 − exp(−λh) for h ≥ 0
F(h; λ) =
(3.9)
0
for h < 0
A random layer thickness with an exponential distribution is generated by

CHAPTER 3. VARIATION OF SITE PROPERTIES

37

0

Cumulative Depth

2.5
5
7.5
10
12.5
0

2

4

6

8

10

Layer Number

Figure 3.3: A layering profile consisting of 8 layers modeled by a homogeneous Poisson process with a rate of 1.
solving Equation 3.9 with respect to thickness (t):
ln [1 − F(h)]
for 0 < F(h) ≤ 1
(3.10)
−λ
By randomly generating probabilities (F(t)) with a uniform distribution
and computing the associated thicknesses with Equation 3.10, a layering
profile was simulated for 8 layers with a rate of 1 shown in Figure 3.3.
Simulating a non-homogeneous Poisson process is achieved by first
simulating a homogeneous Poisson process and then warping the homogeneous layer boundaries (S0n ) into the non-homogeneous equivalent (Sn ).
The process is very similar to the generation of random variables with a
specific distribution where a uniform distribution is warped into a particular distribution using the CDF. For a non-homogeneous Poisson process
with rate λ(d) the cumulative rate (Λ(d)) is defined as:
Z d
λ(s) ds
(3.11)
Λ(d) =
h=

0

Λ(d) represents the expected number of layers up to a depth d. The nonhomogeneous process is simulated by generating exponential distributions

CHAPTER 3. VARIATION OF SITE PROPERTIES

38

with a unit rate (λ = 1) and mapping them to the non-homogeneous space
using the inverse of Equation 3.11 (Λ−1 (d)).
As an example, the warping technique will be used on the homogeneous Poisson process to take a unit rate exponential random variables
and convert them into a non-unit rate
Before Λ−1 (d) can be defined, Λ(d) and λ(d) must be defined.
Toro (1995) proposed the following generic depth dependent rate model:
λ(d) = a · (d + b)c

(3.12)

The coefficients a, b, and c were estimated by Toro (1995) using the method
of maximum likelihood applied to the layering measured at 557 sites,
coming mostly from California. The resulting values of a be 1.98, 10.86,
and -0.89, respectively. The occurrence rate (λ(d)) quickly decreases as
the depth increases (Figure 3.4(a)). This decrease in the occurrence rate
increases the expected thickness of deeper layers. The expected layer
thickness (h) is equal to the inverse of the occurrence rate (λ(d)) and is
shown in Figure 3.4(b). The expected thickness ranges from 4.2 m at the
surface to 59 m at a depth of 200 m.
60
(a)

Expected Thickness, h (m)

Occurance Rate, λ (1/m)

0.25
0.2
0.15
0.1
0.05
0

(b)

40

20

0
0

50

100
Depth (m)

150

200

0

50

100

150

200

Depth (m)

Figure 3.4: Toro (1995) layering model. (a) The occurrence rate (λ) as a
function of depth (d), and (b) the expected layer thickness (h) as a function
of depth.

CHAPTER 3. VARIATION OF SITE PROPERTIES

39

Using Equations 3.11 and 3.12 the cumulative rate for the Toro (1995)
modeled is defined as:
#
"
Z d
(d + b)c+1
bc+1
c
Λ(d) =
−
(3.13)
a · (s + b) ds = a ·
c+1
c+1
0
The inverse cumulative rate function is then defined as:
1
! c+1
d
c
·
d
Λ−1 (d) =
+ + bc+1
−b
a
a

(3.14)

Using this equation the homogeneous Poisson process (Figure 3.3) can be
warped into a non-homogeneous Poisson process as shown in Figure 3.5.
The resulting depth profile is shown in Figure 3.6.
Velocity Model
After the layering of the profile has been established, the shear-wave velocity profile can be generated. In the Toro (1995) model, the shear-wave
velocity at mid-depth of the layer is described by a log-normal distribution.
The standard normal variable (Z) of the ith layer is calculated by:
Zi =

ln(Vi ) − ln[Vmedian (di )]
σln Vs

(3.15)

where Vi is the shear-wave velocity in the ith layer, Vmedian (di ) is the median
shear-wave velocity at mid-depth of the layer and σln Vs is the standard
deviation of the natural logarithm of the shear-wave velocity. Equation 3.15
is then solved for the shear-wave velocity of the ith layer (Vi ):
Vi = exp {σln

V

· Zi + ln [Vmedian (di )]}

(3.16)

Equation 3.16 allows for the calculation of the velocity within a layer for
a given median velocity at the mid-depth of the layer, standard deviation,
and standard normal variable. In the model proposed by Toro (1995),
values for median velocity versus depth (Vmedian (d)) and standard deviation
(σln Vs ) are provided based on site class. However, in the implementation of
the Toro (1995) model in Strata the median shear-wave velocity is defined
by the user. The standard normal variable of the ith layer (Zi ) is correlated
with the layer above it, and this inter-layer correlation is also dependent on

CHAPTER 3. VARIATION OF SITE PROPERTIES

40

12.5

Unit Depth

10
7.5
5
2.5
0
0

50

100

150

200

250

300

Transformed Depth

Figure 3.5: Transformation between a homogeneous Poisson process with
rate 1 to the Toro (1995) non-homogeneous Poisson process.

Cumulative Depth

0

100

200

300

400
0

2

4

6

8

10

Layer Number

Figure 3.6: A layering simulated with the non-homogeneous Poisson process defined by Toro (1995).

CHAPTER 3. VARIATION OF SITE PROPERTIES

41

the site class. The standard normal variable (Zi ) of the shear-wave velocity
in the top layer (i = 1) is independent of all other layers and is defined as:
Z1 = ε1

(3.17)

where ε1 is an independent normal random variable with zero mean and a
unit standard deviation. The standard normal variables of the other layers
in the profile are calculated by a recursive formula, defined as:
q
Zi = ρ · Zi−1 + εi · 1 − ρ2
(3.18)
where Zi−1 is the standard normal variable of the previous layer, εi is a new
normal random variable with zero mean and unit standard deviation, and
ρ is the inter-layer correlation.
Correlation is a measure of the strength and direction of a relationship
between two random variables. The inter-layer correlation between the
shear-wave velocities proposed by Toro (1995) is a function of both the
depth of the layer (d) and the thickness of the layer (h):


ρ(d, h) = 1 − ρd (d) · ρh (h) + ρd (d)
(3.19)
where ρh is the thickness dependent correlation and ρd is the depth dependent correlation. The thickness dependent correlation is defined as:
!
−h
ρh (h) = ρ0 · exp
(3.20)
∆
where ρ0 is the initial correlation and ∆ is a model fitting parameter. As
the thickness of the layer increases, the thickness-dependent correlation
decreases. The depth dependent correlation (ρd ) is defined as a function of
depth (d):

i
h

d+d0 b

 ρ200 · 200+d
for d ≤ 200
0
ρd (d) = 
(3.21)

 ρ200
for d > 200
where ρ200 is the correlation coefficient at 200 m and d0 is an initial depth
parameter. As the depth of the layer increases, the depth-dependent correlation increases. The final layer in a site response model is assumed to
be infinitely thick, therefore the correlation between the last soil layer and
the infinite half-space is only dependent on ρd . (Toro, 1995) evaluated each

CHAPTER 3. VARIATION OF SITE PROPERTIES

42

of the parameters in the correlation models (ρ0 , ρ200 , ∆, d0 , b) for different
generic site classes.
A site class is used to categorize a site based on the shear-wave velocity profile and/or local geology. In the Toro (1995) model, the statistical
properties of the soil profile (the median velocity, standard deviation, and
layer correlation) are provided for two different classifications schemes,
the Geomatrix and USGS classifications. The Geomatrix site classification
classifies sites based on a general description of the geotechnical subsurface
conditions, distinguishing generally between rock, shallow soil, deep soil,
and soft soil (Table 3.1). In contrast, the USGS site classification is based
on the time-weighted average shear-wave velocity (see equation ?? of the
top 30 meters (vs,30 ) (Table 3.2), and requires site specific measurements of
shear-wave velocity.
Toro (1995) computed the statistical properties of the profiles for both
the Geomatrix and USGS classifications using a maximum-likelihood procedure. The procedure used a total of 557 profiles, with 541 profiles for the
USGS classification and only 164 profiles for the Geomatrix classification.
The correlation parameters (ρ0 , ρ200 , ∆, d0 , b) are presented in Table 3.3 and
the median shear-wave velocities in are presented in Table 3.4.
Ten generated shear-wave velocity profiles were created for a deep,
stiff alluvium site using the two previously discussed methods. In the first
method, a generic site profile is generated by using the layering model
coefficients and median shear-wave velocity for a USGS C site class, shown
in Figure 3.7(a). This approach essentially models the site as a generic
USGS class C site, which is the general site classification of the for deep,
stiff alluvium. The second method uses the layer correlation for the USGS
C site class, but the layering and the median shear-wave velocity profile are
defined from field measurements, shown in Figure 3.7(b). The site specific
layering tends to be much thicker than the generic layering as a result of
the field measurements indicating thick layers with the same shear wave
velocity. In general both of the methods show an increase in the shearwave velocity with depth. However, the site-specific shear-wave velocity
values are significantly larger than the generic shear-wave velocity values.
At the surface, the generic site has a median shear-wave velocity of 150 m/s
compared to the site specific shear-wave velocity of 200 m/s. At a depth of
90 m, the difference is even greater, with the generic site having a median
shear-wave velocity of 470 m/s compared to the site specific median shearwave velocity of 690 m/s. The difference in shear-wave velocity is a result of

CHAPTER 3. VARIATION OF SITE PROPERTIES

43

Designation Description
A
Rock
Instrument is found on rock material (vs > 600 m/s) or
a very thin veneer (less than 5 m) of soil overlying rock
material.
B

Shallow (Stiff) Soil
Instrument is founded in/on a soil profile up to 20
m thick overlying rock material, typically a narrow
canyon, near a valley edge, or on a hillside.

C

Deep Narrow Soil
Instrument is found in/on a soil profile at least 20 m
thick overlying rock material in a narrow canyon or
valley no more than several kilometers wide.

D

Deep Broad Soil
Instrument is found in/on a soil profile at least 20 m
thick overlaying rock material in a broad canyon or
valley.

E

Soft Deep Soil
Instrument is found in/on a deep soil profile that exhibits low average shear-wave velocity (vs < 150 m/s).

Table 3.1: The categories of the geotechnical subsurface conditions (third
letter) in the Geomatrix site classification (Toro, 1995).
Designation
A
B
C
D

Average Shear-Wave Velocity (vs,30 )
greater than 750 m/s
360 to 750 m/s
180 to 360 m/s
less than 180 m/s

Table 3.2: The USGS site categories, where vs,30 is the time weighted average
shear-velocity of the top 30 m (Toro, 1995).

CHAPTER 3. VARIATION OF SITE PROPERTIES

44

Site Classification
GeoMatrix
USGS
Parameter A & B C & D A & B C & D
A
B
σlnV
0.46
0.38
0.35
0.36
0.36 0.27
ρ0
0.96
0.99
0.95
0.99
0.95 0.97
1.00
1.00
0.42 1.00
ρ200
0.96
1.00
∆
13.1
8.0
4.2
3.9
3.4
3.8
0.0
0.0
0.0
0.0
d0
0.0
0.0
b
0.095 0.160 0.138 0.293 0.063 0.293
Profiles
Layers

45
243

109
692

204
280

253
1487

35
129

169
750

C
0.31
0.99
0.98
3.9
0.0
0.344

D
0.37
0.00
0.50
5.0
0.0
0.744

226
1349

27
136

Table 3.3: Coefficients for the Toro (1995) model, calculated by maximum
likelihood.
the difference between the site specific information and the generic USGS
site C median shear-wave velocity profile.

CHAPTER 3. VARIATION OF SITE PROPERTIES

Median Shear-Wave Velocity (m/s)
GeoMatrix
USGS
Depth (m) A & B C & D A & B C & D
A
B
0.00
192
144
182
147
314 159
221
164
346 200
1.00
209
159
2.00
230
178
262
178
384 241
297
188
430 275
3.00
253
193
330
193
485 308
4.00
278
204
5.00
303
211
362
196
550 337
6.00
329
217
390
200
624 361
7.20
357
228
412
209
703 382
437
218
789 404
8.64
395
240
468
228
880 433
10.37
443
253
12.44
502
270
504
248
973 467
540
273
1070 501
14.93
575
291
17.92
657
319
578
296
1160 535
615
317
1260 567
21.50
748
357
653
347
1330 605
25.80
825
402
30.96
886
444
702
374
1380 654
37.15
942
474
734
386
1420 687
759
394
1460 711
44.58
998
495
53.20
1060
516
782
403
1500 732
805
427
749
64.20
541
77.04
566
834
459
772
92.44
593
870
488
802
110.93
922
515
847
133.12
983
550
900
159.74
604
682
191.69
230.03
770

45

C
145
163
179
191
200
208
215
226
237
250
269
291
314
336
372
391
401
408
413
433
459
486
513
550
604
676
756

D
176
165
154
142
129
117
109
106
109
117
130
148
170
192
210
229
246
266
289
318
353
392
435

Table 3.4: Median shear-wave velocity (m/s) based on the generic site
classification.

CHAPTER 3. VARIATION OF SITE PROPERTIES

Site Specific Layering and Velocity

0

0

−10

−10

−20

−20

−30

−30

−40

−40
Depth (m)

Depth (m)

Generic Layering and Velocity (USGS C)

−50

−50

−60

−60

−70

−70

−80

−80

−90

−90

−100

46

Input Median
Randomization Median
Randomization

−100
0

200
400
600
800
Shear−wave Velocity (m/s)
(a)

1000

0

200
400
600
800
Shear−wave Velocity (m/s)
(b)

1000

Figure 3.7: Ten generated shear-wave velocity (vs ) profiles for a USGS C
site class. (a) Using generic layering and median vs , (b) using user defined
layering and median vs .

CHAPTER 3. VARIATION OF SITE PROPERTIES

3.3.2

47

Depth to Bedrock Model

...

3.3.3

Non-Linear Soil Properties Model

The Darendeli (2001) empirical model for nonlinear soil properties (G/Gmax
and D) was previously discussed in Section 2.1.3. The Darendeli (2001)
empirical model assumes the variation of the properties is normal in distribution. The standard deviation of the G/Gmax and the D varies with
the magnitude of the property and is calculated with equations (2.21) and
(2.22), respectively. Because the variation of the properties are modeled
with a normal distribution is continuous from −∞ to +∞, the generated
values of the shear-modulus reduction or damping ratio may fall below
zero. The most likely location for the negative values occurs when the
mean value is small, which occurs at large strains for the G/Gmax and at
low strains for D. Negative values for either G/Gmax or D are not physically possible, therefore the normal distributions need to be truncated. To
correct for this problem, minimum values for G/Gmax and D are defined as
0.05 and 0.1%, respectively. These values were chosen as they represent
appropriate minimum values. The influence of the minimum values on the
site response results will be minimal because the truncation occurs only at
the extremes of the strain range.
The G/Gmax and D curves are not independent of each other. Consider
a soil that behaves more linearly, that is to say the G/Gmax is higher than
the median G/Gmax . During a loading cycle, the area inside the hysteresis
would be smaller which is indicative of less damping within the system.
As the linearity of the system increases, the damping decreases. To capture
this effect, the soil properties are assumed to have a negative correlation
with the default value set at -0.5 (i.e. ρ = −0.5). Using a correlation
coefficient of -0.5, the nonlinear properties of sand (PI=0, OCR=0) at a
confining pressure of 1 atm were generated 10 times, shown in Figure 3.8.
Three of the realizations result in large shear modulus reduction curve
relative to the mean. Because of the negative correlation, the relative high
shear modulus reduction corresponds to a relative low damping ratio.

CHAPTER 3. VARIATION OF SITE PROPERTIES

48

25
1

σ0‘ = 1.00 atm
20
Damping (%)

G/Gmax

0.8
0.6
0.4
Input Median
Randomization Median
Randomization

0.2
0
−4
10

−3

10

−2

10
Strain (%)

OCR = 1.00
PI = 0

15
10
5

−1

10

0

10

0
−4
10

−3

10

−2

10
Strain (%)

−1

10

0

10

Figure 3.8: Generated nonlinear properties assuming perfect negative correlation.

Chapter 4
Using Strata
Strata is introduced through two examples that demonstrates the organization and most of the features found in Strata. Before these two examples
some particular differences between Strata and other site response programs are introduced.
To someone that is used to working with text input files, operating
Strata will seem a little foreign. With the exception of acceleration timeseries, all of the input to Strata is entered via the keyboard, or through
copying and pasting from spreadsheets. The input file is not saved in the
typical text format, instead a binary format is used that is only readable by
Strata. Furthermore, when the calculation is complete no output text files
are produced. Instead, the output can be directly view with Strata and
saved, once again to a binary format. There is the option for the data to be
exported to text files that can then be opened with a variety of applications.

4.1
4.1.1

Strata Particulars
Auto-Discretization of Layers

One of the biggest differences between Strata and other site response analysis programs is the fact that the sublayers used in the calculation portion of
the analysis are not defined by the user. Instead, the user defines a velocity
layer that is then subdivided into sublayers by Strata. This fundamental
difference exists because Strata allows for the layering and shear-wave velocity to vary (see Section 3.3.1 and therefore the required thickness of the
49

CHAPTER 4. USING STRATA

50

sublayers changes.
The maximum thickness (hmax ) of the sublayers of i-th velocity layer
is take as a fraction of the minimum wavelength to be captured by the
analysis:
vs,i
(4.1)
hmax,i = λfrac · λmin = λfrac ·
fmax
where λfrac is the wavelength fraction which typically varies between 1/10
and 1/5 (anything greater than 1/3 is not recommended), fmax is the maximum frequency of engineering interest which is typically around 20 Hz,
and Vs,i is the shear-wave velocity of the i-th layer. The actual thickness
of the sublayers is less than the maximum thickness such that the velocity
layer height divided by the sub-layer thickness is a whole number. These
parameters are defined on the General Settings tab. To prevent the layers from being auto-discretized the wavelength fraction can be increased
and the thickness velocity layers defined in the Soil Profile tab can be
reduced to an appropriate level.
In other site response programs, the location of the input motion or
the location of requested output (e.g., acceleration-time history) is generally referenced by a sublayer index. However, because the sublayers are
computed in Strata (and may change for each realization), the location is
defined in terms of the depth within the soil profile or at the top of the
bedrock. When the location is specified as Bedrock then the actual depth of
the location will change if the depth of the bedrock changes. The location
is specified with a drop down list shown in Figure 4.1 where the user can
specify the depth as Bedrock (Figure 4.1a) or a fixed depth (Figure 4.2c).

Figure 4.1: Location selection: (a) top of bedrock, (b) switching to a fixed
depth, and (c) fixed depth of 15.

CHAPTER 4. USING STRATA

4.1.2

51

Interaction with Tables

For the majority of Strata, the input behaves as you would expect. The one
exception may come when editing tables. Table cells are selected with one
click. Once a cell is selected, the cell can be edited by typing. A cell’s edit
mode can be directly entered by double clicking on the cell. In some case
double clicking on a cell will produce a widget to aid specifying input to
the cell.
All tables used in Strata are dynamic meaning the number of rows can
be changed. Rows are added to the bottom on the list with the Add button.
The Insert and Remove buttons are disabled until once a complete row has
been selected, which is most easily achieved by clicking on the number next
to the row of interest. Multiple continuous rows can be selected by pressing
the shift key while selecting the rows. After rows have been selected: Add
will add the same number of rows to the end of the table, Insert will
insert the same number of rows to prior to the currently selected rows, and
Remove will remove the selected rows. All rows in the table can be selected
by click on the button in the upper right portion of the table as shown in
Figure 4.2. Some table has cells that can not be edited and have a light
gray background. An example of these cells is in the Velocity Layers table
(shown Figure 4.2).
Data can copied from spreadsheets and be pasted into tables by:
1. Pressing Ctrl+v,
2. Clicking on the table, and then selecting Paste from the Edit menu,
or
3. Right clicking on the table and selecting paste.
The table will automatically increase the numbers of rows to accommodate
the size of the pasted data.

4.1.3

Non-Linear Curves

The nonlinear (shear-modulus reduction and damping) curves can be specified through three different methods in Strata: (1) fixed models that are
present by default and cannot be removed, (2) user defined curves that can
be used across projects, and (3) temporary models that only exist for the
project.

Figure 4.2: By clicking on the button circled in red all rows in the table are selected.

CHAPTER 4. USING STRATA
52

CHAPTER 4. USING STRATA

53

Fixed Non-Linear Models
The following shear-modulus reduction models are included by default:
• Darendeli (2001)
• Iwasaki et al. (1976) – 0.25 and 1.0 atm
• Seed and Idriss (1970) – Mean and Upper
• Vucetic and Dobry (1991) – PI = 0, 15, 50, 100, and 200
and the following damping curves are included by default:
• Darendeli (2001)
• Seed and Idriss (1970) – Lower and Mean
• Vucetic and Dobry (1991) – PI = 0, 15, 50, 100, and 200
User Defined Models
Non-linear curve models can be defined for use across multiple projects by
adding models to the library. The nonlinear property manager is opened
by selecting Add/Remove Non-Linear Property Curves from the Tools
menu. Using the dialog (Figure 4.3), a new model can be defined by
following these steps:
1. Click the Add button to add a new curve to with the normalized
shear-modulus reduction or damping models list.
2. Rename the model from “Untitled“ to something meaningful.
3. Add the data points to the curve.
A curve can be removed by selecting the curve and then clicking on
the Remove button. Models defined in this manner will be added to the
nonLinearCurves.strd file found in the Strata installation folder.
Temporary Models
If you want to define a curve without adding it to the library of models,
simply select Custom from the drop-down list. Changing to the Custom
model does not clear the previous models data which allows for a model
to be modified.

CHAPTER 4. USING STRATA

Figure 4.3: The nonlinear curve manager.

54

CHAPTER 4. USING STRATA

4.1.4

55

Recorded Motion Dialog

The Recorded Motion Dialog is used to load a recorded motion into Strata
and appears when the Add button is clicked in the Record Motion(s) table.
The dialog, shown in Figure 4.4 allows the user to load a variety of motions
in most formats using the following steps:
1. Click on the File... button and select the acceleration time series
file. If the file is from the NGA database, then the remainder of
the form will be automatically completed as shown in Figure 4.5.
Regardless of the file type, the file is read and loaded into the preview
area.
2. The remaining fields need to be filled to reflect the information in the
file. Information is required for all fields except for the Description
field. Fields can be completed by either typing values in, or selecting
from the file preview and dragging the selected text into the field. The
Start line and Stop line control which lines in the file contain the
data. A zero value for the Stop line will result in the data being read
until the end of the file. The file preview can be colored by clicking
on the Refresh button. The colors have the following meanings:
green text found prior to the acceleration-time series data (ignored).
blue acceleration-time series data
red text after the time series data (ignored).
An example of the colored data is shown in Figure 4.5.
3. The scale factor can be selected at this time or after the motion has
been loaded. The scale factor should result in the motion being in units
of gravity. After the motion has been loaded, the scale factor can also
be adjusted by setting the peak-ground acceleration.
4. After the form has been completed, the time-series can be viewed by
clicking on the Plot button.
5. Click OK to finish loading the file.

Figure 4.4: The initial view of the Recorded Motion Dialog.

CHAPTER 4. USING STRATA
56

Figure 4.5: An example a completed Recorded Motion Dialog.

CHAPTER 4. USING STRATA
57

CHAPTER 4. USING STRATA

4.1.5

58

Output Widget

Once the site response calculation has been completed the view will change
to the Output Widget. The Output Widget is also selected if a Strata output
file is opened, or by selecting Output View from the Window menu.
Strata does not output any data files automatically. The data can also
be exported by selecting Export... from the File menu. The exported
data file is in a comma-separated values (CSV) format that can be easily
opened with Excel or other spreadsheet program. All data (even disabled
results) are included in the files.
Each result is characterized by a site realization and a motion. If the
data is enabled, the check box with be checked, otherwise it is considered
disabled. Results that are enabled will be included when the statistics are
calculated. At the bottom of the table there are two buttons that allow
all motions or sites related to the currently selected result to be enabled or
disabled. In the example (Figure 4.6), both the site and motion are enabled.
Whenever the status of a result is changed (e.g. from disabled to enabled)
the Recompute Statistics button will become enabled allowing the user
to update the median and standard deviation.
The Output Widget shown in Figure 4.6 is used to examine the results
of a calculation. In the figure, the output is acceleration response spectrum
computed at the surface for 20 motions and 20 different site profiles. A
individual result can be selected by either clicking on the corresponding
in the Data Selection table, or by clicking on the result in the plot. In
both cases, the result is colored green if the result is enabled, or red if
the result is disabled. After the status of a result has been changed, the
Recomputed Statistics button will become enabled indicating that the
median and standard deviation (shown on the plot in solid and dash blue
lines, respectively) need to be updated.
The current plot can be printed by selecting Print... or Print to
PDF... from the File menu. The current plot can also be copied by right
clicking on the plot and selecting Copy.

4.2

Examples

The following examples give a basic introduction to using Strata to perform
equivalent linear site response analysis. The examples files are found

Figure 4.6: Using the Output view to examine the results of a calculation.

CHAPTER 4. USING STRATA
59

CHAPTER 4. USING STRATA

60

0

Depth (m)

20
40
60
80

Preferred
Minimum
Maximum

100
0

200

400

600

800

Shear-wave velocity (m/s)

Figure 4.7: The shear-wave velocity profile of the Sylmar County Hospital
Parking Lot site (Chang, 1996).
within the examples folder in the installation path. The examples can be
opened by either double clicking on their file, or selecting them from the
Open... item from the File menu.
All examples use the deep alluvium Sylmar County Hospital Parking
Lot (SCH) site located in Southern California for the site profile. The soil
types and velocity layering of the site was proposed by Chang (1996). The
soil properties are listed in Table 4.1 with a water table at a depth of 46
meters. The nonlinear properties for each of the layers were computed
using the Darendeli (2001) empirical model with PI=0, OCR=1, and the
confining pressures listed in Table 4.1. The corresponding velocity profile
is shown in Figure 4.7. The time-averaged shear-wave velocity over the
top 30 meters was computed as 273 m/s which classifies the site as a USGS
site class C.

4.2.1

Example 1: Basic time domain

In the first example, the site response will be computed for the Sylmar
County Hospital Parking Lot site using seven different recorded mo-

Soil Type
Vs (m/s)
γtot (kN/m3 ) σ0v (kPa)
Alluvium (Sand)
200 (150 to 230)
18
54
Alluvium (Sand)
300 (240 to 350)
18
222
Alluvium (Sand)
460 (370 to 550)
19
562
Alluvium and Older 700 (580 to 750)
22
776
Alluvium (Sand)
Bedrock
760
22

σ0m (atm)
0.36
2.2
5.6
7.7

Table 4.1: Soil profile at the Sylmar County Hospital Parking Lot site (Chang, 1996). The mean effective
stress (σ0m ) is computed assuming a k0 of 1/2 and a water table depth of 46 meters.

91+

Depth Range (m)
0 to 6
6 to 31
31 to 61
61 to 91

CHAPTER 4. USING STRATA
61

CHAPTER 4. USING STRATA

62

tions. This example can be opened from the example-1-td.stri file in
the examples directory. Strata has the option of saving the time-series data
within the project file which is useful when the project files is transferred
to different computers. For this example, the acceleration-times series are
included in the project file. While the different tabs in the project can be
accessed in any order, going from left-to-right is preferred as information
in one tab may depend on tabs to the left of it.
Defining soil types
A soil type is defined by a name, unit weight, an initial damping, shearmodulus reduction model, damping model, and notes about the soil type.
The name is only used to identify the soil type, and the notes are only
present to type maintain sanity. For soil types that use the Darendeli (2001)
model for either of nonlinear properties, addition properties required by
the empirical model must be defined.
For the site profile presented in Table 4.1, four different soil types are
required. Each soil type uses the Darendeli (2001) model for both the
shear-modulus reduction and damping, but the mean effective stress in
the model changes.
Defining soil profile
The soil profile combines the shear-wave velocity profile and soil types
together. For each soil layer in the profile, requires that a thickness, soil
type, and shear-wave velocity must be assigned. The soil types are selected
from a list of soil types defined in the Soil Types tab.
Defining the acceleration-time series
The example provides seven acceleration-time series scaled to a target
response spectrum and standard deviation. A motion is added by clicking
on the Add button in the Recorded Motion(s) table. After clicking add a
dialog will appear to assist in loading the file, information on this dialog is
presented in Section 4.1.4. Once the motion has been loaded, only the scale
factors can be changed by either specifying a peak-ground acceleration or
a scale factor. If the Save motion data with the input file check box
is checked, then the Strata input file will contain all of the information

CHAPTER 4. USING STRATA

63

regarding the time series. This increases the size of the input file, but
allows for a single file to contain all required input.
Defining the output
There are generally four different types of output that can be computed by
Strata:
Response location output involves the response at a given location (e.g.
acceleration response spectrum or a time series).
Ratio output involves the ratio of a response at two different location at a
given location (e.g. transfer function or spectral ratio).
Profile variation of a property with depth (e.g. maximum shear-strain).
The profiles are computed at particular depths because the layering
of the model may change as a result of variation.
Non-linear curves the nonlinear curves.
Note: Only the information requested is stored.
Computing
After the project is fully defined, switch the Compute tab and click on the
Compute button. Information regarding the calculation will be displayed
in the window and once the calculation is complete the view will switch
to the output widget. For more information on using the output widget
see Section 4.1.5. The input widget can be selected by selected Input View
from the Window menu, or by pressing F2.

4.2.2

Example 2: RVT and Site Variation

In the second example, the site response for the SCH site using random
vibration theory and variation of the shear-wave velocity profile. This
example can be opened from the example-2-rvt.stri file in the examples
directory. The site profile is exactly the same as was used in the previous
example. There are only two changes that will be discussed: the change
from time domain to random vibration theory, and how to enable site
variation. Both of these settings are enabled in the General Settings tab.

CHAPTER 4. USING STRATA

64

Figure 4.8 shows the proper settings for enabling RVT and shear-wave
velocity variation.

Figure 4.8: Settings to enable RVT site response and variation of the shearwave velocity.

Varying the shear-wave velocity
The parameters that control the shear-wave velocity variation are found on
the Soil Profile tab. By enabling variation of the site profile a number
of new columns will appear in the velocity layers table, as well as a new
group box containing a number of properties to vary. In this example, only
the shear-wave velocity will be varied using the USGS C parameters for
both the standard deviation and correlation models.

CHAPTER 4. USING STRATA

65

Defining an input acceleration response spectrum
Random vibration theory allows for the input motion to be defined using either a Fourier amplitude spectrum (FAS) or an acceleration response
spectrum. In this example, a response spectrum was computed using Abrahamson and Silva (1997) empirical model for a magnitude 7.0 earthquake
generated by a strike-slip fault with a distance of 20 km. The duration of
the event was computed uses the Abrahamson and Silva (1996) empirical
relationship for the duration using an normalized Arias intensity of 0.75.
Whenever a response spectrum is used for input it is important to check
that the resulting FAS looks appropriate. Strata allows the user to see both
the data and the plots of the data prior to the calculation.

Bibliography
N. A. Abrahamson and W. J. Silva. Empirical ground motion models.
Technical report, Brookhaven National Laboratory, Upton, New York,
1996.
N. A. Abrahamson and W. J. Silva. Empirical response spectral attenuation
relations for shallow crustal earthquakes. Seismological Research Letters,
68(1):94–127, January-Feburary 1997.
A. H.-S. Ang and W. H. Tang. Probability Concepts in Engineering Planning
and Design: Volume I–Basic Principles. John Wiley and Sons, New York,
1975.
A. H.-S. Ang and W. H. Tang. Probability Concepts in Engineering Planning
and Design: Volume II–Descision, Risk, and Reliability. Ang and Tang, 1990.
D. Boore. Simulation of ground motion using the stochastic method. Pure
and Applied Geophysics, 160(3-4):635–676, March 2003.
D. M. Boore and W. B. Joyner. A note on the use of random vibration theory
to predict peak amplitudes of transient signals. BSSA, 74(5):2035–2039,
October 1984.
J. N. Brune. Tectonic stress and the spectra of seismic shear waves from
earthquake. Journal of Geophysics Research, 75(26):4997–5009, 1970.
J. N. Brune. Correction. Journal of Geophysics Research, 76:5002, 1971.
K. W. Campbell. Prediction of strong ground motion using the hybrid
empirical method and its use in the development of ground-motion
(attenuation) relations in the Eastern North America. BSSA, 93(3):1012–
1033, June 2003.
66

BIBLIOGRAPHY

67

D. Cartwright and M. S. Longuet-Higgins. The statistical distribution of the
maxima of a random function. Proceedings of the Royal Soceity of London,
Series A, Mathematical and Physical Sciences, 237(1209):212–232, 1956.
S. W.-Y. Chang. Seismic response of deep stiff soil deposits. PhD thesis, University of California, Berkeley, December 1996.
C. Costantino, E. Heymsfield, and Y. Gu. Site specific estimates of surface
ground motions for the K-reactor site, Savannah River plant. Report
CEERC-91-003, Structural Analysis Division, Nuclear Energy Department, Brookhaven National Laboratory, Upton, NY, 1991.
M. B. Darendeli. Development of a new family of normalized modulus reduction
and material damping curves. PhD thesis, The University of Texas, Austin,
2001.
D. A. Gasparini and E. H. Vanmarcke. Simulated earthquake motions
compatible with prescribed response spectra. Technical Report R76-4,
Massachusetts Institute of Technology, Cambridge, MA, 1976.
B. O. Hardin and V. P. Drnevich. Shear modulus and damping in soils:
design equations and curves. Journal of soil mechanics and foundation
engineering division, ASCE, 98(SM7):667–692, June 1972.
I. M. Idriss and J. I. Sun. SHAKE91: a computer program for conducting
equivalent linear seismic response analyses of horizontally layered soil deposits.
Center for Geotechnical Modelling, Department of Civil and Environmentall Engineering, University of California, Davis, 1992.
T. Iwasaki, F. Tatsuoka, and Y. Takagi. Dynamic shear deformation properties of sand for wide strain range. Technical Report 1085, Civil Engineering Institue, Ministry of Construction, Tokyo, Japan, 1976.
S. L. Kramer. Geotechnical Earthquake Engineering. Prentice Hall, Upper
Saddle River, NJ, 1996.
R. McGuire, G. Toro, T. O’Hara, J. Jacobson, and W. Silva. Probablistic seismic hazard evaluations at nuclear plant sites in the Central and Eastern
United States: Resolution of the Charleston Earthquake issue. Technical Report NP-6395-D, Electric Power Research Institute, Palo Alto, CA,
April 1989.

BIBLIOGRAPHY

68

E. M. Rathje, A. R. Kottke, and C. M. Ozbey. Using inverse random vibration theory to develop input Fourier amplitude spectra for use in site
response. In 16th International Conference on Soil Mechanics and Geotechnical
Eningeering: TC4 Earthquake Geotechnical Engineering Satellite Conference,
pages 160–166, Osaka, Japan, September 2005.
P. B. Schnabel, J. Lysmer, and H. B. Seed. Shake: A computer program for
earthquake response analysis of horizontally-layered sites. Technical Report EERC-72/12, Earthquake Engineering Research Center, University
of California, Berkeley, CA, 1972.
J. F. Schneider, W. J. Silva, S. J. Chiou, and J. C. Stepp. Estimation of ground
motion at close distances using the band-limited-white-noise model. In
Fourth Interntional Conference on Seismic Zonation, volume 4, pages 187–
194, Stanford, CA, 1991. EERI.
H. B. Seed and I. M. Idriss. Soil moduli and damping factors for dynamic
response analyses. Technical Report EERC 70-10, Earthquake Engineering Research Center, University of California, Berkeley, CA, 1970.
W. J. Silva, N. Abrahamson, G. Toro, and C. Costantino. Description and
validation of the stochastic ground motion model. Final Report Contract
No. 770573, Brookhaven National Laboratory, Upton, N.Y., 1997.
G. Toro, W. Silva, R. McGuire, and R. Hermann. Probabilistic seismic
hazard mapping of the Mississippi Embayment. Seismological Research
Letters, 63(3):449–475, July-September 1992.
G. R. Toro. Probablistic models of site velocity profiles for generic and sitespecific ground-motion amplification studies. Technical Report 779574,
Brookhaven National Laboratory, Upton, New York, 1995.
E. Vanmarcke. Random fields, analysis and synthesis. The MIT Press, Cambridge, Massachusetts, 1983.
M. Vucetic and R. Dobry. Effect of soil plasticity on cyclic response. Journal
of the Geoechnical Engineering Division, ASCE, 117(1):89–107, 1991.

Index
Fast Fourier transform, 16
motion type
outcrop, 6
within, 7
Nyquist frequency, 16
random vibration theory
Parseval’s theorem, 19
peak factor, 19
response spectrum inversion, 21
shear modulus
complex, 3
maximum (small strain), 8
transfer function
acceleration, 5
strain, 8
wave equation, 2

69



Source Exif Data:
File Type                       : PDF
File Type Extension             : pdf
MIME Type                       : application/pdf
PDF Version                     : 1.5
Linearized                      : No
Page Count                      : 77
Page Mode                       : UseOutlines
Author                          : 
Title                           : 
Subject                         : 
Creator                         : LaTeX with hyperref package
Producer                        : pdfTeX-1.40.18
Create Date                     : 2018:01:22 08:31:26-08:00
Modify Date                     : 2018:01:22 08:31:26-08:00
Trapped                         : False
PTEX Fullbanner                 : This is pdfTeX, Version 3.14159265-2.6-1.40.18 (TeX Live 2017/Arch Linux) kpathsea version 6.2.3
EXIF Metadata provided by EXIF.tools

Navigation menu