Skip to contents

Newmark ensemble formulation

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 analytic finite-mixture marginals, or from quantile curves when the input distribution is defined only in that form, through 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, fitSaFMixture)

Rock-level spectral ordinates are amplified to the target site Vs30 using the NGA-East ergodic site-response models: "ST20" (default; Stewart et al. 2020 linear model, Hashash et al. 2020 nonlinear model) or "ST17" (the 2017 PEER-report generation), selected via the models argument.

  • OpenQuake PSHA input: fitSaLognormal() preserves every reported mean and uses the seven raw fractiles only to fit sigmaLn(Tn). The raw fractiles do not survive calibration. After the consumer assigns the component key and unit weight, fitSaFMixture() transforms that analytic component through site response.
  • Analytic mixture input (fitSaFMixture): each rock component is conditioned on its same-realization PGA and composed with the site-factor lognormal by Gauss–Hermite quadrature. The output is another analytic finite mixture; no report fractiles or quantile spline are materialized.
  • Quantile input (fitSaF): retained for spectra whose distribution is already defined by quantiles, including the current MCE envelope workflow. It is not an alternate route for raw OpenQuake fractiles.

See ?fitSaF and ?fitSaFMixture 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). The draw table retains the sampled PGA as PGA.s (g) for the same realisation across every active model and ky.

Stage 3 — Coefficient inversion (invertDnDraws)

Each curve is projected monotone and inverted in log-log space to find the yield acceleration km,s such that Dm,s(km,s) = d*. The inversion is performed once for each model m and realisation s.

If a complete AM88 curve lies on its ky / PGA.s = 0.9999 no-sliding limit because PGA.s is below the stored ky grid, the retained PGA.s and plateau displacement define the unique root on 0 < ky / PGA.s < 0.9999. Other flat supports remain explicit errors.

The root distributions use the same epistemic mixture as the displacement draws:

FK(k)=∑mwmFKm(k). F_K(k)=\sum_m w_m F_{K_m}(k).

Every root in model m receives the model weight wm when the weighted mean and quantiles are evaluated. The weights are normalised across the active family. This preserves each model distribution instead of pairing models by the otherwise arbitrary simulation label s.

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().
  • Zero-displacement screening: the flexible-block papers define a discrete probability mass P(D = 0) (negligible displacement) ahead of the lognormal branch. The package intentionally omits it and returns the unconditional lognormal — a conservative convention (small displacements are over- rather than under-predicted).
  • Near-fault pulse behaviour (BM19) is a declared attribute of the record/scenario (pulse argument, default ordinary), not an automatic PGV-threshold switch.

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: invertDnDraws() defaults to the mean plus p ∈ {0.05, 0.10, 0.16, 0.50, 0.84, 0.90, 0.95}; other levels are requested via its p argument.

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.
  • Stewart et al. (2020). Earthquake Spectra 36(1). doi:10.1177/8755293019878185.
  • Hashash et al. (2020). Earthquake Spectra 36(1). doi:10.1177/8755293019878193.
  • Hashash et al. (2017). PEER Report 2017/05.