Source code for openquake.hazardlib.gsim.campbell_bozorgnia_2008

# -*- coding: utf-8 -*-
# vim: tabstop=4 shiftwidth=4 softtabstop=4
#
# Copyright (C) 2013-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:`CampbellBozorgnia2008`, and
:class:'CampbellBozorgnia2008Arbitrary'
"""
import numpy as np
from math import log, exp
from openquake.hazardlib.gsim.base import GMPE, CoeffsTable
from openquake.hazardlib import const
from openquake.hazardlib.imt import PGA, PGV, PGD, CAV, SA


def _compute_basin_response_term(C, z2pt5):
    """
    Returns the basin response term (equation 12, page 146)
    """
    fsed = np.zeros_like(z2pt5, dtype=float)
    idx = z2pt5 < 1.0
    if np.any(idx):
        fsed[idx] = C['c11'] * (z2pt5[idx] - 1.0)

    idx = z2pt5 > 3.0
    if np.any(idx):
        fsed[idx] = (C['c12'] * C['k3'] * exp(-0.75)) *\
            (1.0 - np.exp(-0.25 * (z2pt5[idx] - 3.0)))
    return fsed


def _compute_distance_term(C, ctx):
    """
    Returns the distance scaling factor (equation (3), page 145)
    """
    return (C['c4'] + C['c5'] * ctx.mag) * \
        np.log(np.sqrt(ctx.rrup ** 2. + C['c6'] ** 2.))


def _compute_hanging_wall_term(C, ctx):
    """
    Returns the hanging wall scaling term, the product of the scaling
    coefficient and four separate scaling terms for distance, magnitude,
    rupture depth and dip (equations 6 - 10, page 146). Individual
    scaling terms defined in separate functions
    """
    return (C['c9'] *
            _get_hanging_wall_distance_term(ctx) *
            _get_hanging_wall_magnitude_term(ctx.mag) *
            _get_hanging_wall_depth_term(ctx.ztor) *
            _get_hanging_wall_dip_term(ctx.dip))


def _compute_imt1100(C, ctx, get_pga_site=False):
    """
    Computes the PGA on reference (Vs30 = 1100 m/s) rock.
    """
    # Calculates simple site response term assuming all sites 1100 m/s
    fsite = (C['c10'] + (C['k2'] * C['n'])) * log(1100. / C['k1'])
    # Calculates the PGA on rock
    pga1100 = np.exp(_compute_magnitude_term(C, ctx.mag) +
                     _compute_distance_term(C, ctx) +
                     _compute_style_of_faulting_term(C, ctx) +
                     _compute_hanging_wall_term(C, ctx) +
                     _compute_basin_response_term(C, ctx.z2pt5) +
                     fsite)
    # If PGA at the site is needed then remove factor for rock and
    # re-calculate on correct site condition
    if get_pga_site:
        pga_site = np.exp(np.log(pga1100) - fsite)
        fsite = _compute_shallow_site_response(C, ctx, pga1100)
        pga_site = np.exp(np.log(pga_site) + fsite)
    else:
        pga_site = None
    return pga1100, pga_site


def _compute_intra_event_alpha(C, vs30, pga1100):
    """
    Returns the linearised functional relationship between fsite and
    pga1100, determined from the partial derivative defined on equation 17
    on page 148
    """
    alpha = np.zeros_like(vs30, dtype=float)
    idx = vs30 < C['k1']
    if np.any(idx):
        temp1 = (pga1100[idx] +
                 C['c'] * (vs30[idx] / C['k1']) ** C['n']) ** -1.
        temp1 = temp1 - ((pga1100[idx] + C['c']) ** -1.)
        alpha[idx] = C['k2'] * pga1100[idx] * temp1

    return alpha


def _compute_intra_event_std(C, vs30, pga1100, sigma_pga):
    """
    Returns the intra-event standard deviation at the site, as defined in
    equation 15, page 147
    """
    # Get intra-event standard deviation at the base of the site profile
    sig_lnyb = np.sqrt(C['s_lny'] ** 2. - C['s_lnAF'] ** 2.)
    sig_lnab = np.sqrt(sigma_pga ** 2. - C['s_lnAF'] ** 2.)
    # Get linearised relationship between f_site and ln PGA
    alpha = _compute_intra_event_alpha(C, vs30, pga1100)

    return np.sqrt(
        (sig_lnyb ** 2.) +
        (C['s_lnAF'] ** 2.) +
        ((alpha ** 2.) * (sig_lnab ** 2.)) +
        (2.0 * alpha * C['rho'] * sig_lnyb * sig_lnab))


def _compute_magnitude_term(C, mag):
    """
    Returns the magnitude scaling factor (equation (2), page 144)
    """
    fmag = C['c0'] + C['c1'] * mag
    term = C['c2'] * (mag - 5.5)
    term[mag <= 5.5] = 0.
    term[mag > 6.5] = C['c2'] * (mag[mag > 6.5] - 5.5) + (
        C['c3'] * (mag[mag > 6.5] - 6.5))
    return fmag + term


def _compute_shallow_site_response(C, ctx, pga1100):
    """
    Returns the shallow site response term (equation 11, page 146)
    """
    stiff_factor = C['c10'] + (C['k2'] * C['n'])
    # Initially default all sites to intermediate rock value
    fsite = stiff_factor * np.log(ctx.vs30 / C['k1'])
    # Check for soft soil ctx
    idx = ctx.vs30 < C['k1']
    if np.any(idx):
        pga_scale = np.log(pga1100[idx] +
                           (C['c'] * ((ctx.vs30[idx] / C['k1']) **
                            C['n']))) - np.log(pga1100[idx] + C['c'])
        fsite[idx] = C['c10'] * np.log(ctx.vs30[idx] / C['k1']) + \
            (C['k2'] * pga_scale)
    # Any very hard rock ctx are rendered to the constant amplification
    # factor
    idx = ctx.vs30 >= 1100.
    if np.any(idx):
        fsite[idx] = stiff_factor * log(1100. / C['k1'])

    return fsite


def _compute_style_of_faulting_term(C, ctx):
    """
    Returns the style of faulting factor, depending on the mechanism (rake)
    and top of rupture depth (equations (4) and (5), pages 145 - 146)
    """
    frv, fnm = _get_fault_type_dummy_variables(ctx.rake)
    ffltz = np.zeros_like(ctx.rake)
    # Top of rupture depth term only applies to reverse faults
    ffltz[(frv > 0.) & (ctx.ztor < 1)] = ctx.ztor[(frv > 0.) & (ctx.ztor < 1)]
    ffltz[(frv > 0.) & (ctx.ztor >= 1)] = 1.
    return C['c7'] * frv * ffltz + C['c8'] * fnm


def _get_fault_type_dummy_variables(rake):
    """
    Returns the coefficients FRV and FNM, describing if the rupture is
    reverse (FRV = 1.0, FNM = 0.0), normal (FRV = 0.0, FNM = 1.0) or
    strike-slip/oblique-slip (FRV = 0.0, FNM = 0.0). Reverse faults are
    classified as those with a rake in the range 30 to 150 degrees. Normal
    faults are classified as having a rake in the range -150 to -30 degrees
    :returns:
        FRV, FNM
    """
    frv, fnm = np.zeros_like(rake), np.zeros_like(rake)
    frv[(rake > 30.0) & (rake < 150.)] = 1.
    fnm[(rake > -150.) & (rake < -30.)] = 1.
    return frv, fnm


def _get_hanging_wall_depth_term(ztor):
    """
    Returns the hanging wall depth scaling term (equation 9, page 146)
    """
    return np.where(ztor >= 20.0, 0., (20. - ztor) / 20.0)


def _get_hanging_wall_dip_term(dip):
    """
    Returns the hanging wall dip scaling term (equation 10, page 146)
    """
    return np.where(dip > 70.0, (90.0 - dip) / 20.0, 1.0)


def _get_hanging_wall_distance_term(ctx):
    """
    Returns the hanging wall distance scaling term (equation 7, page 146)
    """
    fhngr = np.ones_like(ctx.rjb, dtype=float)
    idx = (ctx.rjb > 0.) & (ctx.ztor < 1)
    temp_rjb = np.sqrt(ctx.rjb[idx] ** 2. + 1.)
    r_max = np.max(np.column_stack([ctx.rrup[idx], temp_rjb]), axis=1)
    fhngr[idx] = (r_max - ctx.rjb[idx]) / r_max
    idx = (ctx.rjb > 0.) & (ctx.ztor >= 1)
    fhngr[idx] = (ctx.rrup[idx] - ctx.rjb[idx]) / ctx.rrup[idx]
    return fhngr


def _get_hanging_wall_magnitude_term(mag):
    """
    Returns the hanging wall magnitude scaling term (equation 8, page 146)
    """
    return np.clip(2. * (mag - 6.0), 0., 1.)


def _get_stddevs(kind, C, ctx, pga1100, sigma_pga):
    """
    Returns the standard deviations as described in the "ALEATORY
    UNCERTAINTY MODEL" section of the paper. Equations 13 to 19, pages 147
    to 151
    """
    std_intra = _compute_intra_event_std(C, ctx.vs30, pga1100, sigma_pga)
    std_inter = C['t_lny'] * np.ones_like(ctx.vs30)
    return [_get_total_sigma(kind, C, std_intra, std_inter),
            std_inter, std_intra]


def _get_total_sigma(kind, C, std_intra, std_inter):
    """
    Returns the total sigma term as defined by equation 16, page 147
    This method is defined here as the Campbell & Bozorgnia (2008) model
    can also be applied to the "arbitrary" horizontal component
    definition, in which case the total sigma is modified (see
    equation 18, page 150).
    """
    tot2 = std_intra ** 2. + std_inter ** 2.
    if kind == 'arbitrary':
        return np.sqrt(tot2 + C['c_lny'] ** 2.)
    return np.sqrt(tot2)


[docs]class CampbellBozorgnia2008(GMPE): """ Implements GMPE developed by Kenneth W. Campbell and Yousef Bozorgnia, published as "NGA Ground Motion Model for the Geometric Mean Horizontal Component of PGA, PGV, PGD and 5 % Damped Linear Elastic Response Spectra for Periods Ranging from 0.01 to 10s" (2008, Earthquake Spectra, Volume 24, Number 1, pages 139 - 171). This class implements the model for the Geometric Mean of the elastic spectra. Included in the coefficient set are the coefficients for the Campbell & Bozorgnia (2010) GMPE for predicting Cumulative Absolute Velocity (CAV), published as "A Ground Motion Prediction Equation for the Horizontal Component of Cumulative Absolute Velocity (CSV) Based on the PEER-NGA Strong Motion Database" (2010, Earthquake Spectra, Volume 26, Number 3, 635 - 650). """ kind = "base" #: Supported tectonic region type is active shallow crust DEFINED_FOR_TECTONIC_REGION_TYPE = const.TRT.ACTIVE_SHALLOW_CRUST #: Supported intensity measure types are spectral acceleration, peak #: ground velocity, peak ground displacement and peak ground acceleration #: Additional model for cumulative absolute velocity defined in #: Campbell & Bozorgnia (2010) DEFINED_FOR_INTENSITY_MEASURE_TYPES = {PGA, PGV, PGD, CAV, SA} #: Supported intensity measure component is orientation-independent #: average horizontal :attr:`~openquake.hazardlib.const.IMC.GMRotI50` DEFINED_FOR_INTENSITY_MEASURE_COMPONENT = const.IMC.GMRotI50 #: Supported standard deviation types are inter-event, intra-event #: and total, see section "Aleatory Uncertainty Model", page 147. DEFINED_FOR_STANDARD_DEVIATION_TYPES = { const.StdDev.TOTAL, const.StdDev.INTER_EVENT, const.StdDev.INTRA_EVENT} #: Required site parameters are Vs30, Vs30 type (measured or inferred), #: and depth (km) to the 2.5 km/s shear wave velocity layer (z2pt5) REQUIRES_SITES_PARAMETERS = {'vs30', 'z2pt5'} #: Required rupture parameters are magnitude, rake, dip, ztor REQUIRES_RUPTURE_PARAMETERS = {'mag', 'rake', 'dip', 'ztor'} #: Required distance measures are Rrup and Rjb. REQUIRES_DISTANCES = {'rrup', 'rjb'}
[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()] for m, imt in enumerate(imts): C = self.COEFFS[imt] # compute median pga on rock (vs30=1100), needed for site response # term calculation # For spectral accelerations at periods between 0.0 and 0.25 s, # Sa (T) cannot be less than PGA on soil, therefore if the IMT is # in this period range it is necessary to calculate PGA on soil get_pga_site = imt.period > 0.0 and imt.period < 0.25 pga1100, pga_site = _compute_imt1100(C_PGA, ctx, get_pga_site) # Get the median ground motion mean[m] = (_compute_magnitude_term(C, ctx.mag) + _compute_distance_term(C, ctx) + _compute_style_of_faulting_term(C, ctx) + _compute_hanging_wall_term(C, ctx) + _compute_shallow_site_response(C, ctx, pga1100) + _compute_basin_response_term(C, ctx.z2pt5)) # If it is necessary to ensure that Sa(T) >= PGA # (see previous comment) if get_pga_site: idx = mean[m] < np.log(pga_site) mean[m, idx] = np.log(pga_site[idx]) sig[m], tau[m], phi[m] = _get_stddevs( self.kind, C, ctx, pga1100, C_PGA['s_lny'])
COEFFS = CoeffsTable(sa_damping=5, table="""\ imt c0 c1 c2 c3 c4 c5 c6 c7 c8 c9 c10 c11 c12 k1 k2 k3 c n s_lny t_lny s_lnAF c_lny rho cav -4.354 0.942 -0.178 -0.346 -1.309 0.087 7.24 0.111 -0.108 0.362 2.549 0.090 1.277 400 -2.690 1.000 1.88 1.18 0.371 0.196 0.300 0.089 0.735 pgd -5.270 1.600 -0.070 0.000 -2.000 0.17 4.00 0.000 0.000 0.000 -0.820 0.300 1.000 400 0.000 2.744 1.88 1.18 0.667 0.485 0.300 0.290 0.174 pgv 0.954 0.696 -0.309 -0.019 -2.016 0.17 4.00 0.245 0.000 0.358 1.694 0.092 1.000 400 -1.955 1.929 1.88 1.18 0.484 0.203 0.300 0.190 0.691 pga -1.715 0.500 -0.530 -0.262 -2.118 0.17 5.60 0.280 -0.120 0.490 1.058 0.040 0.610 865 -1.186 1.839 1.88 1.18 0.478 0.219 0.300 0.166 1.000 0.010 -1.715 0.500 -0.530 -0.262 -2.118 0.17 5.60 0.280 -0.120 0.490 1.058 0.040 0.610 865 -1.186 1.839 1.88 1.18 0.478 0.219 0.300 0.166 1.000 0.020 -1.680 0.500 -0.530 -0.262 -2.123 0.17 5.60 0.280 -0.120 0.490 1.102 0.040 0.610 865 -1.219 1.840 1.88 1.18 0.480 0.219 0.300 0.166 0.999 0.030 -1.552 0.500 -0.530 -0.262 -2.145 0.17 5.60 0.280 -0.120 0.490 1.174 0.040 0.610 908 -1.273 1.841 1.88 1.18 0.489 0.235 0.300 0.165 0.989 0.050 -1.209 0.500 -0.530 -0.267 -2.199 0.17 5.74 0.280 -0.120 0.490 1.272 0.040 0.610 1054 -1.346 1.843 1.88 1.18 0.510 0.258 0.300 0.162 0.963 0.075 -0.657 0.500 -0.530 -0.302 -2.277 0.17 7.09 0.280 -0.120 0.490 1.438 0.040 0.610 1086 -1.471 1.845 1.88 1.18 0.520 0.292 0.300 0.158 0.922 0.100 -0.314 0.500 -0.530 -0.324 -2.318 0.17 8.05 0.280 -0.099 0.490 1.604 0.040 0.610 1032 -1.624 1.847 1.88 1.18 0.531 0.286 0.300 0.170 0.898 0.150 -0.133 0.500 -0.530 -0.339 -2.309 0.17 8.79 0.280 -0.048 0.490 1.928 0.040 0.610 878 -1.931 1.852 1.88 1.18 0.532 0.280 0.300 0.180 0.890 0.200 -0.486 0.500 -0.446 -0.398 -2.220 0.17 7.60 0.280 -0.012 0.490 2.194 0.040 0.610 748 -2.188 1.856 1.88 1.18 0.534 0.249 0.300 0.186 0.871 0.250 -0.890 0.500 -0.362 -0.458 -2.146 0.17 6.58 0.280 0.000 0.490 2.351 0.040 0.700 654 -2.381 1.861 1.88 1.18 0.534 0.240 0.300 0.191 0.852 0.300 -1.171 0.500 -0.294 -0.511 -2.095 0.17 6.04 0.280 0.000 0.490 2.460 0.040 0.750 587 -2.518 1.865 1.88 1.18 0.544 0.215 0.300 0.198 0.831 0.400 -1.466 0.500 -0.186 -0.592 -2.066 0.17 5.30 0.280 0.000 0.490 2.587 0.040 0.850 503 -2.657 1.874 1.88 1.18 0.541 0.217 0.300 0.206 0.785 0.500 -2.569 0.656 -0.304 -0.536 -2.041 0.17 4.73 0.280 0.000 0.490 2.544 0.040 0.883 457 -2.669 1.883 1.88 1.18 0.550 0.214 0.300 0.208 0.735 0.750 -4.844 0.972 -0.578 -0.406 -2.000 0.17 4.00 0.280 0.000 0.490 2.133 0.077 1.000 410 -2.401 1.906 1.88 1.18 0.568 0.227 0.300 0.221 0.628 1.000 -6.406 1.196 -0.772 -0.314 -2.000 0.17 4.00 0.255 0.000 0.490 1.571 0.150 1.000 400 -1.955 1.929 1.88 1.18 0.568 0.255 0.300 0.225 0.534 1.500 -8.692 1.513 -1.046 -0.185 -2.000 0.17 4.00 0.161 0.000 0.490 0.406 0.253 1.000 400 -1.025 1.974 1.88 1.18 0.564 0.296 0.300 0.222 0.411 2.000 -9.701 1.600 -0.978 -0.236 -2.000 0.17 4.00 0.094 0.000 0.371 -0.456 0.300 1.000 400 -0.299 2.019 1.88 1.18 0.571 0.296 0.300 0.226 0.331 3.000 -10.556 1.600 -0.638 -0.491 -2.000 0.17 4.00 0.000 0.000 0.154 -0.820 0.300 1.000 400 0.000 2.110 1.88 1.18 0.558 0.326 0.300 0.229 0.289 4.000 -11.212 1.600 -0.316 -0.770 -2.000 0.17 4.00 0.000 0.000 0.000 -0.820 0.300 1.000 400 0.000 2.200 1.88 1.18 0.576 0.297 0.300 0.237 0.261 5.000 -11.684 1.600 -0.070 -0.986 -2.000 0.17 4.00 0.000 0.000 0.000 -0.820 0.300 1.000 400 0.000 2.291 1.88 1.18 0.601 0.359 0.300 0.237 0.200 7.500 -12.505 1.600 -0.070 -0.656 -2.000 0.17 4.00 0.000 0.000 0.000 -0.820 0.300 1.000 400 0.000 2.517 1.88 1.18 0.628 0.428 0.300 0.271 0.174 10.00 -13.087 1.600 -0.070 -0.422 -2.000 0.17 4.00 0.000 0.000 0.000 -0.820 0.300 1.000 400 0.000 2.744 1.88 1.18 0.667 0.485 0.300 0.290 0.174 """)
[docs]class CampbellBozorgnia2008Arbitrary(CampbellBozorgnia2008): """ Implements the Campbell & Bozorgnia (2008) GMPE as modified to represent the arbitrary horizontal component of ground motion, instead of the Rotationally Independent Geometric Mean (GMRotI) originally defined in the paper. """ kind = "arbitrary" #: Supported intensity measure component is arbitrary horizontal #: :attr:`~openquake.hazardlib.const.IMC.HORIZONTAL`, DEFINED_FOR_INTENSITY_MEASURE_COMPONENT = const.IMC.HORIZONTAL