Source code for openquake.hazardlib.gsim.nz22.nz_nshm2022_parker
"""New Zealand National Seismic Hazard Model 2022 Revision modification of ParkerEtAl2020 GMMs.Soil nonlinearity sigma model modification of ParkerEtAl2020.Bradley, B. A., S. Bora, R. L. Lee, E. F. Manea, M. C. Gerstenberger, P.J. Stafford, G. M. Atkinson, G. Weatherill, J. Hutchinson, C. A. dela Torre, et al. (2022). Summary of the ground-motion character-isation model for the 2022 New Zealand National Seismic HazardModel, GNS Science Rept. 2022/46, GNS Science, Lower Hutt, NewZealand, 44 pp.Bradley, B., S. Bora, R. Lee, E. Manea, M. Gerstenberger, P. Stafford, G.Atkinson, G. Weatherill, J. Hutchinson, C. de la Torre, et al.(2023). Summary of the ground-motion characterisation modelfor the 2022 New Zealand National Seismic Hazard Model,Bull. Seismol. Soc. Am.Lee, R.L., B.A. Bradley, E.F. Manea, J.A. Hutchinson, S.S.Bora (2022). Evaluation of Empirical Ground-Motion Models for the 2022 NewZealand National Seismic Hazard Model Revision, Bull. Seismol. Soc. Am.Module exports :class:`NZNSHM2022_ParkerEtAl2020SInter` :class:`NZNSHM2022_ParkerEtAl2020SInterB` :class:`NZNSHM2022_ParkerEtAl2020SSlab` :class:`NZNSHM2022_ParkerEtAl2020SSlabB`"""importmathimportnumpyasnpfromscipy.interpolateimportinterp1dfromopenquake.hazardlibimportconstfromopenquake.hazardlib.imtimportPGAfromopenquake.hazardlib.gsim.baseimportCoeffsTable,add_aliasfromopenquake.hazardlib.gsim.nz22.constimport(periods_AG20,rho_Ws,rho_Bs,periods,theta7s,theta8s,)fromopenquake.hazardlib.gsim.parker_2020import(ParkerEtAl2020SInter,_c0,_depth_scaling,_linear_amplification,_magnitude_scaling,_path_term,_get_basin_term,CONSTANTS,)def_non_linear_term(C,imt,vs30,fp,fm,c0,fd=0):""" Non-linear site term. The hard coded fnl = 0 for T >=3 is removed in the NZ version of the model (personal communication with Grace). """# fd for slab onlypgar=np.exp(fp+fm+c0+fd)fnl=(C["f4"]*(np.exp(C["f5"]*(np.minimum(vs30,CONSTANTS["vref_fnl"])-CONSTANTS["Vb"]))-math.exp(C["f5"]*(CONSTANTS["vref_fnl"]-CONSTANTS["Vb"])))*np.log((pgar+CONSTANTS["f3"])/CONSTANTS["f3"]))returnfnl
[docs]defget_stddevs(C,rrup,vs30):""" Returns the standard deviations. Generate tau, phi, and total sigma computed from both total and partitioned phi models. """# define period-independent coefficients for phi modelsv1=200v2=500r1=200r2=500# total Phiphi_rv=np.zeros(len(vs30))fori,vs30iinenumerate(vs30):ifrrup[i]<=r1:phi_rv[i]=C["phi21"]elifrrup[i]>=r2:phi_rv[i]=C["phi22"]else:phi_rv[i]=((C["phi22"]-C["phi21"])/(math.log(r2)-math.log(r1)))*(math.log(rrup[i])-math.log(r1))+C["phi21"]ifvs30i<=v1:phi_rv[i]+=C["phi2V"]*(math.log(r2/max(r1,min(r2,rrup[i])))/math.log(r2/r1))elifvs30i<v2:phi_rv[i]+=(C["phi2V"]*((math.log(v2/min(v2,vs30i)))/math.log(v2/v1))*(math.log(r2/max(r1,min(r2,rrup[i])))/math.log(r2/r1)))phi_tot=np.sqrt(phi_rv)return[np.sqrt(C["Tau"]**2+phi_tot**2),C["Tau"],phi_tot]
[docs]defget_nonlinear_stddevs(C,C_PGA,imt,pgar,rrup,vs30):""" This NZ specific modification. Get the nonlinear tau and phi terms for Parker's model. This routine is based upon Peter Stafford suggested implementation shared on slack, which is based on AG20 implementation. """period=imt.period# Linear Tautau_lin=C["Tau"]*np.ones(vs30.shape)tau_lin_pga=C_PGA["Tau"]*np.ones(vs30.shape)r1=200.0r2=500.0# Linear Phiphi2_rv=np.zeros(len(vs30))phi2_rv_pga=np.zeros(len(vs30))fori,vs30iinenumerate(vs30):ifrrup[i]<=r1:phi2_rv[i]=C["phi21"]phi2_rv_pga[i]=C_PGA["phi21"]elifrrup[i]>=r2:phi2_rv[i]=C["phi22"]phi2_rv_pga[i]=C_PGA["phi22"]else:phi2_rv[i]=((C["phi22"]-C["phi21"])/(math.log(r2)-math.log(r1)))*(math.log(rrup[i])-math.log(r1))+C["phi21"]phi2_rv_pga[i]=((C_PGA["phi22"]-C_PGA["phi21"])/(math.log(r2)-math.log(r1)))*(math.log(rrup[i])-math.log(r1))+C_PGA["phi21"]phi_lin=np.sqrt(phi2_rv)phi_lin_pga=np.sqrt(phi2_rv_pga)# Assume that the site response variability is constant with period.phi_amp=0.3phi_B=np.sqrt(phi_lin**2-phi_amp**2)phi_B_pga=np.sqrt(phi_lin_pga**2-phi_amp**2)rho_W_itp=interp1d(np.log(periods_AG20),rho_Ws)rho_B_itp=interp1d(np.log(periods_AG20),rho_Bs)ifperiod<0.01:rhoW=1.0rhoB=1.0else:rhoW=rho_W_itp(np.log(period))rhoB=rho_B_itp(np.log(period))f2=C["f4"]*(np.exp(C["f5"]*(np.minimum(vs30,CONSTANTS["vref_fnl"])-CONSTANTS["Vb"]))-math.exp(C["f5"]*(CONSTANTS["vref_fnl"]-CONSTANTS["Vb"])))f3=CONSTANTS["f3"]partial_f_pga=f2*pgar/(pgar+f3)partial_f_pga=partial_f_pga*np.ones(vs30.shape)# nonlinear variance componentsphi2_NL=(phi_lin**2+partial_f_pga**2*phi_B_pga**2+2*partial_f_pga*phi_B_pga*phi_B*rhoW)tau2_NL=(tau_lin**2+partial_f_pga**2*tau_lin_pga**2+2*partial_f_pga*tau_lin_pga*tau_lin*rhoB)return[np.sqrt(tau2_NL+phi2_NL),np.sqrt(tau2_NL),np.sqrt(phi2_NL)]
[docs]defget_sigma_epistemic(trt,region,imt):""" This is a NZ-NSHM-2022 specific modification. Currently the epistemic sigma model is applied to Global model only. As for NZ we are using only the global model. Henec below the coefficients are just for the global model. """ifregionisNone:iftrt==const.TRT.SUBDUCTION_INTRASLAB:sigma_epsilon1=0.35sigma_epsilon2=0.22T1=0.15T2=2else:sigma_epsilon1=0.4sigma_epsilon2=0.4T1=0.2T2=0.4period=imt.periodifperiod<0.01:sigma_epi=sigma_epsilon1elif(period>=0.01)&(period<T1):sigma_epi=sigma_epsilon1elif(period>=T1)&(period<T2):per_ratio1=np.log(period/T1)per_ratio2=np.log(T2/T1)sigma_epi=sigma_epsilon1-(sigma_epsilon1-sigma_epsilon2)*(per_ratio1/per_ratio2)else:sigma_epi=sigma_epsilon2returnsigma_epielse:return0.0
[docs]defget_backarc_term(trt,imt,ctx):""" This is a NZ NSHM-2022 specific modification. The backarc correction factors to be applied with the ground motion prediction. In the NZ context, it is applied to only subduction intraslab events. It is essentially the correction factor taken from BC Hydro 2016. Abrahamson et al. (2016) Earthquake Spectra. The correction is applied only for sites in the backarc region as function of distance. """period=imt.periodw_epi_factor=1.008theta7_itp=interp1d(np.log(periods[1:]),theta7s[1:])theta8_itp=interp1d(np.log(periods[1:]),theta8s[1:])# Note that there is no correction for PGV.# Hence, I make theta7 and theta8 as 0 for periods < 0.ifperiod<0:theta7=0.0theta8=0.0elifperiod>=0andperiod<0.02:theta7=1.0988theta8=-1.42else:theta7=theta7_itp(np.log(period))theta8=theta8_itp(np.log(period))dists=ctx.rrupiftrt==const.TRT.SUBDUCTION_INTRASLAB:min_dist=85.0backarc=np.bool_(ctx.backarc)f_faba=np.zeros_like(dists)fixed_dists=dists[backarc]fixed_dists[fixed_dists<min_dist]=min_distf_faba[backarc]=theta7+theta8*np.log(fixed_dists/40.0)returnf_faba*w_epi_factorelse:f_faba=np.zeros_like(dists)returnf_faba
[docs]classNZNSHM2022_ParkerEtAl2020SInter(ParkerEtAl2020SInter):""" Implements NZ NSHM 2022 Soil nonlinearity sigma model modification of ParkerEtAl2020SInter for NZ NSHM 2022. """def__init__(self,region=None,saturation_region=None,basin=None,sigma_mu_epsilon=0.0,modified_sigma=False,):""" Enable setting regions to prevent messy overriding and code duplication. """super().__init__(region=region,saturation_region=saturation_region,basin=basin,sigma_mu_epsilon=sigma_mu_epsilon)self.modified_sigma=modified_sigma
[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. """trt=self.DEFINED_FOR_TECTONIC_REGION_TYPEC_PGA=self.COEFFS[PGA()]form,imtinenumerate(imts):C=self.COEFFS[imt]# Regional Mb factorifself.saturation_regioninself.MB_REGIONS:m_b=self.MB_REGIONS[self.saturation_region]else:m_b=self.MB_REGIONS["default"]c0,c0_pga=_c0(trt,self.region,self.saturation_region,C,C_PGA)fm,fm_pga=_magnitude_scaling(self.SUFFIX,C,C_PGA,ctx.mag,m_b)fp,fp_pga=_path_term(trt,self.region,self.basin,self.SUFFIX,C,C_PGA,ctx.mag,ctx.rrup,m_b,)fp_pga=(fp_pga+get_backarc_term(trt,PGA(),ctx))# backarc term applied to path function for reference rock PGAfd=_depth_scaling(trt,C,ctx)fd_pga=_depth_scaling(trt,C_PGA,ctx)fb=_get_basin_term(self.region,self.basin,C,ctx)flin=_linear_amplification(self.region,C,ctx.vs30)fnl=_non_linear_term(C,imt,ctx.vs30,fp_pga,fm_pga,c0_pga,fd_pga)fba=get_backarc_term(trt,imt,ctx)# backarc correction factor from BC Hydro at individual period# The output is the desired median model prediction in LN units# Take the exponential to get PGA, PSA in g or the PGV in cm/smean[m]=fp+fnl+fb+flin+fm+c0+fd+fbaifself.sigma_mu_epsilon:# Apply an epistmic adjustment factor.# Currently, its applied to only global model.mean[m]+=self.sigma_mu_epsilon*get_sigma_epistemic(trt,self.region,imt)# The default sigma is modified sigma that accounts for soil# nonlinearityifself.modified_sigma:pgar=np.exp(fp_pga+fm_pga+c0_pga+fd_pga)# the backarc correction is already applied in f_pgasig[m],tau[m],phi[m]=get_nonlinear_stddevs(C,C_PGA,imt,pgar,ctx.rrup,ctx.vs30)else:sig[m],tau[m],phi[m]=get_stddevs(C,ctx.rrup,ctx.vs30)
[docs]classNZNSHM2022_ParkerEtAl2020SInterB(NZNSHM2022_ParkerEtAl2020SInter):""" For Cascadia and Japan where basins are defined (also require z2pt5). """REQUIRES_SITES_PARAMETERS={"vs30","z2pt5"}
[docs]classNZNSHM2022_ParkerEtAl2020SSlab(NZNSHM2022_ParkerEtAl2020SInter):""" Modifications for subduction slab. """DEFINED_FOR_TECTONIC_REGION_TYPE=const.TRT.SUBDUCTION_INTRASLAB# slab also requires hypo_depthREQUIRES_RUPTURE_PARAMETERS={"mag","hypo_depth"}# slab also requires backarcREQUIRES_SITES_PARAMETERS={"vs30","backarc"}# constant table suffixSUFFIX="slab"MB_REGIONS={"Aleutian":7.98,"AK":7.2,"Cascadia":7.2,"CAM_S":7.6,"CAM_N":7.4,"JP_Pac":7.65,"JP_Phi":7.55,"SA_N":7.3,"SA_S":7.25,"TW_W":7.7,"TW_E":7.7,"default":7.6,}
[docs]classNZNSHM2022_ParkerEtAl2020SSlabB(NZNSHM2022_ParkerEtAl2020SSlab):""" For Cascadia and Japan where basins are defined (also require z2pt5). """REQUIRES_SITES_PARAMETERS={"vs30","z2pt5"}
