# -*- coding: utf-8 -*-
# vim: tabstop=4 shiftwidth=4 softtabstop=4
#
# Copyright (C) 2015-2018 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:`AbrahamsonEtAl2018SInter`
:class:`AbrahamsonEtAl2018SInterHigh`
:class:`AbrahamsonEtAl2018SInterLow`
:class:`AbrahamsonEtAl2018SSlab`
:class:`AbrahamsonEtAl2018SSlabHigh`
:class:`AbrahamsonEtAl2018SSlabLow`
"""
import numpy as np
from openquake.hazardlib.gsim.base import GMPE, CoeffsTable
from openquake.hazardlib import const
from openquake.hazardlib.imt import PGA, SA
CONSTANTS = {"n": 1.18,
"c": 1.88,
"C4": 10.0,
"a3": 0.10,
"a5": 0.0,
"a9": 0.4,
"a10": 1.73,
"C1slab": 7.2,
"phiamp": 0.3}
def _compute_pga_rock(slab, C_PGA, ctx):
"""
Returns the PGA on rock (vs30 = 1000 m / s)
"""
lpga1000 = (compute_base_term(slab, C_PGA) +
compute_magnitude_term(slab, C_PGA, ctx.mag) +
compute_depth_term(slab, C_PGA, ctx) +
compute_distance_term(slab, C_PGA, ctx.rrup, ctx.mag))
# Get linear site term for the case where vs30 = 1000.0
flin = _get_linear_site_term(
C_PGA, 1000.0 * np.ones_like(ctx.rrup))
return lpga1000 + flin
[docs]def compute_base_term(slab, C):
"""
Returns the base coefficient of the GMPE, which for interface events
is just the coefficient a1 (adjusted regionally)
"""
if slab:
return C["a1"] + C["a4"] * (CONSTANTS["C1slab"] - C["C1inter"]) + \
CONSTANTS["a10"]
return C["a1"]
[docs]def compute_magnitude_term(slab, C, mag):
"""
Returns the magnitude scaling term
"""
if slab:
# Returns the magnitude scaling term, this time using the constant
# "C1slab" as the hinge magnitude
f_mag = C["a13"] * ((10.0 - mag) ** 2.)
return np.where(mag <= CONSTANTS["C1slab"],
C["a4"] * (mag - CONSTANTS["C1slab"]) + f_mag,
# parameter "a5" is zero, so linear term disappears
f_mag)
f_mag = C["a13"] * (10.0 - mag) ** 2.
return np.where(mag <= C["C1inter"],
C["a4"] * (mag - C["C1inter"]) + f_mag,
# C["a5"] is zero so linear term disappears
f_mag)
[docs]def compute_distance_term(slab, C, rrup, mag):
"""
Returns the distance attenuation
"""
scale = _get_magnitude_scale(slab, C, mag)
fdist = scale * np.log(rrup + CONSTANTS["C4"] *
np.exp(CONSTANTS["a9"] * (mag - 6.0)))
return fdist + C["a6"] * rrup
def _get_magnitude_scale(slab, C, mag):
"""
Returns the magnitude scaling term that modifies the distance
attenuation
"""
if slab:
return C["a2"] + C["a14"] + CONSTANTS["a3"] * (mag - 7.8)
return C["a2"] + CONSTANTS["a3"] * (mag - 7.8)
[docs]def compute_depth_term(slab, C, ctx):
"""
No top of rupture depth term for interface events
"""
if slab: # Equation on P11
return np.where(ctx.ztor <= 100.0,
C["a11"] * (ctx.ztor - 60.0),
C["a11"] * (100.0 - 60.0))
return 0.0
[docs]def compute_site_term(C, vs30, pga1000):
"""
Returns the site amplification
"""
vsstar = np.copy(vs30)
vsstar[vs30 >= 1000.0] = 1000.0
f_site = np.zeros_like(vs30)
# Consider the cases of only linear amplification
idx = vs30 >= C["vlin"]
if np.any(idx):
f_site[idx] = _get_linear_site_term(C, vsstar[idx])
# Consider now only the nonlinear amplification cases
idx = np.logical_not(idx)
if np.any(idx):
# Linear term
flin = C["a12"] * np.log(vsstar[idx] / C["vlin"])
# Nonlinear term
fnl = (-C["b"] * np.log(pga1000[idx] + CONSTANTS["c"])) +\
(C["b"] * np.log(pga1000[idx] + CONSTANTS["c"] *
((vsstar[idx] / C["vlin"]) ** CONSTANTS["n"])))
f_site[idx] = flin + fnl
return f_site
def _get_linear_site_term(C, vsstar):
"""
As the linear site scaling is used for both the pga1000 case and the
general case the relevant common part is returned here
"""
return (C["a12"] + C["b"] * CONSTANTS["n"]) *\
np.log(vsstar / C["vlin"])
[docs]def get_stddevs(C, C_PGA, pga1000, vs30):
"""
Returns the standard deviations
"""
dln = _get_dln_amp(C, pga1000, vs30)
tau = get_inter_event_stddev(C, C_PGA, dln)
phi = get_within_event_stddev(C, C_PGA, dln)
return [np.sqrt(tau ** 2 + phi ** 2), tau, phi]
def _get_dln_amp(C, pga1000, vs30):
"""
Returns the partial deriviative of the amplification term with respect
to pga1000
"""
dln = np.zeros(vs30.shape)
idx = vs30 < C["vlin"]
if np.any(idx):
dln[idx] = C["b"] * pga1000[idx] * (
(-1. / (pga1000[idx] + CONSTANTS["c"])) +
(1. / (pga1000[idx] + CONSTANTS["c"] *
((vs30[idx] / C["vlin"]) ** CONSTANTS["n"]))))
return dln
[docs]def get_within_event_stddev(C, C_PGA, dln):
"""
Returns the within-event aleatory uncertainty, phi
"""
phi_amp2 = CONSTANTS["phiamp"] ** 2.
phi_b = np.sqrt(C["phi0"] ** 2. - phi_amp2)
phi_b_pga = np.sqrt((C_PGA["phi0"] ** 2.) - phi_amp2)
phi = (C["phi0"] ** 2.) +\
((dln ** 2.) * (phi_b ** 2.)) +\
(2.0 * dln * phi_b * phi_b_pga * C["rho_w"])
return np.sqrt(phi)
[docs]def get_inter_event_stddev(C, C_PGA, dln):
"""
Returns the between event aleatory uncertainty, tau
"""
tau = (C["tau0"] ** 2.) +\
((dln ** 2.) * (C["tau0"] ** 2.)) +\
(2.0 * dln * C["tau0"] * C_PGA["tau0"] * C["rho_b"])
return np.sqrt(tau)
[docs]class AbrahamsonEtAl2018SInter(GMPE):
"""
Implements the 2018 updated Abrahamson et al. (2018) "BC Hydro" GMPE for
application to subduction earthquakes, for the case of subduction interface
events.
Abrahamson, N. A., Keuhn, N., Gulerce, Z., Gregor, N., Bozognia, Y.,
Parker, G., Stewart, J., Chiou, B., Idriss, I. M., Campbell, K. and
Youngs, R. (2018) "Update of the BC Hydro Subduction Ground-Motion Model
using the NGA-Subduction Dataset", Pacific Earthquake Engineering
Research Center (PEER) Technical Report, PEER 2018/02
Whilst the original model provides coefficients for different regional
variations, these are incomplete for the purpose of a working
implementation. As the authors indicate in the source, the coefficients
and adjustment factors are intended only for application to the
Cascadia region; hence this is the only version implemented here.
Furthermore, scalar adjustments are intended to be applied in order to
define an "upper", "central" and "lower" branch to cover the epistemic
uncertainty of the core model.
"""
#: Supported tectonic region type is subduction interface
DEFINED_FOR_TECTONIC_REGION_TYPE = const.TRT.SUBDUCTION_INTERFACE
#: Supported intensity measure types are spectral acceleration,
#: and peak ground acceleration
DEFINED_FOR_INTENSITY_MEASURE_TYPES = {PGA, SA}
#: Supported intensity measure component is the geometric mean component
DEFINED_FOR_INTENSITY_MEASURE_COMPONENT = const.IMC.GEOMETRIC_MEAN
#: Supported standard deviation types are inter-event, intra-event
#: and total, see section 4.5
DEFINED_FOR_STANDARD_DEVIATION_TYPES = {
const.StdDev.TOTAL, const.StdDev.INTER_EVENT, const.StdDev.INTRA_EVENT}
#: Site amplification is dependent only upon Vs30
REQUIRES_SITES_PARAMETERS = {'vs30'}
#: Required rupture parameters are only magnitude for the interface model
REQUIRES_RUPTURE_PARAMETERS = {'mag'}
#: Required distance measure is closest distance to rupture, for
#: interface events
REQUIRES_DISTANCES = {'rrup'}
#: A "low" and "high" epistemic adjustment factor will be applied to
#: subclasses of this model
EPISTEMIC_ADJUSTMENT = None
#: Adjustment variable to match Cascadia to global average
CASCADIA_ADJUSTMENT = "adj_int"
[docs] def compute(self, ctx: np.recarray, imts, mean, sig, tau, phi):
"""
See :meth:`superclass method
<.base.GroundShakingIntensityModel.compute>`
for spec of input and result values.
"""
C_PGA = self.COEFFS[PGA()]
slab = self.CASCADIA_ADJUSTMENT == "adj_slab"
for m, imt in enumerate(imts):
C = self.COEFFS[imt]
# compute median pga on rock (vs30=1000), needed for site response
# term calculation
pga1000 = np.exp(_compute_pga_rock(slab, C_PGA, ctx) +
C_PGA[self.CASCADIA_ADJUSTMENT])
# Get full model
mean[m] = (compute_base_term(slab, C) +
compute_magnitude_term(slab, C, ctx.mag) +
compute_depth_term(slab, C, ctx) +
compute_distance_term(slab, C, ctx.rrup, ctx.mag) +
compute_site_term(C, ctx.vs30, pga1000))
sig[m], tau[m], phi[m] = get_stddevs(C, C_PGA, pga1000, ctx.vs30)
if self.EPISTEMIC_ADJUSTMENT:
mean[m] += (C[self.CASCADIA_ADJUSTMENT] +
C[self.EPISTEMIC_ADJUSTMENT])
else:
mean[m] += C[self.CASCADIA_ADJUSTMENT]
COEFFS = CoeffsTable(sa_damping=5, table="""\
imt C1inter vlin b a1 a2 a4 a6 a11 a12 a13 a14 adj_int adj_slab phi0 tau0 rho_w rho_b SINTER_LOW SINTER_HIGH SSLAB_LOW SSLAB_HIGH
pga 8.20 865.1 -1.186 2.34047 -1.044 0.59 -0.00705 0.0170 0.818 -0.0135 -0.223 1.04398007004686 0.83429221602170 0.61 0.580 1.0000 1.0000 -0.3 0.3 -0.50 0.50
0.010 8.20 865.1 -1.186 2.34047 -1.044 0.59 -0.00705 0.0170 0.818 -0.0135 -0.223 1.04398007004686 0.83429221602170 0.61 0.580 1.0000 1.0000 -0.3 0.3 -0.50 0.50
0.020 8.20 865.1 -1.219 2.36005 -1.044 0.59 -0.00707 0.0170 0.857 -0.0135 -0.196 1.04603523427650 0.79173401239125 0.61 0.580 1.0000 1.0000 -0.3 0.3 -0.50 0.50
0.030 8.20 907.8 -1.273 2.38396 -1.080 0.59 -0.00710 0.0170 0.921 -0.0135 -0.128 1.23377425709977 0.70576030699887 0.61 0.580 0.9910 0.9910 -0.3 0.3 -0.50 0.50
0.050 8.20 1053.5 -1.346 2.44598 -1.110 0.59 -0.00725 0.0180 1.007 -0.0138 -0.130 1.34342293514146 0.97601032906970 0.61 0.580 0.9730 0.9730 -0.3 0.3 -0.50 0.50
0.075 8.20 1085.7 -1.471 2.75111 -1.110 0.59 -0.00758 0.0180 1.225 -0.0142 -0.130 1.32123133203975 0.99202176921030 0.61 0.580 0.9520 0.9520 -0.3 0.3 -0.50 0.50
0.100 8.20 1032.5 -1.624 3.01943 -1.110 0.59 -0.00788 0.0180 1.457 -0.0145 -0.130 1.32307413547709 0.99825322734373 0.61 0.580 0.9290 0.9290 -0.3 0.3 -0.50 0.50
0.150 8.20 877.6 -1.931 3.34867 -1.084 0.59 -0.00820 0.0175 1.849 -0.0153 -0.156 1.20604738361364 0.91891214881125 0.61 0.560 0.8960 0.8960 -0.3 0.3 -0.50 0.50
0.200 8.20 748.2 -2.188 3.28397 -1.027 0.62 -0.00835 0.0170 2.082 -0.0162 -0.172 1.14232895808973 0.87856284947151 0.61 0.540 0.8740 0.8740 -0.3 0.3 -0.50 0.50
0.250 8.20 654.3 -2.381 3.21131 -0.983 0.64 -0.00835 0.0160 2.240 -0.0172 -0.184 1.04808298000789 0.81234114576059 0.61 0.520 0.8560 0.8560 -0.3 0.3 -0.46 0.46
0.300 8.20 587.1 -2.518 3.14548 -0.947 0.66 -0.00828 0.0152 2.341 -0.0183 -0.194 0.94565364697989 0.74536286357934 0.61 0.505 0.8410 0.8410 -0.3 0.3 -0.42 0.42
0.400 8.20 503.0 -2.657 2.99656 -0.890 0.68 -0.00797 0.0140 2.415 -0.0206 -0.210 0.79394049645376 0.62046390803996 0.61 0.480 0.8180 0.8180 -0.3 0.3 -0.38 0.38
0.500 8.20 456.6 -2.669 2.83921 -0.845 0.68 -0.00770 0.0130 2.359 -0.0231 -0.223 0.66202915411124 0.54491903264494 0.61 0.460 0.7830 0.7830 -0.3 0.3 -0.34 0.34
0.600 8.20 430.3 -2.599 2.65827 -0.809 0.68 -0.00740 0.0122 2.227 -0.0256 -0.233 0.53601960759514 0.45568525725929 0.61 0.450 0.7315 0.7315 -0.3 0.3 -0.30 0.30
0.750 8.15 410.5 -2.401 2.34577 -0.760 0.68 -0.00698 0.0113 1.949 -0.0296 -0.245 0.35555118688641 0.32581112214331 0.61 0.450 0.6800 0.6800 -0.3 0.3 -0.30 0.30
1.000 8.10 400.0 -1.955 1.85131 -0.698 0.68 -0.00645 0.0100 1.402 -0.0363 -0.261 0.24282500000000 0.22809166666667 0.61 0.450 0.6070 0.6070 -0.3 0.3 -0.30 0.30
1.500 8.05 400.0 -1.025 1.21559 -0.612 0.68 -0.00570 0.0082 0.329 -0.0493 -0.285 -0.07680250000000 -0.01107333333333 0.61 0.450 0.5004 0.5040 -0.3 0.3 -0.30 0.30
2.000 8.00 400.0 -0.299 0.64875 -0.550 0.68 -0.00510 0.0070 -0.487 -0.0610 -0.301 -0.20847000000000 -0.08503000000000 0.61 0.450 0.4301 0.4310 -0.3 0.3 -0.30 0.30
2.500 7.95 400.0 0.000 0.08221 -0.501 0.68 -0.00465 0.0060 -0.770 -0.0711 -0.313 -0.20955750000000 -0.05003166666667 0.61 0.450 0.3795 0.3795 -0.3 0.3 -0.30 0.30
3.000 7.90 400.0 0.000 -0.36926 -0.460 0.68 -0.00430 0.0052 -0.700 -0.0798 -0.323 -0.21660250000000 -0.01508333333333 0.61 0.450 0.3280 0.3280 -0.3 0.3 -0.30 0.30
4.000 7.85 400.0 0.000 -1.03439 -0.455 0.68 -0.00390 0.0040 -0.607 -0.0935 -0.282 -0.06406000000000 0.06239166666667 0.61 0.450 0.2505 0.2550 -0.3 0.3 -0.30 0.30
5.000 7.80 400.0 0.000 -1.51967 -0.450 0.73 -0.00370 0.0030 -0.540 -0.0980 -0.250 0.06459750000000 0.06338666666667 0.61 0.450 0.2000 0.2000 -0.3 0.3 -0.30 0.30
6.000 7.80 400.0 0.000 -1.81025 -0.450 0.78 -0.00357 0.0022 -0.479 -0.0980 -0.250 0.09014000000000 0.09661000000000 0.61 0.450 0.2000 0.2000 -0.3 0.3 -0.30 0.30
7.500 7.80 400.0 0.000 -2.17269 -0.450 0.84 -0.00340 0.0013 -0.393 -0.0980 -0.250 0.13793250000000 0.12577500000000 0.61 0.450 0.2000 0.2000 -0.3 0.3 -0.30 0.30
10.00 7.80 400.0 0.000 -2.71182 -0.450 0.93 -0.00327 0.0000 -0.350 -0.0980 -0.250 0.31434000000000 0.23962833333333 0.61 0.450 0.2000 0.2000 -0.3 0.3 -0.30 0.30
""")
[docs]class AbrahamsonEtAl2018SInterHigh(AbrahamsonEtAl2018SInter):
"""
Abrahamson et al (2018) subduction interface GMPE with the positive
epistemic adjustment factor applied
"""
EPISTEMIC_ADJUSTMENT = "SINTER_HIGH"
[docs]class AbrahamsonEtAl2018SInterLow(AbrahamsonEtAl2018SInter):
"""
Abrahamson et al (2018) subduction interface GMPE with the negative
epistemic adjustment factor applied
"""
EPISTEMIC_ADJUSTMENT = "SINTER_LOW"
[docs]class AbrahamsonEtAl2018SSlab(AbrahamsonEtAl2018SInter):
"""
Abrahamson et al. (2018) updated "BC Hydro" subduction GMPE for application
to subduction in-slab earthquakes.
"""
#: Supported tectonic region type is subduction in-slab
DEFINED_FOR_TECTONIC_REGION_TYPE = const.TRT.SUBDUCTION_INTRASLAB
#: Required rupture parameters for the in-slab model are magnitude and top
# of rupture depth
REQUIRES_RUPTURE_PARAMETERS = {'mag', 'ztor'}
#: Cascadia adjustment factor
CASCADIA_ADJUSTMENT = "adj_slab"
[docs]class AbrahamsonEtAl2018SSlabHigh(AbrahamsonEtAl2018SSlab):
"""
Abrahamson et al (2018) subduction in-slab GMPE with the positive
epistemic adjustment factor applied
"""
EPISTEMIC_ADJUSTMENT = "SSLAB_HIGH"
[docs]class AbrahamsonEtAl2018SSlabLow(AbrahamsonEtAl2018SSlab):
"""
Abrahamson et al (2018) subduction in-slab GMPE with the negative
epistemic adjustment factor applied
"""
EPISTEMIC_ADJUSTMENT = "SSLAB_LOW"