Source code for openquake.hazardlib.gsim.gulerce_2017
# -*- 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:`GulerceEtAl2017` :class:`GulerceEtAl2017RegTWN` :class:`GulerceEtAl2017RegITA` :class:`GulerceEtAl2017RegMID` :class:`GulerceEtAl2017RegCHN` :class:`GulerceEtAl2017RegJPN`"""importnumpyasnpfromopenquake.baselib.generalimportCallableDictfromopenquake.hazardlibimportconstfromopenquake.hazardlib.gsim.baseimportCoeffsTable,GMPEfromopenquake.hazardlib.imtimportSA#: equation constants (that are IMT independent)CONSTS={# m1, m2 specified at section "Moderate-to-Large Magnitude Scaling"'m1':6.75,'m2':5.50,# h1, h2, h3 specified at section "Hanging Wall Effects"'h1':+0.25,'h2':+1.50,'h3':-0.75}def_get_basic_term(C,ctx):""" Compute and return basic form, see Equation 11 to 13. """# Fictitious depth calculation, Equation 13. Unlike ASK14, the break in# the c4m function is shifted to M6.0.# The equation for c4m for M4.0-6.0 is different from GKAS16 EQS paper,# but used the supplementary material instead after code verification.c4m=C['c4']-(C['c4']-1.)*(6.-ctx.mag)/2.c4m[ctx.mag>6.]=C['c4']c4m[ctx.mag<4.]=1.# Equation 12R=np.sqrt(ctx.rrup**2.+c4m**2.)# basic form, Equation 11base_term=C['a1']*np.ones_like(ctx.rrup)+C['a17']*ctx.rrup# NB: m1 > m2after=ctx.mag>=CONSTS['m1']within=(ctx.mag<CONSTS['m1'])&(ctx.mag>=CONSTS['m2'])before=ctx.mag<CONSTS['m2']base_term[after]+=(C['a5']*(ctx.mag-CONSTS['m1'])+C['a8']*(8.5-ctx.mag)**2.+(C['a2']+C['a3']*(ctx.mag-CONSTS['m1']))*np.log(R))[after]base_term[within]+=(C['a4']*(ctx.mag-CONSTS['m1'])+C['a8']*(8.5-ctx.mag)**2.+(C['a2']+C['a3']*(ctx.mag-CONSTS['m1']))*np.log(R))[within]base_term[before]+=(C['a4']*(CONSTS['m2']-CONSTS['m1'])+C['a8']*(8.5-CONSTS['m2'])**2.+C['a6']*(ctx.mag-CONSTS['m2'])+(C['a2']+C['a3']*(CONSTS['m2']-CONSTS['m1']))*np.log(R))[before]returnbase_termdef_get_faulting_style_term(C,ctx):""" Compute and return faulting style term, that is the sum of the second and third terms in Equation 1. """# this implements Equations 3 and 4;# f7 is the term for reverse fault mechanisms;# f8 is the term for normal fault mechanisms.f7=np.where(ctx.mag>5.,C['a11'],np.where(ctx.mag>=4.,C['a11']*(ctx.mag-4.),0.))f8=np.where(ctx.mag>5.,C['a12'],np.where(ctx.mag>=4.,C['a12']*(ctx.mag-4.),0.))# ranges of rake values for each faulting mechanism are same with ASK14return(f7*((ctx.rake>30)&(ctx.rake<150))+f8*((ctx.rake>-150)&(ctx.rake<-30)))def_get_hanging_wall_term(C,ctx):""" Compute and return hanging wall model term, see section on "Hanging Wall Effects". """Fhw=np.zeros_like(ctx.rx)Fhw[ctx.rx>0]=1.# Compute dip taper t1, Equation 6T1=np.ones_like(ctx.rx)T1*=np.where(ctx.dip<=30,60/45,(90.-ctx.dip)/45)# Compute magnitude taper t2, Equation 7, with a2hw set to 0.2.T2=np.zeros_like(ctx.rx)a2hw=0.2after=ctx.mag>=6.5within=(ctx.mag>5.5)&(ctx.mag<6.5)before=ctx.mag<=5.5T2[after]+=(1.+a2hw*(ctx.mag[after]-6.5))T2[within]+=(1.+a2hw*(ctx.mag[within]-6.5)-(1.-a2hw)*(ctx.mag[within]-6.5)**2)T2[before]=0.# Compute distance taper t3, Equation 8T3=np.zeros_like(ctx.rx)r1=ctx.width*np.cos(np.radians(ctx.dip))# The r2 term is different here from ASK14 where r2 = 3*r1.r2=4.*r1#idx=ctx.rx<r1T3[idx]=(np.ones_like(ctx.rx)[idx]*CONSTS['h1']+CONSTS['h2']*(ctx.rx[idx]/r1[idx])+CONSTS['h3']*(ctx.rx[idx]/r1[idx])**2)#idx=(ctx.rx>=r1)&(ctx.rx<=r2)T3[idx]=1.-(ctx.rx[idx]-r1[idx])/(r2[idx]-r1[idx])# Compute depth taper t4, Equation 9T4=np.zeros_like(ctx.rx)#T4[ctx.ztor<=10.]+=(1.-ctx.ztor[ctx.ztor<=10.]**2./100.)# Compute off-edge distance taper T5, Equation 10# ry1 computed same as in ASK14T5=np.zeros_like(ctx.rx)ry1=ctx.rx*np.tan(np.radians(20.))#idx=(ctx.ry0-ry1)<=0.0T5[idx]=1.#idx=(((ctx.ry0-ry1)>0.0)&((ctx.ry0-ry1)<5.0))T5[idx]=1.-(ctx.ry0[idx]-ry1[idx])/5.0# Finally, compute the hanging wall term, Equation 5Fhw[ctx.dip==90.0]=0.returnFhw*C['a13']*T1*T2*T3*T4*T5def_get_inter_event_std(region,C,mag):""" Returns inter-event standard deviation, Tau, Equation 20 """tau=_get_tau_regional(region,C,mag)returntaudef_get_intra_event_std(region,C,mag):""" Returns intra-event std dev, Phi, Equation 19. """# Intra-event standard deviation model is simplified since the effect# of nonlinearity of the rock motion is not incorporated# (Equations 27-30 in ASK14 are not used).phi=_get_phi_regional(region,C,mag)returnphi_get_phi_regional=CallableDict()@_get_phi_regional.add("CAL","CHN","ITA","TWN","MID")def_get_phi_regional_1(region,C,mag):""" Returns regional (default) intra-event standard deviation """phi=C['s1']+(C['s2_noJP']-C['s1'])/2.*(mag-4.)phi[mag<4]=C['s1']phi[mag>6]=C['s2_noJP']returnphi@_get_phi_regional.add("JPN")def_get_phi_regional_2(region,C,mag):""" Returns regional intra-event standard deviation (Phi) for Japan """phi=C['s1']+(C['s2_all']-C['s1'])/2.*(mag-4.)phi[mag<4]=C['s1']phi[mag>6]=C['s2_all']returnphi_get_regional_term=CallableDict()@_get_regional_term.add("CAL")def_get_regional_term_CAL(region,C,imt,vs30,rrup):""" As with ASK14, we assume California as the default region, hence here the regional term is assumed = 0. """return0.@_get_regional_term.add("TWN")def_get_regional_term_TWN(region,C,imt,vs30,rrup):""" Compute regional term for Taiwan, see section "Regionalization and Aftershocks" """vs30star=_get_vs30star(vs30,imt)returnC['a31']*np.log(vs30star/C['vlin'])+C['a25']*rrup@_get_regional_term.add("ITA")def_get_regional_term_ITA(region,C,imt,vs30,rrup):""" Compute regional term for Italy, see section "Regionalization and Aftershocks" """# removed regional linear vs30 scaling term since a32=0returnC['a26']*rrup@_get_regional_term.add("MID")def_get_regional_term_MID(region,C,imt,vs30,rrup):""" Compute regional term for Middle East, see section "Regionalization and Aftershocks" """returnC['a27']*rrup@_get_regional_term.add("CHN")def_get_regional_term_CHN(region,C,imt,vs30,rrup):""" Compute regional term for China, see section "Regionalization and Aftershocks" """returnC['a28']*rrup@_get_regional_term.add("JPN")def_get_regional_term_JPN(region,C,imt,vs30,rrup):""" Compute regional term for Japan, see section "Regionalization and Aftershocks" """vs30star=_get_vs30star(vs30,imt)returnC['a35']*np.log(vs30star/C['vlin'])+C['a29']*rrupdef_get_site_response_term(C,imt,vs30):""" Compute and return site response model term; see section "Site Amplification Effects". """# vs30 star, Equation 15vs30_star=_get_vs30star(vs30,imt)# compute the site termsite_resp_term=np.zeros_like(vs30)# Unlike ASK14, the site term here is independent of nonlinear response# parameters b, c, and n.vs30_rat=vs30_star/C['vlin']site_resp_term=C['a10']*np.log(vs30_rat)returnsite_resp_termdef_get_stddevs(region,C,imt,ctx):""" Return standard deviations as described in section "Equations for Standard Deviation". """std_intra=_get_intra_event_std(region,C,ctx.mag)std_inter=_get_inter_event_std(region,C,ctx.mag)return[np.sqrt(std_intra**2+std_inter**2),std_inter,std_intra]_get_tau_regional=CallableDict()@_get_tau_regional.add("CAL","CHN","ITA","TWN","MID")def_get_tau_regional_CAL(region,C,mag):""" Returns regional (default) inter-event standard deviation """tau=C['s3']+(C['s4_noJP']-C['s3'])/2.*(mag-5.)tau[mag<5]=C['s3']tau[mag>7]=C['s4_noJP']returntau@_get_tau_regional.add("JPN")def_get_tau_regional_JPN(region,C,mag):""" Returns regional inter-event standard deviation (Tau) for Japan """tau=C['s3']+(C['s4_all']-C['s3'])/2.*(mag-5.)tau[mag<5]=C['s3']tau[mag>7]=C['s4_all']returntaudef_get_top_of_rupture_depth_term(C,imt,ctx):""" Compute and return top-of-rupture depth term, see section "Deph Scaling Effects". """returnnp.where(ctx.ztor>=20.,C['a15'],C['a15']*ctx.ztor/20.)def_get_vs30star(vs30,imt):""" This computes and returns the tapered Vs30, in Equations 15 and 16. """# compute the limiting v1 value, see Equation 16.t=imt.periodift<=0.50:v1=1500.0elift<3.0:# changed to -0.351 for additional significant figuresv1=np.exp(-0.351*np.log(t/0.5)+np.log(1500.))else:v1=800.0# set the vs30 star value, see Equation 15.vs30_star=np.ones_like(vs30)*vs30vs30_star[vs30>=v1]=v1returnvs30_star
[docs]classGulerceEtAl2017(GMPE):""" Implements the GKAS16 GMPE by Gulerce et al. (2017) for vertical-component ground motions from the PEER NGA-West2 Project. This model follows the same functional form as in ASK14 by Abrahamson et al. (2014) with minor modifications to the underlying parameters. Note that this is a more updated version than the GMPE described in the original PEER Report 2013/24. Reference: Gulerce, Z., Kamai, R., Abrahamson, N., & Silva, W. (2017) "Ground Motion Prediction Equations for the Vertical Ground Motion Component Based on the NGA-W2 Database", Earthquake Spectra, 33(2), 499-528. """region="CAL"#: Supported tectonic region type is active shallow crust, as part of the#: NGA-West2 Database; re-defined here for clarity.DEFINED_FOR_TECTONIC_REGION_TYPE=const.TRT.ACTIVE_SHALLOW_CRUST#: Supported intensity measure type is spectral acceleration#: at T=0.01 to 10.0 s; see Tables 1a and 1b.DEFINED_FOR_INTENSITY_MEASURE_TYPES={SA}#: Supported intensity measure component is the#: :attr:`~openquake.hazardlib.const.IMC.Vertical` direction component;#: see title.DEFINED_FOR_INTENSITY_MEASURE_COMPONENT=const.IMC.VERTICAL#: Supported standard deviation types are inter-event, intra-event#: and total; see the section for "Equations for Standard Deviation".DEFINED_FOR_STANDARD_DEVIATION_TYPES={const.StdDev.TOTAL,const.StdDev.INTER_EVENT,const.StdDev.INTRA_EVENT}#: Required site parameter is Vs30 only. Unlike in ASK14, the nonlinear#: site response and Z1.0 scaling is not incorporated; see the section#: for "Site Amplification Effects".REQUIRES_SITES_PARAMETERS={'vs30'}#: Required rupture parameters are magnitude, rake, dip, ztor, and width;#: see the section for "Functional Form of the Model".REQUIRES_RUPTURE_PARAMETERS={'mag','rake','dip','ztor','width'}#: Required distance measures are Rrup, Rjb, Ry0 and Rx;#: see the section for "Functional Form of the Model".REQUIRES_DISTANCES={'rrup','rjb','rx','ry0'}
[docs]defcompute(self,ctx:np.recarray,imts,mean,sig,tau,phi):""" See :meth:`superclass method for spec of input and result values. <.base.GroundShakingIntensityModel.compute>` """form,imtinenumerate(imts):C=self.COEFFS[imt]# get the mean valuemean[m]=(_get_basic_term(C,ctx)+_get_faulting_style_term(C,ctx)+_get_site_response_term(C,imt,ctx.vs30)+_get_hanging_wall_term(C,ctx)+_get_top_of_rupture_depth_term(C,imt,ctx))mean[m]+=_get_regional_term(self.region,C,imt,ctx.vs30,ctx.rrup)# get standard deviationssig[m],tau[m],phi[m]=_get_stddevs(self.region,C,imt,ctx)
[docs]classGulerceEtAl2017RegTWN(GulerceEtAl2017):""" Implements the GKAS16 GMPE by Gulerce et al. (2017) for vertical-component ground motions from the PEER NGA-West2 Project. Regional corrections for Taiwan """region="TWN"
[docs]classGulerceEtAl2017RegITA(GulerceEtAl2017):""" Implements the GKAS16 GMPE by Gulerce et al. (2017) for vertical-component ground motions from the PEER NGA-West2 Project. Regional corrections for Italy """region="ITA"
[docs]classGulerceEtAl2017RegMID(GulerceEtAl2017):""" Implements the GKAS16 GMPE by Gulerce et al. (2017) for vertical-component ground motions from the PEER NGA-West2 Project. Regional corrections for Middle East """region="MID"
[docs]classGulerceEtAl2017RegCHN(GulerceEtAl2017):""" Implements the GKAS16 GMPE by Gulerce et al. (2017) for vertical-component ground motions from the PEER NGA-West2 Project. Regional corrections for China """region="CHN"
[docs]classGulerceEtAl2017RegJPN(GulerceEtAl2017):""" Implements the GKAS16 GMPE by Gulerce et al. (2017) for vertical-component ground motions from the PEER NGA-West2 Project. Regional corrections for Japan """region="JPN"