# -*- coding: utf-8 -*-
# vim: tabstop=4 shiftwidth=4 softtabstop=4
#
# Copyright (C) 2014-2023 GEM Foundation
#
# OpenQuake is free software: you can redistribute it and/or modify it
# under the terms of the GNU Affero General Public License as published
# by the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# OpenQuake is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Affero General Public License for more details.
#
# You should have received a copy of the GNU Affero General Public License
# along with OpenQuake. If not, see <http://www.gnu.org/licenses/>.
"""
Module exports :class:`TromansEtAl2019`
"""
import numpy as np
from openquake.hazardlib.gsim.base import registry, GMPE, CoeffsTable
from openquake.hazardlib import const
from openquake.hazardlib.imt import from_string
#: Coefficient tables as per annex B of Abrahamson et al. (2014)
ASK_TAU_COEFFS = CoeffsTable(sa_damping=5, table="""\
IMT s3 s4
pga 0.47 0.36
pgv 0.38 0.38
0.010 0.47 0.36
0.020 0.47 0.36
0.030 0.47 0.36
0.050 0.47 0.36
0.075 0.47 0.36
0.100 0.47 0.36
0.150 0.47 0.36
0.200 0.47 0.36
0.250 0.47 0.36
0.300 0.47 0.36
0.400 0.47 0.36
0.500 0.47 0.36
0.750 0.47 0.36
1.000 0.47 0.36
1.500 0.47 0.36
2.000 0.47 0.36
3.000 0.47 0.36
4.000 0.47 0.36
5.000 0.47 0.36
6.000 0.47 0.36
7.500 0.47 0.36
10.00 0.47 0.36
""")
#: Single station intra-event term for the constant and magnitude-dependent
#: model, as provided in Table 4 of Rodriguez-Marek et al (2013)
PHI_SS_COEFFS = CoeffsTable(sa_damping=5, table="""\
IMT const_phiss phi1m phi2m
pgv 0.46 0.49 0.37
pga 0.46 0.49 0.35
0.01 0.46 0.49 0.35
0.1 0.45 0.45 0.43
0.2 0.48 0.51 0.37
0.3 0.48 0.51 0.37
0.5 0.46 0.49 0.37
1.0 0.45 0.46 0.40
3.0 0.41 0.41 0.41
""")
#: Single ststion dfS2S model
DELTA_PHI_S2S = CoeffsTable(sa_damping=5, table="""\
IMT dfs2s
pga 0.0000
0.01000 0.0000
0.99999 0.0000
1.00000 0.0001
2.00000 0.0069
3.00000 0.0131
10.0000 0.0131
""")
[docs]def get_heteroskedastic_tau(imt, mag):
"""
Returns the magnitude dependent inter-event variability using the
model of Abrahamson et al (2014).
:param dict C:
Coefficients dictionary
:param float mag:
Magnitude
"""
C = ASK_TAU_COEFFS[imt]
tau = C['s3'] + (C['s4'] - C['s3']) / 2. * (mag - 5.)
tau[mag < 5] = C['s3']
tau[mag > 7] = C['s4']
return tau
[docs]def get_heteroskedastic_phi(imt, mag):
"""
Returns the heteroskedastic intra-event term, taken as the maximum
of the constant single-station phi and the magnitude dependent
single-station phi provided in Table 4 of Rodriguez-Marek et al (2014)
"""
C = PHI_SS_COEFFS[imt]
mag_phi = C["phi1m"] + (C["phi2m"] - C["phi1m"]) * (mag - 5.0) / 2.0
mag_phi[mag < 5.0] = C["phi1m"]
mag_phi[mag > 7.0] = C["phi2m"]
mag_phi[mag_phi <= C["const_phiss"]] = C["const_phiss"]
return mag_phi
HETEROSKEDASTIC_PHI = {
"upper": lambda imt, mag: 1.16 * get_heteroskedastic_phi(imt, mag),
"central": lambda imt, mag: get_heteroskedastic_phi(imt, mag),
"lower": lambda imt, mag: 0.84 * get_heteroskedastic_phi(imt, mag)}
HOMOSKEDASTIC_PHI = {
"upper": lambda imt: 1.16 * PHI_SS_COEFFS[imt]["const_phiss"],
"central": lambda imt: PHI_SS_COEFFS[imt]["const_phiss"],
"lower": lambda imt: 0.84 * PHI_SS_COEFFS[imt]["const_phiss"]}
HETEROSKEDASTIC_TAU = {
"upper": lambda imt, mag: get_heteroskedastic_tau(imt, mag) + 0.075,
"central": lambda imt, mag: get_heteroskedastic_tau(imt, mag),
"lower": lambda imt, mag: get_heteroskedastic_tau(imt, mag) - 0.075}
SIX = np.array([6.0])
HOMOSKEDASTIC_TAU = {
"upper": lambda imt: get_heteroskedastic_tau(imt, SIX) + 0.075,
"central": lambda imt: get_heteroskedastic_tau(imt, SIX),
"lower": lambda imt: get_heteroskedastic_tau(imt, SIX) - 0.075}
uppernames = '''
DEFINED_FOR_INTENSITY_MEASURE_TYPES
DEFINED_FOR_STANDARD_DEVIATION_TYPES
REQUIRES_SITES_PARAMETERS
REQUIRES_RUPTURE_PARAMETERS
REQUIRES_DISTANCES
'''.split()
[docs]def get_alatik_youngs_sigma_mu(mag, rake, imt):
"""
Implements the statistical uncertainty model of Al Atik & Youngs (2014)
given in equations 9 to 11 in the manuscript.
"""
if str(imt) == "PGA":
period = 0.01
elif str(imt).startswith("SA"):
period = imt.period
else:
raise ValueError("Al Atik & Youngs (2014) Model not supported "
"for %s" % str(imt))
sigma_mu = np.where(mag >= 7., 0.056 * (mag - 7.0) + 0.083, 0.083)
if period >= 1.0:
sigma_mu += 0.0171 * np.log(period)
# Normal faulting case
sigma_mu[(rake >= -135.) & (rake <= -45.)] += 0.038
return sigma_mu
[docs]def get_stddevs(branch, phi_ds2s, homoskedastic_sigma, imt, mag):
"""
Returns the standard deviations as described in Figure 10 and
section 4 of Tromans et al. (2019).
"""
if homoskedastic_sigma:
# Homoskedastic sigma branch
tau = HOMOSKEDASTIC_TAU[branch](imt)
phi = HOMOSKEDASTIC_PHI[branch](imt)
else:
# Heteroskedastic sigma branch
tau = HETEROSKEDASTIC_TAU[branch](imt, mag)
phi = HETEROSKEDASTIC_PHI[branch](imt, mag)
# Add on the delta phi_d2s
if phi_ds2s:
phi = np.sqrt(phi ** 2 + DELTA_PHI_S2S[imt]["dfs2s"] ** 2)
return [np.sqrt(tau ** 2 + phi ** 2), tau, phi]
[docs]class TromansEtAl2019(GMPE):
"""
Implements a modifiable GMPE to apply the standard deviation model and
adjustments described in Tromans et al. (2019), for application to
a nuclear power plant site in the UK:
Tromans, I. J., Aldama-Bustos, G., Douglas, J., Lessi-Cheimariou, A.,
Hunt, S., Davi, M., Musson, R. M. W., Garrard, G., Strasser, F. and
Robertson, C. (2019) "Probabilistic seismic hazard assessment for a
new-build nuclear power plant site in the UK", Bulletin of Earthquake
Engineering, 17: 1- 36
:param gmpe:
The GMPE for calculation of the medeian ground motion model
:param string branch:
The model defines three branches for the different aleatory
uncertainty models "lower", "central" and "upper"
:param float scaling_factor:
Factor to scale the median values of the GMPE to account for, for
example, stress drop uncertainty
:param bool homoskedastic sigma:
Determines whether to use the homoskedastic uncertainty model (True)
or the heteroskedastic model (False)
:param vskappa:
Apply vs-kappa adjustment factors defined using a dictionary organised
by IMT, or else none.
:param phi_ds2s:
Adds the phi_ds2s term to the sigma model (True) or retains the
single station model
"""
#: Supported tectonic region type is 'active shallow crust'
DEFINED_FOR_TECTONIC_REGION_TYPE = const.TRT.ACTIVE_SHALLOW_CRUST
#: Set of :mod:`intensity measure types <openquake.hazardlib.imt>`
#: this GSIM can calculate. A set should contain classes from module
#: :mod:`openquake.hazardlib.imt`.
DEFINED_FOR_INTENSITY_MEASURE_TYPES = set()
#: Supported intensity measure component is the geometric mean of two
#: horizontal components
DEFINED_FOR_INTENSITY_MEASURE_COMPONENT = const.IMC.RotD50
#: Supported standard deviation types are inter-event, intra-event
#: and total
DEFINED_FOR_STANDARD_DEVIATION_TYPES = {
const.StdDev.TOTAL, const.StdDev.INTER_EVENT, const.StdDev.INTRA_EVENT}
#: Required site parameter is not set
REQUIRES_SITES_PARAMETERS = set()
#: Required rupture parameters are magnitude, others will be taken from
#: the GMPE
REQUIRES_RUPTURE_PARAMETERS = {'mag'}
#: Required distance measure will be set by the GMPE
REQUIRES_DISTANCES = set()
def __init__(self, gmpe_name, branch="central", sigma_mu_epsilon=0.0,
homoskedastic_sigma=False, scaling_factor=None,
vskappa=None, phi_ds2s=True, **kwargs):
super().__init__(gmpe_name=gmpe_name, branch=branch,
homoskedastic_sigma=homoskedastic_sigma,
scaling_factor=scaling_factor, vskappa=vskappa,
phi_ds2s=phi_ds2s, **kwargs)
self.gmpe = registry[gmpe_name]()
# Update the required_parameters
for name in uppernames:
setattr(self, name,
frozenset(getattr(self, name) | getattr(self.gmpe, name)))
# If the scaling factor take the natural log
if scaling_factor:
self.scaling_factor = np.log(scaling_factor)
else:
self.scaling_factor = None
# If vs-kappa is passed as a dictionary then transform to CoeffsTable
if isinstance(vskappa, dict):
in_vskappa = {}
for key in vskappa:
in_vskappa[from_string(key)] = {"vskappa":
np.log(vskappa[key])}
self.vskappa = CoeffsTable.fromdict(in_vskappa)
else:
self.vskappa = None
self.branch = branch
self.sigma_mu_epsilon = sigma_mu_epsilon
self.homoskedastic_sigma = homoskedastic_sigma
self.phi_ds2s = phi_ds2s
[docs] def compute(self, ctx: np.recarray, imts, mean, sig, tau, phi):
"""
Returns the mean and standard deviations applying, where specified,
scalar adjustment and vs-kappa adjustment to the mean from the
original GMPE.
"""
# Retrieve the mean values from the GMPE
self.gmpe.compute(ctx, imts, mean, sig, tau, phi)
# Apply scaling factor
if self.scaling_factor:
mean += self.scaling_factor
for m, imt in enumerate(imts):
# Apply sigma_mu epsilon factor
if self.sigma_mu_epsilon:
mean += (self.sigma_mu_epsilon *
get_alatik_youngs_sigma_mu(ctx.mag, ctx.rake, imt))
# Apply vs-kappa correction
if self.vskappa:
mean += self.vskappa[imt]["vskappa"]
# Get stddevs
sig[m], tau[m], phi[m] = get_stddevs(
self.branch, self.phi_ds2s, self.homoskedastic_sigma,
imt, ctx.mag)
[docs]class TromansEtAl2019SigmaMu(TromansEtAl2019):
"""
Extension of the Tromans et al. (2019) to facilitate the application
of the statistical uncertainty (sigma_mu) adjustment using the factors
described by Al Atik & Youngs (2014)
Al Atik, L. and Youngs, R. R. (2014) "Epistemic Uncertainty for
NGA-West 2 Models", Earthquake Spectra, 30(3): 1301 - 1318
"""
#: Required rupture parameters are magnitude and style of faulting, others
#: will be taken from the GMPE
REQUIRES_RUPTURE_PARAMETERS = {'mag', 'rake'}
