Source code for openquake.hazardlib.gsim.zhao_2006
# -*- coding: utf-8 -*-# vim: tabstop=4 shiftwidth=4 softtabstop=4## Copyright (C) 2012-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:`ZhaoEtAl2006Asc`, :class:`ZhaoEtAl2006SInter`,:class:`ZhaoEtAl2006SSlab`, :class:`ZhaoEtAl2006SInterNSHMP2008` and:class:`ZhaoEtAl2006SSlabNSHMP2014`"""importnumpyasnp# standard acceleration of gravity in m/s**2fromscipy.constantsimportgimportcopyfromopenquake.hazardlib.gsim.baseimportGMPE,CoeffsTablefromopenquake.hazardlibimportconstfromopenquake.hazardlib.imtimportPGA,SAfromopenquake.hazardlib.gsim.mgmpe.cb14_basin_termimport_get_cb14_basin_termdef_get_us23_nshm_adjustments(ln_mean,imt,ctx,cb14_basin_term,m9_basin_term):""" If specified get the US 2023 NSHM basin adjustments: 1) The ZhaoEtAl2006 GMM lacks a basin term and therefore uses the Campbell and Bozorgnia (2014) (z2pt5-based) basin term. 2) If Seattle Basin region, the M9 basin term is applied instead of the CB14 basin term if long period ground-motion (T >= 1.9s) and is a deep basin site (z2pt5 >= 6.0 km) """# Set a null basin termfb=np.zeros(len(ln_mean))# Apply cb14 basin term if specifiedifcb14_basin_term:fb=_get_cb14_basin_term(imt,ctx)# Apply m9 basin term if specified (will override# cb14 basin term for basin sites if T >= 1.9 s)ifm9_basin_termandimt.period>=1.9:fb[ctx.z2pt5>=6.0]=np.log(2.0)# Basin sites use m9 basinreturnfbdef_compute_distance_term(C,mag,rrup):""" Compute second and third terms in equation 1, p. 901. """term1=C['b']*rrupterm2=-np.log(rrup+C['c']*np.exp(C['d']*mag))returnterm1+term2def_compute_faulting_style_term(C,rake):""" Compute fifth term in equation 1, p. 901. """# p. 900. "The differentiation in focal mechanism was# based on a rake angle criterion, with a rake of +/- 45# as demarcation between dip-slip and strike-slip."return((rake>45.0)&(rake<135.0))*C['FR']def_compute_focal_depth_term(C,hypo_depth):""" Compute fourth term in equation 1, p. 901. """# p. 901. "(i.e, depth is capped at 125 km)".focal_depth=np.clip(hypo_depth,0,125.)# p. 902. "We used the value of 15 km for the# depth coefficient hc ...".hc=15.0# p. 901. "When h is larger than hc, the depth terms takes# effect ...". The next sentence specifies h>=hc.return(focal_depth>=hc)*C['e']*(focal_depth-hc)def_compute_magnitude_squared_term(P,M,Q,W,mag):""" Compute magnitude squared term, equation 5, p. 909. """returnP*(mag-M)+Q*(mag-M)**2+Wdef_compute_magnitude_term(C,mag):""" Compute first term in equation 1, p. 901. """returnC['a']*magdef_compute_site_class_term(C,vs30):""" Compute nine-th term in equation 1, p. 901. """# map vs30 value to site class, see table 2, p. 901.site_term=np.zeros(len(vs30))# hard rocksite_term[vs30>1100.0]=C['CH']# rocksite_term[(vs30>600)&(vs30<=1100)]=C['C1']# hard soilsite_term[(vs30>300)&(vs30<=600)]=C['C2']# medium soilsite_term[(vs30>200)&(vs30<=300)]=C['C3']# soft soilsite_term[vs30<=200]=C['C4']returnsite_termdef_compute_slab_correction_term(C,rrup):""" Compute path modification term for slab events, that is the 8-th term in equation 1, p. 901. """slab_term=C['SSL']*np.log(rrup)returnslab_termdef_set_stddevs(sig,tau,phi,Cs,Ct):""" Set standard deviations as defined in equation 3 p. 902. """sig[:]=np.sqrt(Cs**2+Ct**2)tau[:]=Ctphi[:]=Cs
[docs]classZhaoEtAl2006Asc(GMPE):""" Implements GMPE developed by John X. Zhao et al. and published as "Attenuation Relations of Strong Ground Motion in Japan Using Site Classification Based on Predominant Period" (2006, Bulletin of the Seismological Society of America, Volume 96, No. 3, pages 898-913). This class implements the equations for 'Active Shallow Crust' (that's why the class name ends with 'Asc'). """#: Supported tectonic region type is active shallow crust, this means#: that factors SI, SS and SSL are assumed 0 in equation 1, p. 901.DEFINED_FOR_TECTONIC_REGION_TYPE=const.TRT.ACTIVE_SHALLOW_CRUST#: Supported intensity measure types are spectral acceleration,#: and peak ground acceleration, see paragraph 'Development of Base Model'#: p. 901.DEFINED_FOR_INTENSITY_MEASURE_TYPES={PGA,SA}#: Supported intensity measure component is geometric mean#: of two horizontal components :#: attr:`~openquake.hazardlib.const.IMC.GEOMETRIC_MEAN`, see paragraph#: 'Development of Base Model', p. 901.DEFINED_FOR_INTENSITY_MEASURE_COMPONENT=const.IMC.GEOMETRIC_MEAN#: Supported standard deviation types are inter-event, intra-event#: and total, see equation 3, p. 902.DEFINED_FOR_STANDARD_DEVIATION_TYPES={const.StdDev.TOTAL,const.StdDev.INTER_EVENT,const.StdDev.INTRA_EVENT}#: Required site parameters is Vs30.#: See table 2, p. 901.REQUIRES_SITES_PARAMETERS={'vs30'}#: Required rupture parameters are magnitude, rake, and focal depth.#: See paragraph 'Development of Base Model', p. 901.REQUIRES_RUPTURE_PARAMETERS={'mag','rake','hypo_depth'}#: Required distance measure is Rrup.#: See paragraph 'Development of Base Model', p. 902.REQUIRES_DISTANCES={'rrup'}#: Reference conditions. See Table 2 at page 901. The hard rock conditions#: is 1100 m/s. Here we force it to 800 to make it compatible with a#: generic site term# DEFINED_FOR_REFERENCE_VELOCITY = 1100DEFINED_FOR_REFERENCE_VELOCITY=800def__init__(self,cb14_basin_term=False,m9_basin_term=False):ifcb14_basin_termorm9_basin_term:self.REQUIRES_SITES_PARAMETERS=frozenset(self.REQUIRES_SITES_PARAMETERS|{'z2pt5'})self.cb14_basin_term=cb14_basin_termself.m9_basin_term=m9_basin_term
[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. """form,imtinenumerate(imts):# extracting dictionary of coefficients specific to required# intensity measure type.C=self.COEFFS_ASC[imt]# mean value as given by equation 1, p. 901, without considering# interface and intraslab terms (that is SI, SS, SSL = 0) and# inter and intra event terms, plus the magnitude-squared term# correction factor (equation 5 p. 909).mean_i=_compute_magnitude_term(C,ctx.mag)+\
_compute_distance_term(C,ctx.mag,ctx.rrup)+\
_compute_focal_depth_term(C,ctx.hypo_depth)+\
_compute_faulting_style_term(C,ctx.rake)+\
_compute_site_class_term(C,ctx.vs30)+\
_compute_magnitude_squared_term(P=0.0,M=6.3,Q=C['QC'],W=C['WC'],mag=ctx.mag)# convert from cm/s**2 to g and back into natural logln_mean=np.log(np.exp(mean_i)*1e-2/g)# Get basin adjustments if specifiedfb=_get_us23_nshm_adjustments(ln_mean,imt,ctx,self.cb14_basin_term,self.m9_basin_term)mean[m]=ln_mean+fb_set_stddevs(sig[m],tau[m],phi[m],C['sigma'],C['tauC'])
[docs]classZhaoEtAl2006SInter(ZhaoEtAl2006Asc):""" Implements GMPE developed by John X. Zhao et al and published as "Attenuation Relations of Strong Ground Motion in Japan Using Site Classification Based on Predominant Period" (2006, Bulletin of the Seismological Society of America, Volume 96, No. 3, pages 898-913). This class implements the equations for 'Subduction Interface' (that's why the class name ends with 'SInter'). This class extends the :class:`openquake.hazardlib.gsim.zhao_2006.ZhaoEtAl2006Asc` because the equation for subduction interface is obtained from the equation for active shallow crust, by removing the faulting style term and adding a subduction interface term. """#: Supported tectonic region type is subduction interface, this means#: that factors FR, SS and SSL are assumed 0 in equation 1, p. 901.DEFINED_FOR_TECTONIC_REGION_TYPE=const.TRT.SUBDUCTION_INTERFACE#: Required rupture parameters are magnitude and focal depth.REQUIRES_RUPTURE_PARAMETERS={'mag','hypo_depth'}
[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. """form,imtinenumerate(imts):# extracting dictionary of coefficients specific to required# intensity measure type.C=self.COEFFS_ASC[imt]C_SINTER=self.COEFFS_SINTER[imt]# mean value as given by equation 1, p. 901, without considering# faulting style and intraslab terms (that is FR, SS, SSL = 0) and# inter and intra event terms, plus the magnitude-squared term# correction factor (equation 5 p. 909)mean_i=_compute_magnitude_term(C,ctx.mag)+\
_compute_distance_term(C,ctx.mag,ctx.rrup)+\
_compute_focal_depth_term(C,ctx.hypo_depth)+\
_compute_site_class_term(C,ctx.vs30)+ \
_compute_magnitude_squared_term(P=0.0,M=6.3,Q=C_SINTER['QI'],W=C_SINTER['WI'],mag=ctx.mag)+C_SINTER['SI']# Convert from cm/s**2 to g then back into natural logln_mean=np.log(np.exp(mean_i)*1e-2/g)# Get basin adjustments if specifiedfb=_get_us23_nshm_adjustments(ln_mean,imt,ctx,self.cb14_basin_term,self.m9_basin_term)mean[m]=ln_mean+fb_set_stddevs(sig[m],tau[m],phi[m],C['sigma'],C_SINTER['tauI'])
[docs]classZhaoEtAl2006SSlab(ZhaoEtAl2006Asc):""" Implements GMPE developed by John X. Zhao et al and published as "Attenuation Relations of Strong Ground Motion in Japan Using Site Classification Based on Predominant Period" (2006, Bulletin of the Seismological Society of America, Volume 96, No. 3, pages 898-913). This class implements the equations for 'Subduction Slab'. (that's why the class name ends with 'SSlab'). This class extends the :class:`openquake.hazardlib.gsim.zhao_2006.ZhaoEtAl2006Asc` because the equation for subduction slab is obtained from the equation for active shallow crust, by removing the faulting style term and adding subduction slab terms. """#: Supported tectonic region type is subduction interface, this means#: that factors FR, SS and SSL are assumed 0 in equation 1, p. 901.DEFINED_FOR_TECTONIC_REGION_TYPE=const.TRT.SUBDUCTION_INTRASLAB#: Required rupture parameters are magnitude and focal depth.REQUIRES_RUPTURE_PARAMETERS={'mag','hypo_depth'}
[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. """form,imtinenumerate(imts):# extracting dictionary of coefficients specific to required# intensity measure type.C=self.COEFFS_ASC[imt]C_SSLAB=self.COEFFS_SSLAB[imt]# to avoid singularity at 0.0 (in the calculation of the# slab correction term), replace 0 values with 0.1d=np.array(ctx.rrup)# make a copyd[d==0.0]=0.1# mean value as given by equation 1, p. 901, without considering# faulting style and intraslab terms (that is FR, SS, SSL = 0) and# inter and intra event terms, plus the magnitude-squared term# correction factor (equation 5 p. 909)mean_i=_compute_magnitude_term(C,ctx.mag)+\
_compute_distance_term(C,ctx.mag,d)+\
_compute_focal_depth_term(C,ctx.hypo_depth)+\
_compute_site_class_term(C,ctx.vs30)+\
_compute_magnitude_squared_term(P=C_SSLAB['PS'],M=6.5,Q=C_SSLAB['QS'],W=C_SSLAB['WS'],mag=ctx.mag)+C_SSLAB['SS']+\
_compute_slab_correction_term(C_SSLAB,d)# Convert from cm/s**2 to g and back into natural logln_mean=np.log(np.exp(mean_i)*1e-2/g)# Get basin adjustments if specifiedfb=_get_us23_nshm_adjustments(ln_mean,imt,ctx,self.cb14_basin_term,self.m9_basin_term)mean[m]=ln_mean+fb_set_stddevs(sig[m],tau[m],phi[m],C['sigma'],C_SSLAB['tauS'])
[docs]classZhaoEtAl2006SInterNSHMP2008(ZhaoEtAl2006SInter):""" Extend :class:`ZhaoEtAl2006SInter` and fix hypocentral depth at 20 km as defined the by National Seismic Hazard Mapping Project for the 2008 US hazard model. The calculation of the total standard deviation is done considering the inter-event standard deviation as defined in table 5, page 903 of Zhao's paper. The class implement the equation as coded in ``subroutine zhao`` in ``hazSUBXnga.f`` Fotran code available at: http://earthquake.usgs.gov/hazards/products/conterminous/2008/software/ """
[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. Call super class method with hypocentral depth fixed at 20 km """# create new rupture context to avoid changing the original onectx=copy.copy(ctx)ctx.hypo_depth=20.super().compute(ctx,imts,mean,sig,tau,phi)
[docs]classZhaoEtAl2006SSlabNSHMP2014(ZhaoEtAl2006SSlab):""" For the 2014 US National Seismic Hazard Maps the magnitude of Zhao et al. (2006) for the subduction inslab events is capped at magnitude Mw 7.8 """
[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. """form,imtinenumerate(imts):# extracting dictionary of coefficients specific to required# intensity measure type.C=self.COEFFS_ASC[imt]C_SSLAB=self.COEFFS_SSLAB[imt]# to avoid singularity at 0.0 (in the calculation of the# slab correction term), replace 0 values with 0.1d=np.array(ctx.rrup)# make a copyd[d==0.0]=0.1rup_mag=np.clip(ctx.mag,0.,7.8)# mean value as given by equation 1, p. 901, without considering# faulting style and intraslab terms (that is FR, SS, SSL = 0) and# inter and intra event terms, plus the magnitude-squared term# correction factor (equation 5 p. 909)mean_i=_compute_magnitude_term(C,rup_mag)+\
_compute_distance_term(C,rup_mag,d)+\
_compute_focal_depth_term(C,ctx.hypo_depth)+\
_compute_site_class_term(C,ctx.vs30)+\
_compute_magnitude_squared_term(P=C_SSLAB['PS'],M=6.5,Q=C_SSLAB['QS'],W=C_SSLAB['WS'],mag=rup_mag)+\
C_SSLAB['SS']+_compute_slab_correction_term(C_SSLAB,d)# Convert from cm/s**2 to g and back into natural logln_mean=np.log(np.exp(mean_i)*1e-2/g)# Get basin adjustments if specifiedfb=_get_us23_nshm_adjustments(ln_mean,imt,ctx,self.cb14_basin_term,self.m9_basin_term)mean[m]=ln_mean+fb_set_stddevs(sig[m],tau[m],phi[m],C['sigma'],C_SSLAB['tauS'])
# Coefficient table taken from Gail Atkinson's "White paper on# Proposed Ground-motion Prediction Equations (GMPEs) for 2015# National Seismic Hazard Maps" (2012, page 16).# Values were interpolated to include all listed periods.# MF is the linear multiplicative factor.COEFFS_SITE_FACTORS=CoeffsTable(sa_damping=5,table="""\ IMT MF pga 0.50 pgv 1.00 0.05 0.44 0.10 0.44 0.15 0.53 0.20 0.60 0.25 0.72 0.30 0.81 0.40 1.00 0.50 1.01 0.60 1.02 0.70 1.02 0.80 1.03 0.90 1.04 1.00 1.04 1.25 1.19 1.50 1.31 2.00 1.51 2.50 1.34 3.00 1.21 4.00 1.09 5.00 1.00 """)
[docs]classZhaoEtAl2006SInterCascadia(ZhaoEtAl2006SInter):""" Implements the interface GMPE developed by John X. Zhao et al modified by the Japan/Cascadia site factors as proposed by Atkinson, G. M. (2012). White paper on proposed ground-motion prediction equations (GMPEs) for 2015 National Seismic Hazard Maps Final Version, Nov. 2012, 50 pp. This class extends the :class:`openquake.hazardlib.gsim.zhao_2006.ZhaoEtAl2006Asc` because the equation for subduction interface is obtained from the equation for active shallow crust, by removing the faulting style term and adding a subduction interface term. """
[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. """form,imtinenumerate(imts):# extracting dictionary of coefficients specific to required# intensity measure type.C=self.COEFFS_ASC[imt]C_SINTER=self.COEFFS_SINTER[imt]C_SF=COEFFS_SITE_FACTORS[imt]# mean value as given by equation 1, p. 901, without considering# faulting style and intraslab terms (that is FR, SS, SSL = 0) and# inter and intra event terms, plus the magnitude-squared term# correction factor (equation 5 p. 909)mean_i=_compute_magnitude_term(C,ctx.mag)+\
_compute_distance_term(C,ctx.mag,ctx.rrup)+\
_compute_focal_depth_term(C,ctx.hypo_depth)+\
_compute_site_class_term(C,ctx.vs30)+ \
_compute_magnitude_squared_term(P=0.0,M=6.3,Q=C_SINTER['QI'],W=C_SINTER['WI'],mag=ctx.mag)+C_SINTER['SI']# multiply by site factor to "convert" Japan values to Cascadia# values then convert from cm/s**2 to gln_mean=np.log((np.exp(mean_i)*C_SF["MF"])*1e-2/g)# Get basin adjustments if specifiedfb=_get_us23_nshm_adjustments(ln_mean,imt,ctx,self.cb14_basin_term,self.m9_basin_term)mean[m]=ln_mean+fb_set_stddevs(sig[m],tau[m],phi[m],C['sigma'],C_SINTER['tauI'])
[docs]classZhaoEtAl2006SSlabCascadia(ZhaoEtAl2006SSlab):""" Implements GMPE developed by John X. Zhao et al modified by the Japan/Cascadia site factors as proposed by Atkinson, G. M. (2012). White paper on proposed ground-motion prediction equations (GMPEs) for 2015 National Seismic Hazard Maps Final Version, Nov. 2012, 50 pp. This class extends the :class:`openquake.hazardlib.gsim.zhao_2006.ZhaoEtAl2006Asc` because the equation for subduction slab is obtained from the equation for active shallow crust, by removing the faulting style term and adding subduction slab terms. """
[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. """form,imtinenumerate(imts):# extracting dictionary of coefficients specific to required# intensity measure type.C=self.COEFFS_ASC[imt]C_SSLAB=self.COEFFS_SSLAB[imt]C_SF=COEFFS_SITE_FACTORS[imt]# to avoid singularity at 0.0 (in the calculation of the# slab correction term), replace 0 values with 0.1d=np.array(ctx.rrup)# make a copyd[d==0.0]=0.1# mean value as given by equation 1, p. 901, without considering# faulting style and intraslab terms (that is FR, SS, SSL = 0) and# inter and intra event terms, plus the magnitude-squared term# correction factor (equation 5 p. 909)mean_i=_compute_magnitude_term(C,ctx.mag)+\
_compute_distance_term(C,ctx.mag,d)+\
_compute_focal_depth_term(C,ctx.hypo_depth)+\
_compute_site_class_term(C,ctx.vs30)+\
_compute_magnitude_squared_term(P=C_SSLAB['PS'],M=6.5,Q=C_SSLAB['QS'],W=C_SSLAB['WS'],mag=ctx.mag)+\
C_SSLAB['SS']+_compute_slab_correction_term(C_SSLAB,d)# multiply by site factor to "convert" Japan values to Cascadia# values then convert from cm/s**2 to g and back into natural logln_mean=np.log((np.exp(mean_i)*C_SF["MF"])*1e-2/g)# Get basin adjustments if specifiedfb=_get_us23_nshm_adjustments(ln_mean,imt,ctx,self.cb14_basin_term,self.m9_basin_term)mean[m]=ln_mean+fb_set_stddevs(sig[m],tau[m],phi[m],C['sigma'],C_SSLAB['tauS'])
[docs]classZhaoEtAl2006AscSGS(ZhaoEtAl2006Asc):""" This class extends the original base class :class:`openquake.hazardlib.gsim.zhao_2006.ZhaoEtAl2006Asc` by introducing a distance filter for the near field, as implemented by SGS for the national PSHA model for Saudi Arabia. """
[docs]defcompute(self,ctx:np.recarray,imts,mean,sig,tau,phi):""" Using a minimum distance of 5km for the calculation. """ctx=copy.deepcopy(ctx)ctx.rrup[ctx.rrup<=5.]=5.super().compute(ctx,imts,mean,sig,tau,phi)