Merge pull request #8872 from thpap/swiss-gmpes

Adjustements to Swiss GMPEs for application in the Earthquake Risk Model of Switzerland

debian/changelog

@@ -1,3 +1,16 @@
+  [Athanasios Papadopoulos]
+  * Adjusted the Swiss-specific implementations of the GMPEs used
+    in the Swiss national seismic hazard (SUIhaz15) and risk (ERM-CH23) models.
+    The new implementation returns all of total, intra- and inter-event sigmas,
+    rather than just the total one.  
+  * Extended the ModifiableGMPE class by adding an 
+    apply_swiss_amplification_sa method that is used in ERM-CH23 to
+    apply site specific adjustments for site effects and intra-event
+    uncertainty. 
+  * Added ch_ampl03, ch_ampl06, ch_phis2s03, ch_phis2s06,
+    ch_phisss03, ch_phis2s06 site parameters to the site.py file.
+    These parameters are required by the apply_swiss_amplification_sa method.
+
   [Paolo Tormene]
   * Upgraded requirements: Shapely to version 2.0.1 and Pandas to version 2.0.3
 

openquake/hazardlib/contexts.py

@@ -405,6 +405,13 @@ def __init__(self, trt, gsims, oq, monitor=Monitor(), extraparams=()):
                     reqset.update(gsim.gmpe.REQUIRES_SITES_PARAMETERS)
                     if 'apply_swiss_amplification' in gsim.params:
                         reqset.add('amplfactor')
+                    if ('apply_swiss_amplification_sa' in gsim.params):
+                        reqset.add('ch_ampl03')
+                        reqset.add('ch_ampl06')
+                        reqset.add('ch_phis2s03')
+                        reqset.add('ch_phis2s06')
+                        reqset.add('ch_phiss03')
+                        reqset.add('ch_phiss06')
             setattr(self, 'REQUIRES_' + req, reqset)
         try:
             self.min_iml = param['min_iml']

openquake/hazardlib/gsim/akkar_bommer_2010.py

@@ -287,7 +287,8 @@ class AkkarBommer2010SWISS01(AkkarBommer2010):
 
     Model implemented by laurentiu.danciu@gmail.com
     """
-    DEFINED_FOR_STANDARD_DEVIATION_TYPES = {const.StdDev.TOTAL}
+    DEFINED_FOR_STANDARD_DEVIATION_TYPES = {
+        const.StdDev.TOTAL, const.StdDev.INTER_EVENT, const.StdDev.INTRA_EVENT}
 
     #: Vs30 value representing typical rock conditions in Switzerland.
     #: confirmed by the Swiss GMPE group
@@ -309,8 +310,8 @@ def compute(self, ctx: np.recarray, imts, mean, sig, tau, phi):
                 AkkarBommer2010.COEFFS, self.COEFFS_FS_ROCK[imt], tau_ss,
                 mean[m], sig[m], tau[m], phi[m], ctx, ctx.rjb, imt, log_phi_ss)
             sig[m] = np.log(10 ** sig[m])
-            #tau[m] = np.log(10 ** tau[m])
-            #phi[m] = np.log(10 ** phi[m])
+            tau[m] = np.log(10 ** tau[m])
+            phi[m] = np.log(10 ** phi[m])
 
     COEFFS_FS_ROCK = COEFFS_FS_ROCK_SWISS01
 
@@ -320,7 +321,8 @@ class AkkarBommer2010SWISS04(AkkarBommer2010SWISS01):
     This class extends :class:`AkkarBommer2010` following same strategy
     as for :class:`AkkarBommer2010SWISS01`
     """
-    DEFINED_FOR_STANDARD_DEVIATION_TYPES = {const.StdDev.TOTAL}
+    DEFINED_FOR_STANDARD_DEVIATION_TYPES = {
+        const.StdDev.TOTAL, const.StdDev.INTER_EVENT, const.StdDev.INTRA_EVENT}   
 
     COEFFS_FS_ROCK = COEFFS_FS_ROCK_SWISS04
 
@@ -331,6 +333,7 @@ class AkkarBommer2010SWISS08(AkkarBommer2010SWISS01):
     as for :class:`AkkarBommer2010SWISS01` to be used for the
     Swiss Hazard Model [2014].
     """
-    DEFINED_FOR_STANDARD_DEVIATION_TYPES = {const.StdDev.TOTAL}
+    DEFINED_FOR_STANDARD_DEVIATION_TYPES = {
+        const.StdDev.TOTAL, const.StdDev.INTER_EVENT, const.StdDev.INTRA_EVENT} 
 
     COEFFS_FS_ROCK = COEFFS_FS_ROCK_SWISS08

openquake/hazardlib/gsim/akkar_bommer_2010_swiss_coeffs.py

@@ -34,6 +34,7 @@
     0.3000    0.7376    1.654794E+00    1.255398E+00    9.970978E-01     -2.863545E-01  6.506672E+02   0.64           0.5            0.37           5      7      11      34   0.48
     0.3500    0.7447    1.855761E+00    1.322266E+00    9.964436E-01     -3.280135E-01  4.672187E+02   0.627929292    0.493964646    0.37           5      7      11      34   0.473964646
     0.4000    0.7507    2.030328E+00    1.383505E+00    9.958719E-01     -3.642816E-01  3.078242E+02   0.617473168    0.488736584    0.37           5      7      11      34   0.468736584
+    0.6000    0.7712    2.030328E+00    1.383505E+00    9.958719E-01     -3.642816E-01  3.078242E+02   0.584217936    0.472108968    0.377891032    5      7      11      34   0.459368292
     1.0000    0.7916    -5.169560E+00   1.000000E+00    1.010650E+00     6.221898E-01   1.000000E+09   0.54           0.45           0.4            5      7      11      34   0.45
     1.0500    0.7908    -6.821261E+00   1.000000E+00    1.016859E+00     8.337131E-01   1.000000E+09   0.539555893    0.447779464    0.4            5      7      11      34   0.448223571
     1.1000    0.7912    -8.396109E+00   1.000000E+00    1.022780E+00     1.035395E+00   1.000000E+09   0.539132449    0.445662247    0.4            5      7      11      34   0.446529797
@@ -93,6 +94,7 @@
     0.3000    0.8829    1.654794E+00    1.255398E+00    9.970978E-01    -2.863545E-01   6.506672E+02    0.64           0.5            0.37           5      7      11      34      0.48
     0.3500    0.8696    1.855761E+00    1.322266E+00    9.964436E-01    -3.280135E-01   4.672187E+02    0.627929292    0.493964646    0.37           5      7      11      34      0.473964646
     0.4000    0.8603    2.030328E+00    1.383505E+00    9.958719E-01    -3.642816E-01   3.078242E+02    0.617473168    0.488736584    0.37           5      7      11      34      0.468736584
+    0.6000    0.8461    2.030328E+00    1.383505E+00    9.958719E-01     -3.642816E-01  3.078242E+02    0.584217936    0.472108968    0.377891032    5      7      11      34      0.459368292
     1.0000    0.8376    -5.169560E+00   1.000000E+00    1.010650E+00    6.221898E-01    1.000000E+09    0.54           0.45           0.4            5      7      11      34      0.45
     1.0500    0.8345    -6.821261E+00   1.000000E+00    1.016859E+00    8.337131E-01    1.000000E+09    0.539555893    0.447779464    0.4            5      7      11      34      0.448223571
     1.1000    0.8329    -8.396109E+00   1.000000E+00    1.022780E+00    1.035395E+00    1.000000E+09    0.539132449    0.445662247    0.4            5      7      11      34      0.446529797
@@ -152,6 +154,7 @@
     0.3000    0.9039     1.654794E+00    1.255398E+00    9.970978E-01    -2.863545E-01   6.506672E+02   0.64           0.5            0.37           5      7      11      34    0.48
     0.3500    0.8790     1.855761E+00    1.322266E+00    9.964436E-01    -3.280135E-01   4.672187E+02   0.627929292    0.493964646    0.37           5      7      11      34    0.473964646
     0.4000    0.8608     2.030328E+00    1.383505E+00    9.958719E-01    -3.642816E-01   3.078242E+02   0.617473168    0.488736584    0.37           5      7      11      34    0.468736584
+    0.6000    0.8292     2.030328E+00    1.383505E+00    9.958719E-01    -3.642816E-01   3.078242E+02   0.584217936    0.472108968    0.377891032    5      7      11      34    0.459368292
     1.0000    0.8312     -5.169560E+00   1.000000E+00    1.010650E+00    6.221898E-01    1.000000E+09   0.54           0.45           0.4            5      7      11      34    0.45
     1.0500    0.8317     -6.821261E+00   1.000000E+00    1.016859E+00    8.337131E-01    1.000000E+09   0.539555893    0.447779464    0.4            5      7      11      34    0.448223571
     1.1000    0.8341     -8.396109E+00   1.000000E+00    1.022780E+00    1.035395E+00    1.000000E+09   0.539132449    0.445662247    0.4            5      7      11      34    0.446529797

openquake/hazardlib/gsim/cauzzi_faccioli_2008_swiss.py

@@ -57,8 +57,9 @@ class CauzziFaccioli2008SWISS01(CauzziFaccioli2008):
 
     Model implemented by laurentiu.danciu@gmail.com
     """
-    #: Supported standard deviation type is only total
-    DEFINED_FOR_STANDARD_DEVIATION_TYPES = {const.StdDev.TOTAL}
+    #: Supported standard deviation type is total, inter-event and intra-event
+    DEFINED_FOR_STANDARD_DEVIATION_TYPES = {
+        const.StdDev.TOTAL, const.StdDev.INTER_EVENT, const.StdDev.INTRA_EVENT}
 
     #: Supported intensity measure types are spectral acceleration, peak
     #: ground acceleration and peak ground velocity.
@@ -89,6 +90,8 @@ def compute(self, ctx: np.recarray, imts, mean, sig, tau, phi):
                 mean[m], sig[m], tau[m], phi[m], ctx, ctx.rhypo, imt,
                 log_phi_ss)
             sig[m] = np.log(10 ** sig[m])
+            tau[m] = np.log(10 ** tau[m])
+            phi[m] = np.log(10 ** phi[m])
 
     COEFFS_FS_ROCK = COEFFS_FS_ROCK_SWISS01
 

openquake/hazardlib/gsim/chiou_youngs_2008_swiss.py

@@ -74,7 +74,9 @@ class ChiouYoungs2008SWISS01(ChiouYoungs2008):
 
     Model implemented by laurentiu.danciu@gmail.com
     """
-    DEFINED_FOR_STANDARD_DEVIATION_TYPES = {const.StdDev.TOTAL}
+    #: Supported standard deviation type is total, inter-event and intra-event
+    DEFINED_FOR_STANDARD_DEVIATION_TYPES = {
+        const.StdDev.TOTAL, const.StdDev.INTER_EVENT, const.StdDev.INTRA_EVENT}
 
     #: Vs30 value representing typical rock conditions in Switzerland.
     #: confirmed by the Swiss GMPE group

openquake/hazardlib/gsim/edwards_fah_2013a.py

@@ -157,9 +157,9 @@ class EdwardsFah2013Alpine10Bars(GMPE):
     #: :attr:`~openquake.hazardlib.const.IMC.GEOMETRIC_MEAN`
     DEFINED_FOR_INTENSITY_MEASURE_COMPONENT = const.IMC.GEOMETRIC_MEAN
 
-    #: Supported standard deviation type is total,
-    #: Carlo Cauzzi - Personal Communication
-    DEFINED_FOR_STANDARD_DEVIATION_TYPES = {const.StdDev.TOTAL}
+    #: Supported standard deviation type is total, inter-event and intra-event
+    DEFINED_FOR_STANDARD_DEVIATION_TYPES = {
+        const.StdDev.TOTAL, const.StdDev.INTER_EVENT, const.StdDev.INTRA_EVENT}
 
     #: Required site parameter is only Vs30 (used to distinguish rock
     #: and deep soil).

openquake/hazardlib/gsim/mgmpe/modifiable_gmpe.py

@@ -85,6 +85,27 @@ def apply_swiss_amplification(ctx, imt, me, si, ta, ph):
     me[:] += ctx.amplfactor
 
 
+def apply_swiss_amplification_sa(ctx, imt, me, si, ta, ph):
+    """
+    Adjust Swiss GMPEs to add amplification and correct intra-event residuals
+    """
+
+    if imt.period == 0.3:
+        phis2s = ctx.ch_phis2s03
+        phiss = ctx.ch_phiss03
+        me[:] += ctx.ch_ampl03
+    elif imt.period == 0.6:
+        phis2s = ctx.ch_phis2s06
+        phiss = ctx.ch_phiss06
+        me[:] += ctx.ch_ampl06
+
+    phi_star = np.sqrt(phis2s**2 + phiss**2)
+    total_stddev_star = np.sqrt(ta**2 + phi_star**2)
+
+    ph[:] = phi_star
+    si[:] = total_stddev_star
+
+
 def set_between_epsilon(ctx, imt, me, si, ta, ph, epsilon_tau):
     """
     :param epsilon_tau:
@@ -195,7 +216,7 @@ class ModifiableGMPE(GMPE):
 
     def __init__(self, **kwargs):
         super().__init__(**kwargs)
- 
+
         # Create the original GMPE
         [(gmpe_name, kw)] = kwargs.pop('gmpe').items()
         self.params = kwargs  # non-gmpe parameters

openquake/hazardlib/gsim/utils_swiss_gmpe.py

@@ -82,9 +82,11 @@ def _corr_sig(sig, C, tau_ss, phi_ss, NL=None, tau_value=None):
     s = phi_ss ** 2
     if tau_value is not None and NL is not None:
         s += tau_value * tau_value * ((1 + NL) ** 2)
+        t = tau_value * tau_value * ((1 + NL) ** 2)
     else:
         s += C[tau_ss] * C[tau_ss]
-    return np.sqrt(s)
+        t = C[tau_ss] * C[tau_ss]
+    return np.sqrt(s), np.sqrt(t), phi_ss
 
 
 def _apply_adjustments(COEFFS, C_ADJ, tau_ss, mean, sig, tau, phi, ctx,
@@ -103,4 +105,5 @@ def _apply_adjustments(COEFFS, C_ADJ, tau_ss, mean, sig, tau, phi, ctx,
         _compute_small_mag_correction_term(C_ADJ, ctx.mag, dist)
 
     mean[:] = np.log(mean_corr)
-    sig[:] = _corr_sig(sig, COEFFS[imt], tau_ss, phi_ss, NL, tau_value)
+    sig[:], tau[:], phi[:] = _corr_sig(sig, COEFFS[imt], tau_ss,
+                                       phi_ss, NL, tau_value)

openquake/hazardlib/gsim/zhao_2006_swiss.py

@@ -58,9 +58,11 @@ class ZhaoEtAl2006AscSWISS05(ZhaoEtAl2006Asc):
     Model implemented by laurentiu.danciu@gmail.com
     """
 
-    # Supported standard deviation type is only total, but reported as a
-    # combination of mean and magnitude/distance single station sigma
-    DEFINED_FOR_STANDARD_DEVIATION_TYPES = {const.StdDev.TOTAL}
+    # Supported standard deviation type originally only total, but reported as a
+    # combination of mean and magnitude/distance single station sigma.
+    # updated to return inter and intra event component of st.dev (June 2023)                                                                         
+    DEFINED_FOR_STANDARD_DEVIATION_TYPES = {
+        const.StdDev.TOTAL, const.StdDev.INTER_EVENT, const.StdDev.INTRA_EVENT}
 
     DEFINED_FOR_INTENSITY_MEASURE_TYPES = {PGA, SA}
 

openquake/hazardlib/site.py

@@ -179,6 +179,12 @@ def _extract(array_or_float, indices):
     'h800': numpy.float64,
     'geology': (numpy.string_, 20),
     'amplfactor': numpy.float64,
+    'ch_ampl03': numpy.float64,
+    'ch_ampl06': numpy.float64,
+    'ch_phis2s03': numpy.float64,
+    'ch_phis2s06': numpy.float64,
+    'ch_phiss03': numpy.float64,
+    'ch_phiss06': numpy.float64,
     'fpeak': numpy.float64,
     # Fundamental period and and amplitude of HVRSR spectra
     'THV': numpy.float64,
@@ -391,7 +397,7 @@ def set_global_params(
             self._set('z2pt5', oq.reference_depth_to_2pt5km_per_sec)
         if 'backarc' in req_site_params:
             self._set('backarc', oq.reference_backarc)
-        
+
     def filtered(self, indices):
         """
         :param indices:
@@ -631,8 +637,8 @@ def assoc(self, site_model, assoc_dist, ignore=()):
                 continue
             dt = site_param_dt[param]
             if dt is numpy.float64 and (self.array[param] == 0.).all():
-                raise ValueError('The site parameter %s is always zero: please '
-                                 'check the site model' % param)
+                raise ValueError('The site parameter %s is always zero: please'
+                                 ' check the site model' % param)
         return site_model
 
     def extend(self, lons, lats):

openquake/hazardlib/tests/gsim/mgmpe/modifiable_gmpe_test.py

@@ -25,6 +25,10 @@
 from openquake.hazardlib.tests.gsim.mgmpe.dummy import Dummy
 from openquake.hazardlib.gsim.mgmpe.modifiable_gmpe import (
     ModifiableGMPE, _dict_to_coeffs_table)
+from openquake.hazardlib.gsim.boore_atkinson_2008 import BooreAtkinson2008
+from openquake.hazardlib.gsim.boore_2014 import BooreEtAl2014
+from openquake.hazardlib.imt import from_string
+
 
 
 class ModifiableGMPEAlAtik2015SigmaTest(unittest.TestCase):
@@ -289,7 +293,10 @@ def setUp(self):
         ctx.hypo_depth = 10.
         ctx.occurrence_rate = .001
         sites = Dummy.get_site_collection(
-            4, amplfactor=[-1.0, 1.5, 0.00, -1.99])
+            4, amplfactor=[-1.0, 1.5, 0.00, -1.99],
+            ch_ampl03=[-0.2, 0.4, 0.6, 0],
+            ch_phis2s03=[0.3, 0.4, 0.5, 0.2],
+            ch_phiss03=[0.2, 0.1, 0.3, 0.4])
         for name in sites.array.dtype.names:
             setattr(ctx, name, sites[name])
         ctx.rhypo = ctx.rrup = ctx.repi = np.array([1., 10., 30., 70.])
@@ -315,3 +322,28 @@ def test_get_mean_std(self):
 
             # Check the computed mean + amplification
             np.testing.assert_almost_equal(mean, exp_mean)
+
+        for gmpe_name in ['EdwardsFah2013Alpine10Bars',
+                          'EdwardsFah2013Foreland60Bars',
+                          'ChiouYoungs2008SWISS01']:
+
+            self.imt = from_string('SA(0.3)')
+            gmm = ModifiableGMPE(gmpe={gmpe_name: {}},
+                                 apply_swiss_amplification_sa={})
+            mean, intra_stdev = gmm.get_mean_and_stddevs(
+                                 self.ctx, self.ctx, self.ctx,
+                                 self.imt, [const.StdDev.INTRA_EVENT])
+
+            gmpe = registry[gmpe_name]()
+            emean = get_mean_stds(gmpe, self.ctx, [self.imt],
+                                  truncation_level=0)[0]
+
+            exp_mean = emean + np.array([-0.2, 0.4, 0.6, 0])
+            exp_stdev = np.sqrt(np.array([0.3, 0.4, 0.5, 0.2])**2 +
+                                np.array([0.2, 0.1, 0.3, 0.4])**2)
+
+            # Check the computed mean + amplification
+            np.testing.assert_almost_equal(mean, exp_mean[0])
+
+            # Check the computed intra-event stdev
+            np.testing.assert_almost_equal(intra_stdev[0], exp_stdev)