# -*- coding: utf-8 -*-
# vim: tabstop=4 shiftwidth=4 softtabstop=4
#
# Copyright (C) 2013-2026 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
from openquake.hazardlib.gsim.chiou_youngs_2008 import _get_z1_ref
METRES_PER_KM = 1000.
def _get_cb07_z2pt5_ref(vs30):
"""
Estimate unknown z2pt5 using the relationships described
in Campbell and Bozorgnia (2007). Returns z2pt5 in km.
"""
# First get z1pt0 (in metres) using Chiou and Youngs (2008)
z1pt0 = _get_z1_ref(vs30)
# Now use Campbell and Bozorgnia (2007) to get z2pt5 from z1pt0
z2pt5 = 519 + (3.595 * z1pt0)
return z2pt5 / METRES_PER_KM
def _get_basin_term(C, ctx, region=None):
"""
Returns the basin response term (equation 12, page 146)
"""
z2pt5 = ctx.z2pt5.copy()
# Use GMM's vs30 to z2pt5 for non-measured values
mask = z2pt5 == -999
z2pt5[mask] = _get_cb07_z2pt5_ref(ctx.vs30[mask])
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) +
_get_basin_term(C, ctx) +
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) +
_get_basin_term(C, ctx))
# 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
[docs]class CampbellBozorgnia2008BCHydro(CampbellBozorgnia2008):
"""
CB08 with BCHydro (IMT-dependent) sigma-mu scaling and
single station sigma.
"""
experimental = True
def __init__(self, single_stat_sigma=False, sigma_mu_epsilon=0., **kwargs):
super().__init__(**kwargs)
self.single_stat_sigma = single_stat_sigma
self.sigma_mu_epsilon = sigma_mu_epsilon
[docs] def compute(self, ctx: np.recarray, imts, mean, sig, tau, phi):
super().compute(ctx, imts, mean, sig, tau, phi)
for m, imt in enumerate(imts):
C = self.COEFFS_BCHYDRO[imt]
mean[m] += self.sigma_mu_epsilon * C['fe']
if self.single_stat_sigma:
sig[m] *= C['s2sss']
COEFFS_BCHYDRO = CoeffsTable(sa_damping=5, table="""\
imt s2sss fe
pga 0.8530 0.2500
pgv 0.8700 0.2500
pgd 0.8700 0.2500
cav 0.8700 0.2500
0.010 0.8530 0.2500
0.020 0.8555 0.2500
0.030 0.8570 0.2500
0.050 0.8589 0.2500
0.075 0.8604 0.2500
0.100 0.8615 0.2500
0.150 0.8630 0.2500
0.200 0.8640 0.2500
0.250 0.8649 0.2500
0.300 0.8655 0.2500
0.400 0.8666 0.2500
0.500 0.8674 0.2500
0.750 0.8689 0.2500
1.000 0.8700 0.2500
1.500 0.8715 0.2813
2.000 0.8726 0.3125
3.000 0.8741 0.3750
4.000 0.8751 0.4375
5.000 0.8760 0.5000
7.500 0.8775 0.5000
10.00 0.8785 0.5000
""")
