Skip to contents

Getting started with newmark

newmark-quickstart.Rmd

Overview

newmark implements a four-workflow pipeline for probabilistic seismic performance analysis of slopes and embankments:

  1. Dynamic Site Response — compute the fundamental period Ts and site-amplified UHS from the soil profile (getSiteProperties, getCylinderRoots, fitSaF).
  2. Hazard import — import PSHA output from OpenQuake (buildGMDP).
  3. Displacement curves — Monte Carlo Newmark ensemble Dn(ky) (getDnKy, fitDnCurve).
  4. Seismic coefficient — invert ensemble draws to kmax(d*) (invertDnDraws).

This vignette demonstrates Workflows 3 and 4 using the bundled example dataset. For the full pipeline see the pipeline vignette (vignette("pipeline", package = "newmark")).

Example dataset

uhs.csv is a site-amplified uniform-hazard spectrum included with the package: NBCC hazard model, Vs30 = 560 m/s, TR = 10 000 yr, with quantile levels 0.05–0.95 and mean. It exercises the distinct quantile-table contract. The canonical OpenQuake PSHA path instead calibrates an analytic distribution with fitSaLognormal() and discards the raw OpenQuake fractiles.

library(newmark)
library(data.table)
#> 
#> Attaching package: 'data.table'
#> The following object is masked from 'package:base':
#> 
#>     %notin%

uhs <- fread(system.file("extdata", "uhs.csv", package = "newmark"))
uhs[p %in% c("0.16", "mean", "0.84") & Tn <= 0.2, .(Tn, p, Sa)]
#>        Tn      p      Sa
#>     <num> <char>   <num>
#>  1: 0.000   0.16 0.12513
#>  2: 0.000   0.84 0.21655
#>  3: 0.000   mean 0.18228
#>  4: 0.050   0.16 0.39859
#>  5: 0.050   0.84 0.69551
#>  6: 0.050   mean 0.54639
#>  7: 0.075   0.16 0.43086
#>  8: 0.075   0.84 0.73093
#>  9: 0.075   mean 0.57871
#> 10: 0.100   0.16 0.41504
#> 11: 0.100   0.84 0.72926
#> 12: 0.100   mean 0.56767
#> 13: 0.120   0.16 0.39547
#> 14: 0.120   0.84 0.67237
#> 15: 0.120   mean 0.53474
#> 16: 0.150   0.16 0.36396
#> 17: 0.150   0.84 0.61739
#> 18: 0.150   mean 0.49354
#> 19: 0.180   0.16 0.32485
#> 20: 0.180   0.84 0.54772
#> 21: 0.180   mean 0.44060
#> 22: 0.200   0.16 0.29781
#> 23: 0.200   0.84 0.49987
#> 24: 0.200   mean 0.41868
#>        Tn      p      Sa
#>     <num> <char>   <num>

Parameters 

# Ts: fundamental period of the sliding mass (s).
# In production, derived from getSiteProperties() + getCylinderRoots()
# using the soil USCS profile and slope geometry (Ishihara 1996,
# Gazetas & Dakoulas 1985). Here set to a representative value.
Ts <- 0.60

# Mw: moment magnitude from PSHA disaggregation.
Mw <- 6.8

# Ensemble weights (0 = model inactive).
weights <- c(AM88 = 1, JB07 = 0, BT07 = 1, SR08 = 1, BM17 = 0, BM19 = 1)

# Displacement targets (cm).
Da <- c(0.5, 2.5, 5.0, 25.0)

Workflow 3 — Displacement curve Dn(ky) 

ky     <- getDnKy(uhs = uhs, Ts = Ts)
result <- fitDnCurve(
  uhs     = uhs,
  ky      = ky,
  Ts      = Ts,
  Mw      = Mw,
  NS      = 200,
  weights = weights
)

result$curve[IDn == "ensemble" & p == "mean", .(ky, Dn)]
#>             ky           Dn
#>          <num>        <num>
#>  1: 0.01000000 86.112909559
#>  2: 0.01163118 73.753716887
#>  3: 0.01352844 62.795587132
#>  4: 0.01573517 53.094674718
#>  5: 0.01830186 44.526771614
#>  6: 0.02128722 36.984578509
#>  7: 0.02475955 30.375921204
#>  8: 0.02879828 24.622141350
#>  9: 0.03349580 19.655972385
#> 10: 0.03895957 15.418501983
#> 11: 0.04531458 11.855299242
#> 12: 0.05270620  8.912317146
#> 13: 0.06130353  6.532554143
#> 14: 0.07130325  4.654469890
#> 15: 0.08293409  3.212596456
#> 16: 0.09646214  2.139639920
#> 17: 0.11219685  1.369480074
#> 18: 0.13049819  0.841540406
#> 19: 0.15178479  0.499554323
#> 20: 0.17654363  0.289467154
#> 21: 0.20534109  0.166422447
#> 22: 0.23883592  0.099282452
#> 23: 0.27779437  0.063764255
#> 24: 0.32310764  0.042335893
#> 25: 0.37581233  0.027899861
#> 26: 0.43711410  0.018198085
#> 27: 0.50841530  0.011749885
#> 28: 0.59134701  0.007510662
#> 29: 0.68780637  0.004753403
#> 30: 0.80000000  0.002978862
#>             ky           Dn
#>          <num>        <num>

Workflow 4 — Seismic coefficient kmax(d*) 

kmax <- invertDnDraws(
  draws = result$draws,
  Da = Da,
  weights = weights
)
kmax[p %in% c("0.16", "mean", "0.84")]
#>        Da      p       kmax beyond
#>     <num> <char>      <num>  <num>
#>  1:   0.5   mean 0.13336054      0
#>  2:   0.5   0.16 0.08537750      0
#>  3:   0.5   0.84 0.18307529      0
#>  4:   2.5   mean 0.07778503      0
#>  5:   2.5   0.16 0.04655504      0
#>  6:   2.5   0.84 0.10946700      0
#>  7:   5.0   mean 0.05869568      0
#>  8:   5.0   0.16 0.03183874      0
#>  9:   5.0   0.84 0.08488157      0
#> 10:  25.0   mean 0.02303988      0
#> 11:  25.0   0.16 0.00000000      0
#> 12:  25.0   0.84 0.03993320      0

kmax is in g. The normalised pseudostatic coefficient is Kh = kmax / PGA_rock x 100 %. Model root distributions are aggregated with their epistemic weights. result$draws$PGA.s retains the sampled PGA in g for each realization, including AM88 draws whose roots lie below the stored yield-acceleration grid.

Next steps