Source code for openquake.hazardlib.gsim.yenier_atkinson_2015

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


[docs]class YenierAtkinson2015BSSA(GMPE): """ Implements the GMM of Yenier and Atkinson (2015) as described in the paper titled "Regionally Adjustable Generic Ground-Motion Prediction Equation to Central and Eastern North America" published on BSSA, vol 105. Note that this model does not provide a standard deviation, hence, in order to use it for PSHA calculations if must be combined with a model for ground motion aleatory uncertainty such as, for example, the one proposed by Al Atik (2014). :param focal_depth: A float defining focal depth [km]. :param region: A string specifying a region. Admitted values are 'CENA' (Central and East North America) and 'CA' (California). Default is 'CENA' """ #: Supported tectonic region type is active shallow crust, see title! DEFINED_FOR_TECTONIC_REGION_TYPE = const.TRT.ACTIVE_SHALLOW_CRUST #: Supported intensity measure types are spectral acceleration, peak #: ground velocity and peak ground acceleration, see tables 4 #: pages 1036 DEFINED_FOR_INTENSITY_MEASURE_TYPES = set([PGA, PGV, SA]) #: Supported intensity measure component is orientation-independent #: average horizontal :attr:`~openquake.hazardlib.const.IMC.RotD50`, #: see page 1025. DEFINED_FOR_INTENSITY_MEASURE_COMPONENT = const.IMC.AVERAGE_HORIZONTAL #: Supported standard deviation types are inter-event, intra-event #: and total, see paragraph "Equations for standard deviations", page #: 1046. DEFINED_FOR_STANDARD_DEVIATION_TYPES = set([const.StdDev.TOTAL]) #: Required site parameter is Vs30 REQUIRES_SITES_PARAMETERS = {'vs30'} #: Required rupture parameters are magnitude and hypocenter depth REQUIRES_RUPTURE_PARAMETERS = {'mag', 'hypo_depth'} #: Required distance measures is Rrup REQUIRES_DISTANCES = {'rrup'} def __init__(self, focal_depth=None, region='CENA', **kwargs): super().__init__(focal_depth=focal_depth, region=region, **kwargs) self.focal_depth = focal_depth self.region = region
[docs] def get_mean_and_stddevs(self, sctx, rctx, dctx, imt, stddev_types): # Compute focal depth if not set at the initialization level if self.focal_depth is None: self.focal_depth = rctx.hypo_depth mean = self._get_mean_on_soil(sctx, rctx, dctx, imt, stddev_types) stddevs = np.zeros_like(sctx.vs30) return mean, stddevs
def _get_mean_on_soil(self, sctx, rctx, dctx, imt, stddev_types): # Get PGA on rock tmp = PGA() pga_rock = self._get_mean_on_rock(sctx, rctx, dctx, tmp, stddev_types) pga_rock = np.exp(pga_rock) # Site-effect model f_s = get_fs_SeyhanStewart2014(imt, pga_rock, sctx.vs30) # Compute the mean on soil mean = self._get_mean_on_rock(sctx, rctx, dctx, imt, stddev_types) mean += f_s return mean def _get_mean_on_rock(self, sctx, rctx, dctx, imt, stddev_types): # Get coefficients C2 = self.TAB2[imt] C3 = self.TAB3[imt] C4 = self.TAB4[imt] # Magnitude effect f_m = self._get_f_m(C2, imt, rctx.mag) # Stress adjustment f_delta_sigma = self._get_stress_drop_adjstment(C3, imt, rctx.mag) # Geometrical spreading f_z = self._get_f_z(C2, imt, dctx.rrup, rctx.mag) # Anelastic attenuation function f_gamma = self._get_f_gamma(C4, imt, dctx.rrup) # Regional term for stress drop c_e = self._get_c_e(imt) # Regional term for path duration c_p = self._get_c_p(imt, dctx.rrup, rctx.mag) # Compute mean using equation 26 mean = f_m + f_delta_sigma + f_z + f_gamma + c_e + c_p return mean def _get_f_m(self, C, imt, m): """ Implements eq. 3 at page 1991 """ mh = C['Mh'] if m <= mh: return C['e0'] + C['e1']*(m-mh) + C['e2']*(m-mh)**2 else: return C['e0'] + C['e3']*(m-mh) def _get_edelta(self, C, m, stress_drop): if stress_drop <= 100: edelta = (C['s0'] + C['s1']*m + C['s2']*m**2 + C['s3']*m**3 + C['s4']*m**4) else: edelta = (C['s5'] + C['s6']*m + C['s7']*m**2 + C['s8']*m**3 + C['s9']*m**4) return edelta def _get_stress_drop_adjstment(self, C, imt, m): """ Implements eq. 4 at page 1991 and eq. 17 at page 1994. For CENA we use eq. 21 at page 2001 """ region = self.region if region == 'CENA': d = self.focal_depth t1 = min([0, 0.290*(d - 10.)]) t2 = min([0, 0.229*(m - 5.)]) delta_sigma = np.exp(5.704 + t1 + t2) edelta = self._get_edelta(C, m, delta_sigma) return edelta * np.log(delta_sigma/100.) else: fmt = '{:s} is a region not supported for stress drop adjustment' msg = fmt.format(region) raise ValueError(msg) def _get_f_z(self, C, imt, rrup, m): """ Implements eq. 7 and eq. 8 at page 1991 """ # Pseudo depth - see eq. 6 at page 1991 pseudo_depth = 10**(-0.405+0.235*m) # Effective distance - see eq. 5 at page 1991 reff = (rrup**2+pseudo_depth**2)**0.5 # The transition_distance is 50 km as defined just below eq. 8 transition_dst = 50. # Geometrical spreading rates b1 = -1.3 b2 = -0.5 # Geometrical attenuation z = reff**b1 ratio_a = reff / transition_dst z[reff > transition_dst] = (transition_dst**b1 * (ratio_a[reff > transition_dst])**b2) # Compute geometrical spreading function ratio_b = reff / (1.+pseudo_depth**2)**0.5 return np.log(z) + (C['b3'] + C['b4']*m)*np.log(ratio_b) def _get_f_gamma(self, C, imt, rrup): """ Implements """ region = self.region if region == 'CENA': return rrup * C['gCENA'] elif region == 'CA': return rrup * C['gCalifornia'] else: fmt = '{:s} is a key not supported for region definition' msg = fmt.format(region) raise ValueError(msg) def _get_c_e(self, imt): """ Implements the Ce calibration term. - For 'CENA' See eq. 23 at page 2003 """ region = self.region if region == 'CENA': # See equation 23 page 2003 of Yenier and Atkinson if str(imt) == 'PGA': return -0.25 elif str(imt) == 'PGV': return -0.21 elif imt.period <= 10.: return -0.25 + np.max([0, 0.39*np.log(imt.period/2)]) else: fmt = 'This IMT is not supported by the Ce calibration term' msg = fmt.format(region) raise ValueError(msg) else: fmt = '{:s} region does not have Ce calibration term' msg = fmt.format(region) raise ValueError(msg) def _get_c_p(self, imt, rrup, m): """ Implements the Cp calibration term """ region = self.region if region == 'CENA': # See equations 24 and 25 page 2003 of Yenier and Atkinson if str(imt) == 'PGA': delta_b3 = 0.030 elif str(imt) == 'PGV': delta_b3 = 0.052 elif imt.period <= 10.: tmp = 0.095*np.log(imt.period/0.065) delta_b3 = np.min([0.095, 0.030+np.max([0, tmp])]) else: msg = 'This region is not supported by the Ce calibration term' raise ValueError(msg) # Compute the calibration term pseudo_depth = 10**(-0.405+0.235*m) reff = (rrup**2+pseudo_depth**2)**0.5 cp = np.zeros_like(rrup) # cp[rrup <= 150] = delta_b3 * np.log(rrup[rrup <= 150]/150.) cp[reff <= 150] = delta_b3 * np.log(reff[reff <= 150]/150.) return cp else: fmt = '{:s} region does not have Cp calibration term' msg = fmt.format(region) raise ValueError(msg) TAB2 = CoeffsTable(sa_damping=5, table="""\ imt Mh e0 e1 e2 e3 b3 b4 0.01 5.85 2.227000 0.6874 -0.13630 0.7643 -0.6209 0.060570 0.013 5.90 2.281000 0.6855 -0.12900 0.7617 -0.6259 0.061290 0.016 5.85 2.272000 0.6971 -0.12320 0.7594 -0.6308 0.061910 0.02 5.90 2.378000 0.6999 -0.10660 0.7488 -0.6377 0.062510 0.025 6.00 2.564000 0.6840 -0.09416 0.7413 -0.6311 0.060970 0.03 6.15 2.806000 0.6607 -0.09087 0.7389 -0.6028 0.056410 0.04 5.75 2.731000 0.7034 -0.10860 0.7383 -0.5484 0.048200 0.05 5.35 2.559000 0.7193 -0.16360 0.7545 -0.5096 0.042790 0.065 5.75 2.997000 0.6842 -0.15470 0.7553 -0.4665 0.036400 0.08 5.20 2.576000 0.7651 -0.24340 0.7865 -0.4210 0.030710 0.1 5.45 2.777000 0.7118 -0.26190 0.7941 -0.3774 0.024720 0.13 5.35 2.641000 0.7346 -0.33210 0.8116 -0.3551 0.022240 0.16 5.25 2.466000 0.8088 -0.38710 0.8407 -0.3265 0.019180 0.2 5.45 2.549000 0.8194 -0.38600 0.8426 -0.2868 0.013760 0.25 5.60 2.517000 0.8671 -0.37750 0.8785 -0.2429 0.009209 0.3 5.85 2.635000 0.8471 -0.36310 0.8763 -0.2117 0.005164 0.4 6.15 2.674000 0.8501 -0.34690 0.8966 -0.1927 0.004847 0.5 6.25 2.544000 0.8856 -0.34860 0.9182 -0.2079 0.008540 0.65 6.60 2.617000 0.8758 -0.31600 0.9251 -0.2277 0.013710 0.8 6.85 2.664000 0.9053 -0.28880 0.8944 -0.2523 0.019060 1 6.45 1.986000 1.3400 -0.24560 0.9829 -0.2974 0.027650 1.3 6.75 2.011000 1.3860 -0.20570 1.0000 -0.3503 0.037770 1.6 6.75 1.753000 1.5640 -0.16780 1.0540 -0.3849 0.044300 2 6.65 1.251000 1.7480 -0.13160 1.1920 -0.4353 0.053610 2.5 6.70 0.930800 1.8240 -0.10870 1.2880 -0.4787 0.061430 3 6.65 0.515600 1.9080 -0.08979 1.4180 -0.5129 0.067610 4 6.85 0.343900 1.9310 -0.07471 1.5060 -0.5515 0.074290 5 6.85 -0.079230 1.9800 -0.06211 1.5850 -0.5800 0.078960 6.5 7.15 -0.006674 1.9730 -0.05453 1.6300 -0.5961 0.081170 8 7.50 0.256100 1.9420 -0.05230 1.5930 -0.6090 0.082990 10 7.45 -0.276300 1.9720 -0.04633 1.7230 -0.6196 0.084200 PGA 5.85 2.216000 0.6859 -0.13920 0.7656 -0.6187 0.060290 PGV 5.90 5.960000 1.0300 -0.16510 1.0790 -0.5785 0.057370 """) TAB3 = CoeffsTable(sa_damping=5, table="""\ imt s0 s1 s2 s3 s4 s5 s6 s7 s8 s9 0.01 -2.0480 1.8810 -0.49010 0.056680 -0.002433 -1.437000 1.24200 -0.28920 0.030880 -0.001252 0.013 -1.9220 1.8020 -0.47130 0.054710 -0.002357 -1.348000 1.19500 -0.27990 0.030060 -0.001225 0.016 -1.7110 1.6630 -0.43650 0.050870 -0.002199 -1.079000 1.04100 -0.24660 0.026900 -0.001114 0.02 -1.1600 1.2740 -0.33440 0.039110 -0.001700 -1.272000 1.25400 -0.31710 0.036240 -0.001550 0.025 -1.5350 1.5950 -0.42930 0.051030 -0.002242 -1.454000 1.36600 -0.33720 0.037270 -0.001537 0.03 -1.0560 1.2050 -0.31320 0.036150 -0.001550 -2.243000 1.98100 -0.50860 0.057820 -0.002439 0.04 -0.8571 1.0440 -0.26770 0.030820 -0.001328 -3.310000 2.66300 -0.66830 0.074150 -0.003056 0.05 -0.9628 0.9826 -0.21560 0.020800 -0.000742 -4.228000 3.29300 -0.83160 0.093030 -0.003873 0.065 -2.2250 1.9480 -0.49000 0.054860 -0.002293 -3.960000 2.87100 -0.66750 0.068830 -0.002650 0.08 -3.6850 2.9620 -0.75100 0.084210 -0.003509 -3.139000 2.17700 -0.46740 0.044660 -0.001598 0.1 -4.0510 3.1000 -0.76250 0.083280 -0.003393 -2.452000 1.56900 -0.28900 0.023000 -0.000657 0.13 -4.1740 3.0920 -0.74380 0.079820 -0.003205 -1.384000 0.62640 -0.01161 -0.010920 0.000828 0.16 -3.9650 2.8200 -0.64990 0.067200 -0.002614 -0.199700 -0.33700 0.25700 -0.042520 0.002176 0.2 -2.7070 1.7290 -0.33020 0.028160 -0.000906 0.819700 -1.08300 0.43950 -0.061050 0.002846 0.25 -1.7670 0.9826 -0.13140 0.005998 -0.000012 1.780000 -1.76700 0.60660 -0.078340 0.003498 0.3 -0.3182 -0.1386 0.17040 -0.028500 0.001421 2.245000 -2.00300 0.63260 -0.076990 0.003268 0.4 2.0180 -1.8570 0.61170 -0.076740 0.003341 2.422000 -1.93800 0.55580 -0.061740 0.002390 0.5 3.9560 -3.2880 0.98850 -0.119600 0.005142 0.855500 -0.45280 0.06459 0.005220 -0.000830 0.65 3.6450 -2.8220 0.79320 -0.089260 0.003555 -0.667100 0.92770 -0.37080 0.061830 -0.003430 0.8 2.4040 -1.6520 0.40880 -0.037100 0.001051 -2.124000 2.15200 -0.73010 0.105300 -0.005287 1 1.0660 -0.4552 0.03739 0.010330 -0.001084 -4.473000 4.05100 -1.27400 0.171000 -0.008137 1.3 -2.5080 2.5230 -0.84460 0.120500 -0.006024 -5.494000 4.76600 -1.43900 0.184900 -0.008458 1.6 -5.2640 4.7380 -1.47600 0.196300 -0.009284 -5.880000 4.97800 -1.46500 0.183200 -0.008156 2 -6.6420 5.7670 -1.74200 0.224100 -0.010280 -6.010000 4.98500 -1.43300 0.174800 -0.007587 2.5 -8.0780 6.8350 -2.01900 0.253800 -0.011410 -4.883000 3.94700 -1.08900 0.126400 -0.005173 3 -7.9800 6.6430 -1.92400 0.236600 -0.010390 -4.180000 3.31700 -0.88620 0.098880 -0.003853 4 -7.1230 5.7770 -1.61500 0.190400 -0.007982 -2.627000 1.96300 -0.46160 0.042400 -0.001176 5 -6.3880 5.0830 -1.38100 0.157600 -0.006361 -1.377000 0.90930 -0.14220 0.001323 0.000710 6.5 -4.8000 3.6840 -0.93680 0.097630 -0.003474 -0.392900 0.09826 0.09528 -0.027780 0.001959 8 -3.4160 2.5120 -0.58030 0.051530 -0.001344 -0.006872 -0.18880 0.16920 -0.035340 0.002200 10 -2.1940 1.5110 -0.28710 0.015340 0.000238 0.268400 -0.38620 0.21690 -0.039670 0.002304 PGA -2.1320 1.9370 -0.50400 0.058240 -0.002498 -1.444000 1.23500 -0.28510 0.030210 -0.001217 PGV -2.2460 1.9510 -0.51810 0.061390 -0.002725 -1.758000 1.37900 -0.32560 0.035000 -0.001425 """) TAB4 = CoeffsTable(sa_damping=5, table="""\ imt gCENA gCalifornia 0.01 -0.004661 -0.009823 0.013 -0.004693 -0.009833 0.016 -0.004687 -0.009834 0.02 -0.004668 -0.009817 0.025 -0.004884 -0.009881 0.03 -0.005113 -0.010140 0.04 -0.005266 -0.010760 0.05 -0.005471 -0.011320 0.065 -0.005714 -0.011940 0.08 -0.005794 -0.012390 0.1 -0.005640 -0.012500 0.13 -0.005236 -0.012170 0.16 -0.004771 -0.011660 0.2 -0.004203 -0.010920 0.25 -0.003648 -0.010190 0.3 -0.003121 -0.009435 0.4 -0.002438 -0.008259 0.5 -0.002041 -0.007363 0.65 -0.001638 -0.006452 0.8 -0.001426 -0.005849 1 -0.001259 -0.005130 1.3 -0.001063 -0.004346 1.6 -0.001171 -0.003900 2 -0.001016 -0.003359 2.5 -0.001060 -0.003012 3 -0.001093 -0.002725 4 -0.001301 -0.002123 5 -0.000935 -0.001698 6.5 -0.000787 -0.001306 8 -0.000643 -0.001061 10 -0.000365 -0.000849 PGA -0.004667 -0.009808 PGV -0.002792 -0.006313 """)
[docs]def get_fs_SeyhanStewart2014(imt, pga_rock, vs30): """ Implements eq. 11 and 12 at page 1992 in Yenier and Atkinson (2015) :param pga_rock: Median peak ground horizontal acceleration for reference :param vs30: """ # Get coefficients gmm = BooreEtAl2014() C = gmm.COEFFS[imt] # Linear term flin = vs30 / gmm.CONSTS['Vref'] flin[vs30 > C['Vc']] = C['Vc'] / gmm.CONSTS['Vref'] fl = C['c'] * np.log(flin) # Non-linear term v_s = np.copy(vs30) v_s[vs30 > 760.] = 760. # parameter (equation 8 of BSSA 2014) f_2 = C['f4'] * (np.exp(C['f5'] * (v_s - 360.)) - np.exp(C['f5'] * 400.)) fnl = gmm.CONSTS['f1'] + f_2 * np.log((pga_rock + gmm.CONSTS['f3']) / gmm.CONSTS['f3']) return fl + fnl