Skip to contents

Newmark ensemble formulation

Source: vignettes/ensemble-formulation.Rmd
ensemble-formulation.Rmd

Overview

newmark answers one question: given a site-specific seismic hazard and a slope with known dynamic properties, what horizontal seismic coefficient k_max is required to limit permanent co-seismic displacement to a specified target d*?

The performance-based design problem is formulated as: for a given target tolerable displacement d* > 0, target return period T_R, and target exceedance probability p ∈ (0, 1), choose the smallest yield coefficient k_y such that the displacement induced at T_R remains below d* with probability at least 1 − p:

kmax(p∣d*,TR)=inf⁡{ky:P[DN(ky;TR)>d*]≤p}. k_\text{max}(p \mid d^*, T_R) \;=\; \inf\!\Bigl\{\,k_y \,:\, P\!\bigl[\,D_N(k_y; T_R) > d^*\,\bigr] \,\leq\, p\,\Bigr\}.

The right-hand side has no closed form — neither the joint distribution of intensity measures nor the convolution with the empirical regression residual admits one. The framework evaluates the inverse mapping numerically via Monte Carlo, propagating two sources of uncertainty:

  1. Hazard uncertainty — spectral accelerations are drawn from the quantile curves of the uniform-hazard spectrum via a Gaussian copula that preserves the full inter-period correlation matrix (Baker & Jayaram 2008); no star/hub simplification.
  2. Model uncertainty — six empirical Newmark models are combined in a logic-tree ensemble with user-assigned epistemic weights w_i, ∑w_i = 1.

For the full mathematical derivation and calibration discussion, see:

Verri Kozlowski, A. (2026). Probabilistic estimation of Newmark displacements and seismic coefficients under hazard uncertainty. Working paper.

Three-stage structure

Stage 1 — Site amplification (fitSaF)

Rock-level spectral ordinates are amplified to the target site Vs30 using the Seyhan & Stewart (2014) nonlinear site-factor model (ST17).

  • Mode A (mean-only input): rock Sa is deterministic; dispersion comes from ST17 sigma alone.
  • Mode B (quantile input): Sa is sampled jointly across periods via Gaussian copula; dispersion reflects combined hazard and site-factor uncertainty.

See ?fitSaF for the function-level interface.

Stage 2 — Displacement curve (fitDnCurve)

For each Monte Carlo realisation, one ground-motion scenario is drawn (a consistent tuple of PGA, Sa(1.3Ts), Sa(1.5Ts)) and passed through all active displacement models simultaneously using a shared aleatory residual. The result is a family of coherent displacement curves D_N(k_y).

Stage 3 — Coefficient inversion (invertDnDraws)

Each curve is projected monotone and inverted in log-log space to find the yield acceleration k_max such that D_N(k_max) = d*. Mean and quantiles of k_max over all realisations constitute the output.

Displacement models

ID Authors Type Spectral reference
AM88 Ambraseys & Menu (1988) rigid block PGA
JB07 Jibson (2007) rigid block PGA, Arias intensity
SR08 Saygili & Rathje (2008) rigid block PGA, Arias intensity
BT07 Bray & Travasarou (2007) flexible block Sa(1.5 Ts)
BM17 Bray, Macedo & Travasarou (2018) flexible block, subduction Sa(1.5 Ts)
BM19 Bray & Macedo (2019, corr. 2023) flexible block, crustal Sa(1.3 Ts)

Activate/deactivate models via the weights argument (0 = inactive).

For BM19, near-fault pulse motions are flagged when PGV > 115 cm s⁻¹ and the equation switches to the Bray-Macedo (2023) D100 / D50 form controlled by the NFC argument ("D100" is the maximum-component default, conservative for slopes within ±45° of fault-normal; "D50" is the median-component case for other orientations). A sub-regime split at PGV = 150 cm s⁻¹ inside the pulse equation captures the saturation of seismic displacement at very high PGV.

Key assumptions

  • Shared aleatory residual across models per realisation: a single z^(n) ∼ 𝒩(0,1) is drawn and applied to all six models scaled by each model’s own σ_lnD. This corresponds to ρ = 1 cross-model residual correlation. Independent residuals per model (ρ = 0) underestimates the joint ensemble variability; intermediate cross-model correlations are an open refinement direction.
  • Independent site-factor and Sa uncertainty within each realisation.
  • Log-log linear inversion with boundary extrapolation. Per-model calibration ranges (getKyLimits()) are reported as a diagnostic but not enforced as a clamp inside invertDnDraws().

Operational defaults

  • Monte Carlo sample size N_S: routine practice is N_S ∈ 10³,10⁴10³, 10⁴. Convergence diagnostics (band-width stability under N_S × 2) accompany each application.
  • k_y grid points: N_k = 30, log-spaced over 0.01,max(PGA,0.80)0.01, max(PGA, 0.80) g (default of getDnKy()).
  • Reported quantiles: mean plus p ∈ {0.16, 0.50, 0.84}. Extreme-tail reporting (p ∈ {0.05, 0.95}) is requested explicitly via the p argument to invertDnDraws().

References

  • Baker & Jayaram (2008). Earthquake Spectra 24(1):299–317.
  • Bray & Travasarou (2007). J. Geotech. Geoenviron. Eng. 133(4):381–392.
  • Bray, Macedo & Travasarou (2018). J. Geotech. Geoenviron. Eng. 144(3):04017124.
  • Bray & Macedo (2019, corr. 2023). J. Geotech. Geoenviron. Eng. 145(12); Soil Dyn. Earthq. Eng. 168:107835.
  • Jibson (2007). Engineering Geology 91(2–4):209–218.
  • Saygili & Rathje (2008). J. Geotech. Geoenviron. Eng. 134(6):790–803.
  • Seyhan & Stewart (2014). Earthquake Spectra 30(3):1241–1256.