Source code for openquake.hazardlib.gsim.campbell_bozorgnia_2014
# -*- coding: utf-8 -*-# vim: tabstop=4 shiftwidth=4 softtabstop=4## Copyright (C) 2014-2025 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:`CampbellBozorgnia2014` :class:`CampbellBozorgnia2014HighQ` :class:`CampbellBozorgnia2014LowQ` :class:`CampbellBozorgnia2019` :class:`CampbellBozorgnia2019HighQ` :class:`CampbellBozorgnia2019LowQ`"""importnumpyasnpfromnumpyimportexp,radians,cosfromopenquake.hazardlib.gsim.baseimportGMPE,CoeffsTable,add_aliasfromopenquake.hazardlib.gsim.abrahamson_2014importget_epistemic_sigmafromopenquake.hazardlibimportconstfromopenquake.hazardlib.imtimportPGA,PGV,SA,IA,CAVfromopenquake.hazardlib.gsim.utils_usgs_basin_scalingimport \
_get_z2pt5_usgs_basin_scalingCONSTS={"h4":1.0,"c":1.88,"n":1.18}# CyberShake basin adjustments for CB14 (only applied above # 1.9 seconds so don't provide dummy values listed below 2 s)# Taken from https://code.usgs.gov/ghsc/nshmp/nshmp-lib/-/blob/main/src/main/resources/gmm/coeffs/CB14.csv?ref_type=headsCOEFFS_CY=CoeffsTable(sa_damping=5,table="""\IMT slope_cy2.0 0.7643.0 1.2794.0 1.7265.0 2.080 7.5 3.00010 3.000""")def_get_alpha(C,vs30,pga_rock):""" Returns the alpha, the linearised functional relationship between the site amplification and the PGA on rock. Equation 31. """alpha=np.zeros(len(pga_rock))idx=vs30<C["k1"]ifnp.any(idx):af1=pga_rock[idx]+\
CONSTS["c"]*((vs30[idx]/C["k1"])**CONSTS["n"])af2=pga_rock[idx]+CONSTS["c"]alpha[idx]=C["k2"]*pga_rock[idx]*((1.0/af1)-(1.0/af2))returnalphadef_get_anelastic_attenuation_term(C,rrup):""" Returns the anelastic attenuation term defined in equation 25 """f_atn=np.zeros(len(rrup))idx=rrup>=80.0f_atn[idx]=(C["c20"]+C["Dc20"])*(rrup[idx]-80.0)returnf_atndef_basin_term(C,imt,z2pt5,SJ,cy):""" Calculate the basin term """# Get CyberShake adjustment if required and SA(T > 1.9)ifcyandimt.period>1.9:basin_coeffs=COEFFS_CY[imt]['slope_cy']# Otherwise use regular coefficientselse:basin_coeffs=C["c16"]*C["k3"]fb=np.zeros(len(z2pt5))idx=z2pt5<1.0fb[idx]=(C["c14"]+C["c15"]*SJ)*(z2pt5[idx]-1.0)idx=z2pt5>3.0fb[idx]=basin_coeffs*exp(-0.75)*(1.-np.exp(-0.25*(z2pt5[idx]-3.)))returnfbdef_select_basin_model(SJ,vs30):""" Select the preferred basin model (California or Japan) to scale basin depth with respect to Vs30 """ifSJ:# Japan Basin Model - Equation 34 of Campbell & Bozorgnia (2014)returnnp.exp(5.359-1.102*np.log(vs30))else:# California Basin Model - Equation 33 of# Campbell & Bozorgnia (2014)returnnp.exp(7.089-1.144*np.log(vs30))def_get_basin_term(C,ctx,region,imt,SJ,a1100,usgs_bs=False,cy=False):""" Returns the basin response term defined in equation 20 and apply any required adjustments. """# Get reference basin depthz_ref=_select_basin_model(SJ,1100.0)*np.ones_like(ctx.vs30)z_ref_term=_basin_term(C,imt,z_ref,SJ,False)# Get basin termifisinstance(a1100,np.ndarray):# Site model definedz2pt5=ctx.z2pt5else:z2pt5=z_refz2pt5_term=_basin_term(C,imt,z2pt5,SJ,cy)# Apply USGS basin scaling model if requiredifusgs_bs:# Get the scaling factor per siteusgs_baf=_get_z2pt5_usgs_basin_scaling(ctx.z2pt5,imt.period)z_scaled=z_ref_term*(1.0-usgs_baf)+z2pt5_term*usgs_baf# Apply additional CyberShake (CY_CSIM) adjustment if requiredifcyandimt.period>1.9:z_scaled+=0.1# If period is 0.075 seconds scale by factor of 0.585ifimt.period==0.075:returnz_scaled*0.585# Otherwise return basin term adjusted except for sites shallower than# upper z2pt5 value in USGS basin model (use z_ref_term instead here)ifimt.period>0.5:idx=z2pt5>1.0# Upper z2pt5 depth in USGS basin model is 1 kmz_scaled[idx]=z_ref_term[idx]returnz_scaledreturnz2pt5_termdef_get_f1rx(C,r_x,r_1):""" Defines the f1 scaling coefficient defined in equation 9 """rxr1=r_x/r_1returnC["h1"]+C["h2"]*rxr1+C["h3"]*rxr1**2def_get_f2rx(C,r_x,r_1,r_2):""" Defines the f2 scaling coefficient defined in equation 10 """drx=(r_x-r_1)/(r_2-r_1)returnCONSTS["h4"]+C["h5"]*drx+C["h6"]*drx**2def_get_fault_dip_term(C,ctx):""" Returns the fault dip term, defined in equation 24 """res=C["c19"]*(5.5-ctx.mag)*ctx.dipres[ctx.mag<4.5]=C["c19"]*ctx.dip[ctx.mag<4.5]res[ctx.mag>5.5]=0.0returnresdef_get_geometric_attenuation_term(C,mag,rrup):""" Returns the geometric attenuation term defined in equation 3 """return(C["c5"]+C["c6"]*mag)*np.log(np.sqrt(rrup**2+C["c7"]**2))def_get_hanging_wall_coeffs_dip(dip):""" Returns the hanging wall dip term defined in equation 16 """return(90.0-dip)/45.0def_get_hanging_wall_coeffs_mag(C,mag):""" Returns the hanging wall magnitude term defined in equation 14 """res=(mag-5.5)*(1.0+C["a2"]*(mag-6.5))res[mag<5.5]=0.0res[mag>6.5]=1.0+C["a2"]*(mag[mag>6.5]-6.5)returnresdef_get_hanging_wall_coeffs_rrup(ctx):""" Returns the hanging wall rrup term defined in equation 13 """fhngrrup=np.ones(len(ctx.rrup))idx=ctx.rrup>0.0fhngrrup[idx]=(ctx.rrup[idx]-ctx.rjb[idx])/ctx.rrup[idx]returnfhngrrupdef_get_hanging_wall_coeffs_rx(C,ctx):""" Returns the hanging wall r-x caling term defined in equation 7 to 12 """r_x=ctx.rx# Define coefficients R1 and R2r_1=ctx.width*cos(radians(ctx.dip))r_2=62.0*ctx.mag-350.0fhngrx=np.zeros(len(r_x))# Case when 0 <= Rx <= R1idx=np.logical_and(r_x>=0.,r_x<r_1)fhngrx[idx]=_get_f1rx(C,r_x[idx],r_1[idx])# Case when Rx > R1idx=r_x>=r_1f2rx=_get_f2rx(C,r_x[idx],r_1[idx],r_2[idx])f2rx[f2rx<0.0]=0.0fhngrx[idx]=f2rxreturnfhngrxdef_get_hanging_wall_coeffs_ztor(ztor):""" Returns the hanging wall ztor term defined in equation 15 """res=1.-0.06*ztorres[ztor>16.66]=0.returnresdef_get_hanging_wall_term(C,ctx):""" Returns the hanging wall scaling term defined in equations 7 to 16 """return(C["c10"]*_get_hanging_wall_coeffs_rx(C,ctx)*_get_hanging_wall_coeffs_rrup(ctx)*_get_hanging_wall_coeffs_mag(C,ctx.mag)*_get_hanging_wall_coeffs_ztor(ctx.ztor)*_get_hanging_wall_coeffs_dip(ctx.dip))def_get_hypocentral_depth_term(C,ctx):""" Returns the hypocentral depth scaling term defined in equations 21 - 23 """fhyp_h=np.clip(ctx.hypo_depth-7.0,0.,13.)fhyp_m=C["c17"]+(C["c18"]-C["c17"])*(ctx.mag-5.5)fhyp_m[ctx.mag<=5.5]=C["c17"]fhyp_m[ctx.mag>6.5]=C["c18"]returnfhyp_h*fhyp_mdef_get_magnitude_term(C,mag):""" Returns the magnitude scaling term defined in equation 2 """f_mag=C["c0"]+C["c1"]*magaround5=(mag>4.5)&(mag<=5.5)around6=(mag>5.5)&(mag<=6.5)beyond=mag>6.5f_mag[around5]+=C["c2"]*(mag[around5]-4.5)f_mag[around6]+=(C["c2"]*(mag[around6]-4.5)+C["c3"]*(mag[around6]-5.5))f_mag[beyond]+=(C["c2"]*(mag[beyond]-4.5)+C["c3"]*(mag[beyond]-5.5)+C["c4"]*(mag[beyond]-6.5))returnf_magdef_get_philny(C,mag):""" Returns the intra-event random effects coefficient (phi) Equation 28. """res=C["phi2"]+(C["phi1"]-C["phi2"])*(5.5-mag)res[mag<=4.5]=C["phi1"]res[mag>=5.5]=C["phi2"]returnresdef_get_shallow_site_response_term(SJ,C,vs30,pga_rock):""" Returns the shallow site response term defined in equations 17, 18 and 19 """vs_mod=vs30/C["k1"]# Get linear global site response termf_site_g=C["c11"]*np.log(vs_mod)idx=vs30>C["k1"]f_site_g[idx]=f_site_g[idx]+(C["k2"]*CONSTS["n"]*np.log(vs_mod[idx]))# Get nonlinear site response termidx=np.logical_not(idx)ifnp.any(idx):f_site_g[idx]=f_site_g[idx]+C["k2"]*(np.log(pga_rock[idx]+CONSTS["c"]*(vs_mod[idx]**CONSTS["n"]))-np.log(pga_rock[idx]+CONSTS["c"]))ifSJ:fsite_j=(C["c13"]+C["k2"]*CONSTS["n"])* \
np.log(vs_mod)# additional term activated for soft ctx (Vs30 <= 200m/s)# in Japan dataidx=vs30<=200.0add_soft=(C["c12"]+C["k2"]*CONSTS["n"])* \
(np.log(vs_mod)-np.log(200.0/C["k1"]))# combine termsfsite_j[idx]+=add_soft[idx]returnf_site_g+fsite_jelse:returnf_site_gdef_get_style_of_faulting_term(C,ctx):""" Returns the style-of-faulting scaling term defined in equations 4 to 6 """frv=np.zeros_like(ctx.rake)fnm=np.zeros_like(ctx.rake)frv[(ctx.rake>30.)&(ctx.rake<150.)]=1.fnm[(ctx.rake>-150.)&(ctx.rake<-30.)]=1.fflt_f=C["c8"]*frv+C["c9"]*fnmfflt_m=ctx.mag-4.5fflt_m[ctx.mag<=4.5]=0.fflt_m[ctx.mag>5.5]=1.returnfflt_f*fflt_mdef_get_taulny(C,mag):""" Returns the inter-event random effects coefficient (tau) Equation 28. """res=C["tau2"]+(C["tau1"]-C["tau2"])*(5.5-mag)res[mag<=4.5]=C["tau1"]res[mag>=5.5]=C["tau2"]returnres
[docs]defget_mean_values(SJ,C,ctx,imt,usgs_bs=False,cy=False,a1100=None):""" Returns the mean values for a specific IMT """ifisinstance(a1100,np.ndarray):# Site model definedtemp_vs30=ctx.vs30else:# Default site and basin modeltemp_vs30=1100.0*np.ones(len(ctx))return(_get_magnitude_term(C,ctx.mag)+_get_geometric_attenuation_term(C,ctx.mag,ctx.rrup)+_get_style_of_faulting_term(C,ctx)+_get_hanging_wall_term(C,ctx)+_get_shallow_site_response_term(SJ,C,temp_vs30,a1100)+_get_basin_term(C,ctx,None,imt,SJ,a1100,usgs_bs,cy)+_get_hypocentral_depth_term(C,ctx)+_get_fault_dip_term(C,ctx)+_get_anelastic_attenuation_term(C,ctx.rrup))
def_update_ctx(gsim,ctx):""" Use the ztor, width and hypo_depth formula to estimate if the estimate attribute is set. """ifgsim.estimate_ztor:# Equation 4 and 5 of Chiou & Youngs 2014ctx.ztor=np.where((ctx.rake>30.)&(ctx.rake<150.),np.maximum(2.704-1.226*np.maximum(ctx.mag-5.849,0),0)**2,np.maximum(2.673-1.136*np.maximum(ctx.mag-4.970,0),0)**2)ifgsim.estimate_width:# width estimation requires Zbot# where Zbot is the depth to the bottom of the seismogenic crustifnothasattr(ctx,"zbot"):raiseKeyError('Zbot is required if width is unknown.')# Equation 39 of Campbell & Bozorgnia 2014mask=np.absolute(np.sin(np.radians(ctx.dip)))>0ctx.width=np.sqrt(10**((ctx.mag-4.07)/0.98))ctx.width[mask]=np.minimum(ctx.width[mask],(ctx.zbot[mask]-ctx.ztor[mask])/np.sin(np.radians(ctx.dip[mask])))ifgsim.estimate_hypo_depth:# Equation 36 of Campbell & Bozorgnia 2014fdz_m=np.where(ctx.mag<6.75,-4.317+0.984*ctx.mag,2.325)# Equation 37 of Campbell & Bozorgnia 2014fdz_d=np.where(ctx.dip<=40,0.0445*(ctx.dip-40),0)# The depth to the bottom of the rupture planezbor=ctx.ztor+ctx.width*np.sin(np.radians(ctx.dip))# Equation 35 of Campbell & Bozorgnia 2014mask=zbor>ctx.ztordz=np.zeros_like(ctx.ztor)dz[mask]=np.exp(np.minimum(fdz_m[mask]+fdz_d[mask],np.log(0.9*(zbor[mask]-ctx.ztor[mask]))))ctx.hypo_depth=ctx.ztor+dzdef_get_rholnpga(C,mag):""" Returns the magnitude-dependent correlation coefficient (rho) — Equation 5 of CB19. """rho_ln_pga=C["rho2pga"]+(C["rho1pga"]-C["rho2pga"])*(5.5-mag)rho_ln_pga[mag<=4.5]=C["rho1pga"]rho_ln_pga[mag>=5.5]=C["rho2pga"]returnrho_ln_pga
[docs]classCampbellBozorgnia2014(GMPE):""" Implements NGA-West 2 GMPE developed by Kenneth W. Campbell and Yousef Bozorgnia, published as "NGA-West2 Ground Motion Model for the Average Horizontal Components of PGA, PGV, and 5 % Damped Linear Acceleration Response Spectra" (2014, Earthquake Spectra, Volume 30, Number 3, pages 1087 - 1115). """#: Supported tectonic region type is active shallow crustDEFINED_FOR_TECTONIC_REGION_TYPE=const.TRT.ACTIVE_SHALLOW_CRUST#: Supported intensity measure types are spectral acceleration, peak#: ground velocity and peak ground accelerationDEFINED_FOR_INTENSITY_MEASURE_TYPES={PGA,PGV,SA,IA,CAV}#: 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, 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, rupture#: width and hypocentral depthREQUIRES_RUPTURE_PARAMETERS={'mag','rake','dip','ztor','width','hypo_depth'}#: Required distance measures are Rrup, Rjb and RxREQUIRES_DISTANCES={'rrup','rjb','rx'}def__init__(self,coeffs={},SJ=False,sigma_mu_epsilon=0.0,estimate_ztor=False,estimate_width=False,estimate_hypo_depth=False,usgs_basin_scaling=False,cybershake_basin_adj=False):# tested in logictree/case_71ifcoeffs:# extra coefficients by IMTself.COEFFS|=CoeffsTable.fromdict(coeffs)self.SJ=SJ# flag for Japanself.sigma_mu_epsilon=sigma_mu_epsilonself.estimate_ztor=estimate_ztorself.estimate_width=estimate_widthself.estimate_hypo_depth=estimate_hypo_depthifself.estimate_width:# To estimate a width, the GMPE needs Zbotself.REQUIRES_RUPTURE_PARAMETERS|={"zbot"}self.usgs_basin_scaling=usgs_basin_scalingself.cybershake_basin_adj=cybershake_basin_adj
[docs]defcompute(self,ctx:np.recarray,imts,mean,sig,tau,phi):""" See :meth:`superclass method <.base.GroundShakingIntensityModel.compute>` for spec of input and result values. """if(self.estimate_ztororself.estimate_widthorself.estimate_hypo_depth):ctx=ctx.copy()_update_ctx(self,ctx)C_PGA=self.COEFFS[PGA()]# Get mean and standard deviation of PGA on rock (Vs30 1100 m/s^2)pga1100=np.exp(get_mean_values(self.SJ,C_PGA,ctx,PGA(),self.usgs_basin_scaling,self.cybershake_basin_adj))form,imtinenumerate(imts):C=self.COEFFS[imt]# Get mean and standard deviations for IMTmean[m]=get_mean_values(self.SJ,C,ctx,imt,self.usgs_basin_scaling,self.cybershake_basin_adj,pga1100)mean[m]+=(self.sigma_mu_epsilon*get_epistemic_sigma(ctx))ifimt.string[:2]=="SA"andimt.period<0.25:# According to Campbell & Bozorgnia (2013) [NGA West 2 Report]# If Sa (T) < PGA for T < 0.25 then set mean Sa(T) to mean PGA# Get PGA on soilpga=get_mean_values(self.SJ,C_PGA,ctx,imt,self.usgs_basin_scaling,self.cybershake_basin_adj,pga1100)idx=mean[m]<=pgamean[m,idx]=pga[idx]mean[m]+=self.sigma_mu_epsilon*get_epistemic_sigma(ctx)# Get stddevs for PGA on basement rocktau_lnpga_b=_get_taulny(C_PGA,ctx.mag)phi_lnpga_b=np.sqrt(_get_philny(C_PGA,ctx.mag)**2.-C_PGA["philnAF"]**2.)# Get tau_lny on the basement rocktau_lnyb=_get_taulny(C,ctx.mag)# Get phi_lny on the basement rockphi_lnyb=np.sqrt(_get_philny(C,ctx.mag)**2.-C["philnAF"]**2.)# Get site scaling termalpha=_get_alpha(C,ctx.vs30,pga1100)ifimt.stringin['CAV','IA']:# Use formula in CB19 supplementary spreadsheett=np.sqrt(tau_lnyb**2+alpha**2*tau_lnpga_b**2+2.*alpha*_get_rholnpga(C,ctx.mag)*tau_lnyb*tau_lnpga_b)p=np.sqrt(phi_lnyb**2+C["philnAF"]**2+alpha**2*phi_lnpga_b**2+2.0*alpha*_get_rholnpga(C,ctx.mag)*phi_lnyb*phi_lnpga_b)else:# Use formula in CB14 supplementary spreadsheett=np.sqrt(tau_lnyb**2+alpha**2*tau_lnpga_b**2+2.0*alpha*C["rholny"]*tau_lnyb*tau_lnpga_b)p=np.sqrt(phi_lnyb**2+C["philnAF"]**2+alpha**2*phi_lnpga_b**2+2.0*alpha*C["rholny"]*phi_lnyb*phi_lnpga_b)sig[m]=np.sqrt(t**2+p**2)tau[m]=tphi[m]=p
