Source code for openquake.hazardlib.gsim.si_2020

# -*- coding: utf-8 -*-
# vim: tabstop=4 shiftwidth=4 softtabstop=4
#
# Copyright (C) 2014-2022 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:`SiEtAl2020SInter`
               :class:`SiEtAl2020SSlab`
"""
import numpy as np
from openquake.hazardlib.gsim.base import GMPE, CoeffsTable
from openquake.hazardlib import const
from openquake.hazardlib.imt import PGA, PGV, SA


[docs]def get_base_term(trt, C): """ Returns the constant term of the GMM. Depends on tectonic type """ return C["e"] +\ (C["d0"] if trt == const.TRT.SUBDUCTION_INTERFACE else C["d1"])
[docs]def get_magnitude_scaling_term(C, imt, mag): """ Returns the magnitude scaling term (Equation 3.16) """ if imt.string in ("PGA", "PGV"): m_ref = 8.3 else: if imt.period < 2.0: m_ref = 8.3 else: m_ref = 7.5 fmag = np.where( mag < m_ref, C["a1"] * mag, C["a1"] * mag + (C["a2"] - C["a1"]) * (mag - m_ref) ) return fmag
[docs]def get_depth_scaling_term(C, hypo_depth): """ Returns the depth scaling term (Eqation 3.16) """ return C["h"] * hypo_depth
[docs]def get_anelastic_attenuation_term(C, rrup): """ Returns the anelastic attenuation term (Eq. 3.15) The period dependent coefficients are calculated and added to the coefficients table. """ return C["c_attn"] * rrup
[docs]def get_moho_depth(rup): """ get Moho depth dependent on hypocenter location for now, return 30km everywhere """ return 30.0
[docs]def get_geometric_attenuation_term(C, rup, rrup): """ Returns the geometric attenuation term (Eq. 3.13/3.14) Period dependent coefficients are calculated and added to the coefficients table. """ mref = np.where(rup.mag < 8.3, rup.mag, 8.3) c = C["c_gs"] * 10.0 ** (0.5 * mref) zmoho = get_moho_depth(rup) idx = (rup.hypo_depth > zmoho) & (rrup >= (1.7 * rup.hypo_depth)) fgs = -1.0 * np.log10(rrup + c) fgs[idx] = 0.6 * np.log10(1.7 * rup.hypo_depth + c) -\ 1.6 * np.log10(rrup[idx] + c) return fgs
[docs]def get_shallow_site_response_term(C, vs30, pga760): """ Returns the shallow site response term (Eq. 3.2 tp 3.4) """ vs30_comp = np.clip(vs30, -np.inf, 760.0) # Nonlinear site scaling term f2 = 0.5 * C["f4"] * (np.exp(C["f5"] * (vs30_comp - 360.0)) - np.exp(C["f5"] * (760.0 - 360.0))) # Linear site term f_site_lin = np.where(vs30 <= C["Vc"], C["c"] * np.log(vs30 / C["Vref"]), C["c"] * np.log(C["Vc"] / C["Vref"]) ) f_site_nl = np.full_like(vs30, C["f1"]) idx = f2 != 0 f_site_nl[idx] = f_site_nl[idx] + f2[idx] * np.log( (pga760[idx] + C["f3"]) / C["f3"] ) return (f_site_lin + f_site_nl) / np.log(10.0)
[docs]def get_basin_response_term(C, z_value): """ Returns the basin response term (Eq. 3.10) """ return C["Cd"] + C["Dd"] * z_value
def _get_pga_rock(C, trt, imt, sites, rup, dists): """ Returns the PGA on rock for site response """ mean = (get_base_term(trt, C) + get_magnitude_scaling_term(C, imt, rup.mag) + get_geometric_attenuation_term(C, rup, dists.rrup) + get_basin_response_term(C, sites.z2pt5) + get_depth_scaling_term(C, rup.hypo_depth) + get_anelastic_attenuation_term(C, dists.rrup)) return 10.0 ** mean
[docs]def get_mean_values(C, trt, imt, sites, rup, dists, a760): """ Returns the mean values for a specific IMT """ mean = np.log(10.0) * (get_base_term(trt, C) + get_magnitude_scaling_term(C, imt, rup.mag) + get_depth_scaling_term(C, rup.hypo_depth) + get_geometric_attenuation_term(C, rup, dists.rrup) + get_anelastic_attenuation_term(C, dists.rrup) + get_basin_response_term(C, sites.z2pt5) + get_shallow_site_response_term(C, sites.vs30, a760)) return mean
[docs]class SiEtAl2020SInter(GMPE): """ Implements NGA Subduction model of Si, Midorikawa, Kishida (2020) for interface events Si H, Midorikawa S, Kishida T (2020) "Development of NGA-Sub Ground-Motion Model of 5%-Damped Psuedo-Spectral Acceleration Based on Database for Subduction Earthquakes in Japan" PEER Report No. 2020/06 Implementation is based on preliminary PEER report and R implementation obtained from T. Kishida on 09/16/2020 """ #: 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, PGV, SA} #: Supported intensity measure component is orientation-independent #: average horizontal :attr:`~openquake.hazardlib.const.IMC.GMRotI50` DEFINED_FOR_INTENSITY_MEASURE_COMPONENT = const.IMC.RotD50 #: Supported standard deviation types are inter-event, intra-event #: and total, see section "Aleatory Variability Model", page 1094. DEFINED_FOR_STANDARD_DEVIATION_TYPES = { const.StdDev.TOTAL, const.StdDev.INTER_EVENT, const.StdDev.INTRA_EVENT} #: Required site parameters are Vs30 and Z2.5 REQUIRES_SITES_PARAMETERS = {'vs30', 'z2pt5'} #: Required rupture parameters are magnitude and hypocentral deph REQUIRES_RUPTURE_PARAMETERS = {'mag', 'hypo_depth'} #: Required distance measure is Rrup REQUIRES_DISTANCES = {'rrup'}
[docs] def compute(self, ctx, imts, mean, sig, tau, phi): """ See :meth:`superclass method <.base.GroundShakingIntensityModel.compute>` for spec of input and result values. """ trt = self.DEFINED_FOR_TECTONIC_REGION_TYPE # extract dictionaries of coefficients specific to PGA # intensity measure type and for PGA C_PGA = self.COEFFS[PGA()] # Get mean PGA on rock (Vs30 760 m/s) pga760 = _get_pga_rock(C_PGA, trt, PGA(), ctx, ctx, ctx) for m, imt in enumerate(imts): # Get the coefficients for the IMT C = self.COEFFS[imt] # Get mean and standard deviations for IMT mean[m] = get_mean_values(C, trt, imt, ctx, ctx, ctx, pga760) tau[m] = C["tau"] phi[m] = C["phi"] sig[:] = np.sqrt(tau ** 2.0 + phi ** 2.0)
COEFFS = CoeffsTable(sa_damping=5, table="""\ imt e a1 d0 d1 a2 h Cd Dd c Vc Vref f1 f3 f4 f5 phi tau sigma Mb c_gs c_attn pgv -1.895274 0.638198 0.0 0.194161 0.000000 0.005350 0.000 0.000 -0.84000 1300.00 760.0 0.0 0.1 -0.100000 -0.00844 0.557 0.445 0.713 8.3 0.00550 -0.003000 pga -2.685641 0.485385 0.0 0.224236 0.000000 0.006963 0.000 0.000 -0.60000 1500.00 760.0 0.0 0.1 -0.150000 -0.00701 0.646 0.593 0.877 8.3 0.00550 -0.003000 0.010 -2.646637 0.493277 0.0 0.225570 0.000000 0.007054 0.000 0.000 -0.60372 1500.20 760.0 0.0 0.1 -0.148330 -0.00701 0.632 0.565 0.848 8.3 0.00550 -0.003000 0.020 -2.646100 0.494915 0.0 0.232876 0.000000 0.007012 0.000 0.000 -0.57388 1500.40 760.0 0.0 0.1 -0.147100 -0.00728 0.630 0.564 0.845 8.3 0.00550 -0.003000 0.030 -2.629178 0.496553 0.0 0.240181 0.000000 0.007138 0.000 0.000 -0.53414 1503.00 760.0 0.0 0.1 -0.154850 -0.00735 0.633 0.578 0.857 8.3 0.00550 -0.003000 0.050 -2.567020 0.499829 0.0 0.250252 0.000000 0.007786 0.000 0.000 -0.45795 1501.40 760.0 0.0 0.1 -0.196330 -0.00647 0.664 0.634 0.918 8.3 0.00550 -0.003000 0.075 -2.484857 0.503105 0.0 0.258333 0.000000 0.008437 0.000 0.000 -0.44411 1494.00 760.0 0.0 0.1 -0.228660 -0.00573 0.714 0.688 0.992 8.3 0.00550 -0.003000 0.100 -2.456607 0.508019 0.0 0.265100 0.000000 0.008569 0.000 0.000 -0.48724 1479.10 760.0 0.0 0.1 -0.249160 -0.00560 0.743 0.666 0.998 8.3 0.00550 -0.003000 0.150 -2.485473 0.516208 0.0 0.266446 0.000000 0.008411 0.000 0.000 -0.57962 1442.90 760.0 0.0 0.1 -0.257130 -0.00585 0.748 0.563 0.936 8.3 0.00550 -0.003000 0.200 -2.545081 0.524393 0.0 0.261442 0.000000 0.007341 0.000 0.000 -0.68762 1392.60 760.0 0.0 0.1 -0.246580 -0.00614 0.736 0.532 0.908 8.3 0.00550 -0.003000 0.250 -2.636311 0.532568 0.0 0.255371 0.000000 0.006546 0.000 0.000 -0.77177 1356.20 760.0 0.0 0.1 -0.235740 -0.00644 0.713 0.496 0.868 8.3 0.00550 -0.003000 0.300 -2.757882 0.540729 0.0 0.245784 0.000000 0.006321 0.000 0.000 -0.84165 1308.50 760.0 0.0 0.1 -0.219120 -0.00670 0.694 0.469 0.838 8.3 0.00550 -0.002996 0.400 -3.026741 0.556994 0.0 0.232084 0.000000 0.005857 0.000 0.000 -0.91092 1252.70 760.0 0.0 0.1 -0.195820 -0.00713 0.669 0.440 0.801 8.3 0.00438 -0.002581 0.500 -3.273455 0.573162 0.0 0.219238 0.000000 0.005350 0.000 0.000 -0.96930 1203.90 760.0 0.0 0.1 -0.170410 -0.00744 0.654 0.417 0.776 8.3 0.00351 -0.002259 0.750 -3.784828 0.613050 0.0 0.204750 0.000000 0.004290 0.000 0.000 -1.01540 1147.60 760.0 0.0 0.1 -0.138660 -0.00812 0.659 0.419 0.781 8.3 0.00280 -0.002000 1.000 -4.168096 0.651986 0.0 0.195413 0.044000 0.003904 0.008 0.056 -1.05000 1109.90 760.0 0.0 0.1 -0.105210 -0.00844 0.670 0.425 0.794 8.3 0.00280 -0.002000 1.500 -4.952279 0.725900 0.0 0.188847 0.134000 0.003848 0.030 0.067 -1.04540 1072.40 760.0 0.0 0.1 -0.067941 -0.00771 0.701 0.448 0.832 8.3 0.00280 -0.002000 2.000 -5.573430 0.792692 0.0 0.188157 0.397000 0.004260 0.037 0.081 -1.03920 1009.50 760.0 0.0 0.1 -0.036136 -0.00479 0.714 0.469 0.854 7.5 0.00280 -0.002000 3.000 -6.603707 0.898490 0.0 0.180567 0.442428 0.005716 0.022 0.108 -1.01120 922.43 760.0 0.0 0.1 -0.013577 -0.00183 0.699 0.498 0.858 7.5 0.00280 -0.002000 4.000 -7.279363 0.967563 0.0 0.171225 0.484387 0.005425 -0.021 0.142 -0.96938 844.48 760.0 0.0 0.1 -0.003212 -0.00152 0.669 0.547 0.864 7.5 0.00280 -0.002000 5.000 -7.708043 1.008058 0.0 0.163298 0.521126 0.005093 -0.072 0.181 -0.91954 793.13 760.0 0.0 0.1 -0.000255 -0.00144 0.641 0.564 0.854 7.5 0.00280 -0.002000 7.500 -8.163116 1.040264 0.0 0.151111 0.574112 0.004597 -0.114 0.198 -0.77665 771.01 760.0 0.0 0.1 -0.000055 -0.00137 0.638 0.581 0.863 7.5 0.00280 -0.002000 10.00 -8.366713 1.058093 0.0 0.132831 0.554450 0.001817 -0.133 0.190 -0.65575 775.00 760.0 0.0 0.1 0.000000 -0.00136 0.615 0.614 0.869 7.5 0.00280 -0.002000 """)
[docs]class SiEtAl2020SSlab(SiEtAl2020SInter): """ Implements NGA Subduction model of Si, Midorikawa, Kishida (2020) For Intraslab events. """ #: Supported tectonic region type is subduction interface DEFINED_FOR_TECTONIC_REGION_TYPE = const.TRT.SUBDUCTION_INTRASLAB