4.3. Regional Site Response Assessment

4.3.1. Problem Statement

Site response analysis is commonly performed to analyze the propagation of seismic wave through soil. As shown in Fig. 2.3.5.1, one-dimensional response analyses, as a simplified method, assume that all boundaries are horizontal and that the response of a soil deposit is predominately caused by SH-waves propagating vertically from the underlying bedrock. Ground surface response is usually the major output from these analyses, together with profile plots such as peak horizontal acceleration along the soil profile. When liquefiable soils are presenting, maximum shear strain and excess pore pressure ratio plots are also important.

../../../../../_images/siteResponse.png

Fig. 4.3.1.1 Simplified 1D site response analysis (courtesy of Pedro Arduino)

This example demonstrates a small-scale regional earthquake risk assessment which performs site response simulation for a group of 20 sites located in the San Francisco Bay area. The sites are distributed in space in a 4x5 grid, within a 3x3 grid of event sites. At each event site, 5 ground motion records of similar intensity are assigned.

../../../../../_images/example.png

Fig. 4.3.1.2 Location of testbed sites

4.3.2. Methodology

In this example, the coordinates of a list of sites are known and peak ground accelerations at these sites are of interest. A summary of the RDT procedure for the site response analysis is summarized here:

  1. Assemble required information for interested sites. These data include site locations, input motions, and site geological data, e.g., Vs30 (the time-averaged shear-wave velocity to 30 m depth) and depth to bedrock

  2. If the input motions are not provided, use OpenSHA to obtain proper input motions

  3. If the Vs30 data are not known, use SimCenter EQHazard to obtain proper Vs30 data

  4. If the depth to bedrock data are not known, they can be obtained from SoilGrid250 database

  5. Use Sediment Velocity Model (SVM) to generate soil shear wave velocity profile

  6. Generate material input for site response analysis based on predefined model calibration methods

  7. Perform site response analysis using RDT and obtain acceleration response at the surface of soil columns defined for each site

In this example, a total of twenty sites and their input motions have been pre-selected. Details of selected sits are summarized in tab_inputs. The directory structure for inputs is:

srt_example
├── rWHALE_config_SRT.json
└── input_data
    ├── model
        ├── getVs30.py
        ├── Scenario.json
        ├── FreeField3D_Dry.tcl
    ├── records
        ├── EventGrid.csv
        ├── RSN30.json
        ├── RSN63.json
        .
        .
        .
        ├── site0.csv
        ├── site1.csv
        .
        .
        .
        └── site8.csv
    └── input_params.csv

4.3.2.1. Obtaining Vs30 Data

Vs30 data can be obtained using EQHazard tool. EQHazard is a console application developed by the SimCenter to facilitate performing seismic hazard analysis(SHA) for regional earthquake simulations. The application makes use of OpenSHA to perform deterministic SHA. This application is developed using Java to facilitate the use of OpenSHA Java library. It is worth noting that this tool is most suitable for research applications where SHA needs to be performed for many sites given a single earthquake scenario, such as using a grid of sites for regional simulations. Several site data sources have been included in EQHazard, such as CGS/Wills VS30 Map (2015), Global Vs30 from Topographic Slope (Wald & Allen 2008), and more.

Using EQHazard is simple, the EQHazard tool is distributed as a jar file and can executed directly using:

java -jar EQHazard Scenario.json Output.json

Scenario.json is the name of the input file and output.json is the name to be used for the resulting output file.

Note

Java is required to run EQHazard

Python script getVs30.py provides a wrapper for EQHazard. It uses the site coordinates from input_params.csv and a template scenario file, Scenario.json, to obtain Vs30 data. An updated input file, input_params_srt.csv, is then created with Vs30 data populated as shown in Table 4.3.2.1.1. In this case, CGS/Wills VS30 Map (2015) is used.

Table 4.3.2.1.1 Input Parameters with Vs30 Data Obtained Using EQHazard

id

longitude

latitude

stories

yearbuilt

occupancy

structure

areafootprint

replacementCost

W

f_yield

T1

Vs30

DepthToRock

0

-122.182

37.422

0

1957

RES1

W1

1293.253476

1

55.25

118.75

0.29375

386.6

7.61

1

-122.18575

37.422

0

1927

RES1

W1

2568.835761

1

86.75

148.75

0.15625

386.6

8.79

2

-122.1895

37.422

0

2020

RES1

W1

1127.07738

1

114.75

123.75

0.39375

351.9

9.28

3

-122.19325

37.422

0

1951

RES1

W1

2085.809768

1

107.75

126.25

0.24375

386.6

8.65

4

-122.197

37.422

0

1960

RES1

W1

1481.567067

1

76.25

101.25

0.30625

386.6

7.74

5

-122.182

37.42366667

0

1986

RES1

W1

1240.719532

1

90.25

116.25

0.21875

386.6

7.65

6

-122.18575

37.42366667

0

1984

RES1

W1

2437.970509

1

93.75

141.25

0.35625

386.6

9.35

7

-122.1895

37.42366667

0

1983

RES1

W1

2668.60615

1

111.25

133.75

0.36875

351.9

9.88

8

-122.19325

37.42366667

0

1957

RES1

W1

1553.47035

1

104.25

108.75

0.19375

386.6

9.71

9

-122.197

37.42366667

0

1984

RES1

W1

1625.164708

1

58.75

106.25

0.23125

386.6

9.32

10

-122.182

37.42533333

0

2018

RES1

W1

1787.865419

1

62.25

143.75

0.34375

386.6

7.53

11

-122.18575

37.42533333

0

2018

RES1

W1

1821.174819

1

65.75

138.75

0.20625

386.6

9.02

12

-122.1895

37.42533333

0

1930

RES1

W1

1624.532488

1

69.25

146.25

0.31875

386.6

9.79

13

-122.19325

37.42533333

0

1922

RES1

W1

1134.673241

1

79.75

103.75

0.25625

351.9

9.84

14

-122.197

37.42533333

0

1947

RES1

W1

1625.097354

1

51.75

128.75

0.16875

386.6

10.67

15

-122.182

37.427

0

1969

RES1

W1

2289.96328

1

83.25

131.25

0.33125

386.6

7.66

16

-122.18575

37.427

0

2006

RES1

W1

2288.360712

1

100.75

136.25

0.38125

386.6

9.77

17

-122.1895

37.427

0

1959

RES1

W1

2291.605647

1

97.25

111.25

0.18125

351.9

9.8

18

-122.19325

37.427

0

1979

RES1

W1

2770.811412

1

72.75

113.75

0.26875

351.9

10.03

19

-122.197

37.427

0

1976

RES1

W1

2088.450773

1

118.25

121.25

0.28125

386.6

9.68

4.3.2.2. Obtaining Depth to Bedrock

The depth to bedrock can be obtained using several datasets. For example, SoilGrid250 map provides absolute depth to bedrock (in cm) predicted using the global compilation of soil ground observations. Accuracy assessment of the maps is available in Hengl [] and can be downloaded here. In this dataset, the depth to bedrock is considered as “depth (in cm) from the ground surface to the contact with coherent (continuous) bedrock.” It should be noted that the definition of bedrock in this database is based on lithology rather than shear wave velocity. Hence, depending on the geological condition, bedrock may imply different shear wave velocity and it might vary dramatically, e.g., western U.S. vs. eastern U.S.. Assumption is required for modeling purposes. In this example, the shear wave velocity of bedrock is assumed to be 760 m/s, as the sites are located in western U.S..

../../../../../_images/US_Depth_To_Rock.png

Fig. 4.3.2.2.1 Depth to bedrock (Cropped from SoilGrid250 to show continental U.S., 10 km resolution)

Example python script to obtain the data is shown below:

import rioxarray as rxr
import pandas as pd
from pyproj import Transformer

soil250 = rxr.open_rasterio('BDTICM_M_250m_ll.tif')
inputs = pd.read_csv('input_params_srt.csv')
for index, row in inputs.iterrows():
   lon = row['longitude']
   lat = row['latitude']
   transformer = Transformer.from_crs("EPSG:4326", soil250.rio.crs, always_xy=True)
   xx, yy = transformer.transform(lon, lat)
   depth_To_Rock = soil250.sel(x=xx, y=yy, method="nearest").values[0]
   print('depth to rock is {:.2f} m'.format(depth_To_Rock/100))

4.3.2.3. Generating Soil Profile Using SVM

Once the depth to rock is obtained, the next step is to generate soil profile based on Vs30. Amplitude, frequency content and duration of ground shaking are important parameters in ground motion modeling and are heavily dependent on crustal material properties. The Sediment Velocity Model (SVM) “translates VS30 and a proxy which describes stiffness of the near surface sediments into a 1D velocity profile” (Shi and Asimaki, 2018, []).

Datasets of soil parameters only exist for specific locations and there is no way to collect adequate data on all regions. For this reason, there is a need for a method which can determine these crustal material properties such as wave velocity and mass density for those regions which we do not have physically collected data. One such method, called the Community Velocity Models (CVMs), which was developed by the Southern California Earthquake Center (SCEC), is used to provide 3D seismic velocity information for Southern California. However, CVMs have their limitations. There exist two sub-models of the CVM which are the CVM-SCEC (CVM-S) and CVM-Harvard (CVM-H). Both sub-models rely on the spatial resolution of the available datasets and thus result in crude velocity profiles. Moreover, the CVM-H has an optional geotechnical layer (GTL) option which replaces the soft sediments in the top 350 m of soil, smoothing the velocity profile and obscuring basin edges which cause strong velocity impedance.

The SVM method was developed to get more accurate soil profiles as a function of VS30 which can be easily obtained from maps. It further allows for determination of basin depths based on crustal velocity models and preserves basin edges.

The Shi and Asimaki [] Sediment Velocity Model (SVM) is used in this example because of its ability to predict a soil velocity profile based solely on VS30. From this profile we are able to determine ground motion wavelengths and set a constraint for a minimum number of points needed to fully describe a soil column. Gmax at depth is also obtained from this shear wave velocity profile. The results of the test case are presented in Fig. 4.3.2.3.1.

../../../../../_images/Vs.png

Fig. 4.3.2.3.1 Generated shear wave velocity profile based on SVM

4.3.2.4. Including Soil Nonlinearity

By default, elastic isotropic material is assigned to the soil column due to limited available information. However, Nonlinearity can be incorporated into workflow by switching the constitutive model to other nonlinear model, for example, the multiaxial cyclic plasticity model proposed by Borja and Amies []. This simple total stress-based bounding surface plasticity model for clays was developed to accommodate multiaxial stress reversals. The model was constructed based on the idea of a vanishing elastic region undergoing pure translation inside a bounding surface, and an interpolation function for hardening modulus which varies with stress distance of the elastic region from the unloading point. The model can be applied to cohesive soils undergoing undrained stress reversals and cyclic loading. Fig. 4.3.2.4.1 illustrates behaviors of this model under cyclic loading. More information regarding the implementation of this model in OpenSees can be found here.

../../../../../_images/CTC_Fig3_Comparison.png

Fig. 4.3.2.4.1 Example behavior of Borja and Amies model

Table 4.3.2.4.1 Input for elastic material

id

longitude

latitude

Model

0

-122.182

37.422

Elastic

The example input parameters for using Borja and Amies are summarized in the following table:

Table 4.3.2.4.2 Input for multiaxial cyclic plasticity model for clays (Borja and Amies, [])

id

longitude

latitude

Model

Su_rat

Den

h/G

m

h0

chi

0

-122.182

37.422

BA

0.26

2

1

1

0.2

0

where Su_rat is the ratio between undrain shear strength and initial vertical stress.

4.3.2.5. Adding Surface Structural Model

A structure model can also be added to the workflow. In this case, the structure model uses the OpenSeesPyInput modeling application as discussed in Earthquake Assessment. The buildings are modeled as elastic-perfectly plastic single-degree-of-freedom (SDOF) systems defined by three input model parameters: the weight W, yield strength f_yield, and fundamental period T1. Functions are included which record the peak response as EDPs for each of the EDP types specified in the EDP_specs.json file.

"Modeling": {
"Application": "OpenSeesInputSRT",
"ApplicationData": {
   "mainScript": "FreeField3D_Dry.tcl",
   "modelPath": "model/",
   "ndm": 3,
   "includeStructureModel": true,
   "SimulationApplication": "OpenSeesPy-Simulation",
   "structureModel": "cantilever.py"
}

4.3.2.6. Performing Analysis

Configurations for RDT workflow are defined in rWHALE_config_SRT.json. In this example, the application for Modeling and Simulation is OpenSeesInputSRT and OpenSees-Simulation-SRT, respectively.

For each site, two additional input files for OpenSees are created during analysis: freefield_config.tcl and freefield_material.tcl. freefield_config.tcl defines the discretization of the profile and applied input motion(s).

 1# site response configuration file
 2set soilThick 99.0
 3set numLayers 198
 4# layer thickness - bottom to top
 5set layerThick(1) 0.50
 6set nElemY(1) 1
 7set sElemY(1) 0.500
 8set layerThick(2) 0.50
 9set nElemY(2) 1
10set sElemY(2) 0.500
11
12# ...
13# ...
14
15# motion file (used if the input arguments do not include motion)
16set accFile  xInput.acc
17set dispFile xInput.disp
18set velFile  xInput.vel
19set timeFile xInput.time
20set numEvt 2
21set accFile2  yInput.acc
22set dispFile2 yInput.disp
23set velFile2  yInput.vel
24set rockVs 760.0
25set omega1 8.59
26set omega2 42.95

freefield_material.tcl defines material inputs based on predefined calibration procedure based on Vs profile created from previous step. Input motion(s) is converted from records for each evaluation, for example, xInput.vel and xInput.time defiles velocity and time series, respectively, for one component of the motion.

 1set rho(1) 2.0
 2set shearG(1) 112692.20
 3set nu(1) 0.30
 4set E(1) 292999.71
 5
 6set mat(1) "ElasticIsotropic 1 $E(1) $nu(1) $rho(1)"
 7
 8
 9set rho(2) 2.0
10set shearG(2) 112692.20
11set nu(2) 0.30
12set E(2) 292999.71
13
14set mat(2) "ElasticIsotropic 2 $E(2) $nu(2) $rho(2)"
15
16# ... 
17# ...

The working/results directory structure is:

results
   └── 0
      ├── templatedir
         ├── bim.j
         ├── edb.j
         |── evt.j
         ├── sam.j
         ├── workflow_driver.bat
      ├── workdir.1
         ├── FreeField3D_Dry.tcl
         ├── freefield_config.tcl
         ├── freefield_material.tcl
         ├── xInput.vel
         ├── xInput.time
         ├── yInput.vel
         ├── yInput.time
         .
         .
         .
      ├── workdir.2
      .
      .
      ├── workdir.5
      └── response.csv
   └── 1
   .
   .
   └── 19

EDP_1-19.csv: reports statistics on the EDP results from simulating 5 ground motions for each building asset. The statistics reported are the median and lognormal standard deviation of peak ground acceleration (shown as PFA) in two directions.

Table 4.3.2.6.1 EDP Results

type

PGA

PGA

PGA

PGA

loc

0

0

0

0

dir

1

1

2

2

stat

median

beta

median

beta

0

0.6628029999999999

0.13600449843043236

0.736927

0.1275636943053479

1

0.633193

0.20304404095647471

0.717512

0.24904764429936244

2

0.721433

0.38880051595831655

0.715938

0.3514059337674102

3

0.713946

0.14054739273872968

0.764877

0.3915961386798949

4

0.743542

0.2171580362193277

0.49029799999999996

0.2407229856374786

5

0.6505810000000001

0.15337053057961506

0.56929

0.19106712482954594

6

0.667912

0.0821461264634056

0.76647

0.2854761912873225

7

0.8115100000000001

0.3580416351456673

0.7321770000000001

0.10374091963523832

8

0.5919449999999999

0.484515336289696

0.767825

0.2795311183002603

9

0.8206610000000001

0.169420903616953

0.679726

0.31719489352530283

10

0.6627770000000001

0.1160975874275244

0.623636

0.1383857442038908

11

0.664215

0.3068794966334912

0.681191

0.17780768288420262

12

0.7623260000000001

0.3342121915707141

0.709516

0.10729744618054217

13

0.540041

0.2794970883918737

0.702222

0.07850396436049084

14

0.628897

0.3682471518218779

0.6626449999999999

0.08292658813669537

15

0.63932

0.3279279372030125

0.574796

0.2115260225991987

16

0.74392

0.3603316789411304

0.690448

0.21550453311227016

17

0.8111510000000001

0.2871611406223856

0.714084

0.20898183490987293

18

0.6592899999999999

0.3766289667049084

0.742772

0.11470562563412998

19

0.67061

0.17423255167006754

0.740329

0.3700796708017528

Under Windows system, the commands can be summarized as:

4.3.3. Appendix

Configurations for RDT workflow:

 1{
 2  "Name": "rWHALE_",
 3  "Author": "Adam Zsarnóczay",
 4  "WorkflowType": "Parametric Study",
 5  "buildingFile":"buildings.json",
 6  "runDir": "...",
 7  "units": {
 8      "force": "kips",
 9      "length": "in",
10      "time": "sec"
11   },
12   "outputs": {
13    "EDP": true,
14    "DM": false,
15    "DV": false,
16    "every_realization": false
17  },
18  "Applications": {
19    "Building": {
20      "Application": "CSV_to_BIM",
21      "ApplicationData": {
22        "Min": "0",
23        "Max": "19",
24        "buildingSourceFile":"input_params_srt.csv"
25      }
26    },
27    "RegionalMapping": {
28      "Application": "NearestNeighborEvents",
29      "ApplicationData": {
30        "filenameEVENTgrid": "records/EventGrid.csv",
31        "samples": 5,
32        "neighbors": 4
33      }
34    },
35    "Events": [{
36      "EventClassification": "Earthquake",
37      "Application": "SimCenterEvent",
38      "ApplicationData": {
39        "pathEventData": "records/"
40      }
41    }],
42    "Modeling": {
43      "Application": "OpenSeesInputSRT",
44      "ApplicationData": {
45        "mainScript": "FreeField3D_Dry.tcl",
46        "modelPath": "model/",
47        "ndm": 3
48      }
49    },
50    "EDP": {
51      "Application": "UserDefinedEDP_R",
52      "ApplicationData": {
53        "EDPspecs": "EDP_specs.json"
54      }
55    },
56    "Simulation": {
57       "Application": "OpenSees-Simulation-SRT",
58       "ApplicationData": {
59       }
60    },
61    "UQ": {
62       "Application": "Dakota-UQ",
63       "ApplicationData": {
64          "method": "LHS",
65          "samples": 5,
66          "type": "UQ",
67          "concurrency": 5,
68          "keepSamples": true
69       }
70    }
71    },
72	"runType": "runningLocal",
73	"runDir": "path_to_runDir",
74	"localAppDir": "path_to_SimCenterBackendApplications"
75}