# Newmark ensemble formulation

## 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:

\[ 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`, `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.

- Mode A (mean-only input): rock Sa is deterministic; dispersion comes from the site-model sigma alone.
- Mode B (quantile input): Sa is sampled per period from the empirical quantile function, with the PGA–Sa(Tn) dependence represented by a per-period two-variable Gaussian copula (Baker & Jayaram 2008 correlation); dispersion reflects combined hazard and site-factor uncertainty. The full inter-period correlation matrix is used downstream in `fitDnCurve()`, not in this per-period stage.
- Mixture input (`fitSaFMixture`): when the rock hazard is an analytic lognormal mixture (DSHA scenario engine), the composition with the site-factor lognormal is evaluated exactly by quadrature — no sampling, byte-stable output.

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).

### 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()`.
- 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⁴\]. 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)\] g (default of `getDnKy()`).
- Reported quantiles: `invertDnDraws()` defaults to the mean plus p ∈ {0.05, 0.10, 0.16, 0.84, 0.90, 0.95}; the median and any 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.