# Dynamic site response: fundamental period and stiffness profile

## Why this matters

The flexible-block Newmark displacement models (BT07, BM17, BM19) evaluate spectral acceleration at a period proportional to the fundamental period of the sliding mass T_s. Estimating T_s therefore precedes any flexible-block analysis. The base stiffness parameters required for the estimation — the maximum shear modulus G_o and the shear-wave velocity v^o_S at the embankment base — are established here from the soil profile via Ishihara’s (1996) shear-modulus model and the Gazetas–Dakoulas (1985) inhomogeneous truncated shear-beam formulation.

This methodology was originally distributed as the `dsra` package and was merged into `newmark` 1.1.0 (see `NEWS.md`). The functions `getSiteProperties()`, `geSiteTable()`, `getCylinderRoots()`, `fitModel.Ts()`, `Vs30toSID()` and `SIDtoVs30()`, together with the reference datasets `CylinderRoots`, `ShearModelParameters`, `SiteClass` and `SiteTable`, all originate from `dsra`.

## Fundamental period of an inhomogeneous truncated shear beam

For an embankment whose shear stiffness increases with depth according to a power-law profile

\[ G(z) \;\approx\; G_o \left(\frac{z}{H_\text{max}}\right)^{m_o}, \]

the j-th modal period is

\[ T_s^{(j)} \;\simeq\; \frac{4\,\pi\,H_\text{max}}{a^{(j)}\,(2 - m_o)\,v_S^{o}}, \]

where

- H_max is the maximum height of the embankment;
- a^(j) is the j-th eigenvalue of the characteristic equation for the laterally constrained embankment under harmonic base excitation, depending jointly on m_o and the truncation ratio λ_o;
- v^o_S = √(G_o / ρ) is the shear-wave velocity at the base, with G_o the maximum shear modulus at depth z = H_max and ρ the bulk mass density.

For practical evaluation, the fundamental mode (j = 1) dominates and the expression reduces to a function of the base shear-wave velocity, the maximum height, and the two dimensionless geometric ratios m_o and λ_o.

The truncation ratio λ_o captures the geometric influence of the berm width B_o on the dynamic response:

\[ \lambda_o \;=\; \frac{B_o}{B_\text{max}} \;=\; \frac{B_o}{2\,H_\text{max}\,s + B_o}, \]

where s is the reciprocal of the average slope gradient (tan β = 1 / s). A larger λ_o lengthens the fundamental period.

The homogeneity ratio m_o reflects the power-law variation of the shear modulus with depth. A larger m_o implies a steeper stiffness gradient from surface to base and lengthens the fundamental period relative to the uniform-stiffness case (m_o = 0).

The eigenvalue a^(1) is computed from a precomputed table of characteristic roots derived from the Gazetas–Dakoulas (1985) characteristic equation:

    library(newmark)

    an <- getCylinderRoots(
      mo          = 0.45,                 # homogeneity ratio
      lo          = 0.20,                 # truncation ratio
      no          = 1,                    # mode number
      model       = "nlm",                # interpolation: "lm", "nlm", "dt", "rf"
      extrapolate = TRUE                  # warn (do not clamp) outside the table
    )

The reference table is the `CylinderRoots` dataset; its support in (m_o, λ_o) is m_o ∈ [0, 0.95] and λ_o ∈ [0, 0.495] across modes 1..8.

When the analyst already has a discrete V_S(z) profile (for instance from CPT-derived correlations or downhole logging), `fitModel.Ts()` computes T_s directly via Rayleigh’s method without going through the modal eigenvalue:

    Ts <- fitModel.Ts(
      VSm = c(220, 280, 340),               # shear-wave velocities by layer (m/s)
      hs  = c(10, 10, 10),                  # layer thicknesses (m)
      zm  = c(5, 15, 25)                    # layer mid-depths (m)
    )

## Ishihara stiffness model

The maximum (small-strain) shear modulus at depth z is estimated as a function of void ratio e_o and octahedral effective stress σ′_o(z) using Ishihara’s (1996) form:

\[ G(z,\,e_o) \;\approx\; A\,\frac{(C_e - 1)^{2}}{1 + e_o}\,\left(\frac{\sigma'_o(z)}{p_\text{ref}}\right)^{n}, \]

where

- A is a dimensionless material-dependent amplitude constant;
- C_e is a void-ratio shape parameter controlling the sensitivity of G to void ratio;
- p_ref is a reference pressure (typically atmospheric, 101.3 kPa);
- n is the stress exponent governing the power-law dependence of stiffness on confining pressure.

Material parameters (A, C_e, n) are tabulated in the package as `ShearModelParameters`, indexed by USCS soil class. Void-ratio ranges by USCS class are also packaged (used internally by the synthetic-profile generator):

    data(ShearModelParameters)            # 20 rows, parameters by ModelID
    data(SiteTable)                       # ~256k rows, sampled site profiles
    data(CylinderRoots)                   # ~307k rows, eigenvalues a(mo, lo, n)
    data(SiteClass)                       # 9 rows, ASCE 7-22 site-class table

The maximum shear modulus G_o corresponds to z = H_max, where the octahedral stress reaches its largest value under the self-weight of the fill. The base shear-wave velocity follows directly:

\[ v_S^{o} \;=\; \sqrt{\frac{G_o}{\rho}}. \]

## Synthetic profile generation

For preliminary studies where laboratory data are limited, the package generates Monte Carlo realisations of synthetic soil profiles consistent with a USCS classification. Each realisation samples void ratio, unit weight and plasticity from USCS-conditioned distributions, computes the stiffness profile G(z) via the Ishihara expression, fits the power-law form to extract (G_o, m_o, v^o_S), and computes T_s via the modal expression above.

    SiteProps <- getSiteProperties(
      Hs     = 30,                         # total height to hard ground (m)
      USCS   = c("GC", "CL", "ML"),        # USCS codes (top to bottom)
      POP    = 100,                        # pre-consolidation pressure (kPa)
      Water  = 0,                          # water-table depth as fraction of Hs
      NR     = 100,                        # Monte Carlo realisations
      h      = 1.00,                       # layer thickness (m)
      levels = c(0.05, 0.5, "mean", 0.95),
      Vref   = 760
    )
    # returns a data.table with one row per (Hs, POP, Hw, level) combination
    # columns: USCS, Go, mo, Ts, VSo, VS30, Z500, Z1000, level

`geSiteTable()` is the lower-level entry point that produces a single realisation; `getSiteProperties()` wraps it in a Monte Carlo loop and returns quantile summaries.

## Site-class utilities

The eight-class site identifier of ASCE 7-22 / NBCC is mapped bidirectionally with the equivalent V_S30:

    Vs30toSID(c(120, 360, 760, 1500))
    #> [1] "E"  "CD" "BC" "A"

    SIDtoVs30(c("BC", "C", "CD", "D"))
    #> [1] 760 540 370 255

The canonical class boundaries are:

Class | Vs30 range (m/s) | Representative Vs30
--- | --- | ---
A | ≥ 1500 | 1500
B | 910 – 1500 | 1200
BC | 640 – 910 | 760
C | 440 – 640 | 540
CD | 300 – 440 | 370
D | 210 – 300 | 255
DE | 150 – 210 | 180
E | < 150 | 150

## Site amplification models

Rock-to-site amplification of spectral ordinates uses the NGA-East ergodic site-response family, selected with the `models` argument of `fitSaF()` and `fitSaFMixture()`:

- `"ST20"` (default) — the final published generation: linear model and F760 impedance/gradient weighting per Stewart et al. (2020), nonlinear model per Hashash et al. (2020) with a period-dependent nonlinear reference velocity (760 m/s below 0.4 s, 3000 m/s above) and a PGA-dependent nonlinear sigma. Coefficients are regenerated from the USGS nshmp-lib machine-readable table.
- `"ST17"` — the 2017 PEER-report generation (PEER Report 2017/05), retained for backward comparison with earlier studies.

Both models are referenced to very hard rock (VS30 = 3000 m/s). With `vref = 760` (the usual reference-rock hazard condition) the reported mean factor is the exact difference of the composite model, ln F = F(VS30) − F(760), and the input PGA is converted to the 3000 m/s anchor with the PEER 2017/05 factor PGA3000 = PGA760/2.275. The dispersion σln F is evaluated at the target VS30 only: the reference-rock hazard already carries its own dispersion, so no reference-side sigma is added. At `vs30 == vref` the factor is deterministically 1 (zero dispersion).

When the rock hazard is available as an analytic lognormal mixture (one component per rupture-variant and ground-motion-model branch, as produced by deterministic scenario engines), `fitSaFMixture()` composes the mixture with the site-factor lognormal exactly — Gauss–Hermite quadrature over the PGA channel and monotone inversion of the pooled distribution — with no Monte Carlo sampling; results are reproducible byte-for-byte.

## Application

Most TSF and waste-rock dump assessments call `getSiteProperties()` once per (geometry, material) pair to obtain T_s and m_o, then pass T_s to the Newmark displacement workflow (`fitDnCurve`). The function returns a distribution because each input parameter (void ratio, unit weight, G profile coefficients) is sampled from USCS-conditioned ranges; the analyst chooses a quantile of T_s for the design.

For site amplification of the rock UHS to the target V_S30, see `?fitSaF` and the methodology in `vignette("ensemble-formulation", package = "newmark")`.

## References

- Gazetas, G. & Dakoulas, P. (1985). Soil Dyn. Earthq. Eng. 4(4):166–182.
- Hashash, Y.M.A. et al. (2017). Recommendation for Ergodic Nonlinear Site Amplification in Central and Eastern North America. PEER Report 2017/05.
- Hashash, Y.M.A. et al. (2020). Earthquake Spectra 36(1). doi:10.1177/8755293019878193.
- Ishihara, K. (1996). Soil Behaviour in Earthquake Geotechnics. Oxford University Press.
- Stewart, J.P. et al. (2020). Earthquake Spectra 36(1). doi:10.1177/8755293019878185.
- ASCE/SEI 7-22 (2022). Minimum Design Loads and Associated Criteria for Buildings and Other Structures. American Society of Civil Engineers.