Merge pull request #10131 from gem/parker2020_epistemic

Adding epistemic scaling factors for ParkerEtAl2020 GMM

debian/changelog

@@ -1,3 +1,7 @@
+  [Christopher Brooks]
+  * Added epistemic uncertainty scaling capabilities to
+    the Parker et al. (2020) NGASUB GMM.
+
   [Michele Simionato]
   * Added parameter `with_betw_ratio`
   * Added command `oq info peril`

openquake/hazardlib/gsim/parker_2020.py

@@ -22,22 +22,63 @@
                :class:`ParkerEtAl2020SSlab`
                :class:`ParkerEtAl2020SSlabB`
 """
+import os
 import math
-
 import numpy as np
+import pandas as pd
 from scipy.special import erf
 
 from openquake.baselib.general import CallableDict
 from openquake.hazardlib import const
 from openquake.hazardlib.gsim.base import GMPE, CoeffsTable, add_alias
 from openquake.hazardlib.imt import PGA, SA, PGV
 
+EPI_ADJS = os.path.join(os.path.dirname(__file__),
+                        "parker_2020_epi_adj_table.csv")
+
 CONSTANTS = {"b4": 0.1, "f3": 0.05, "Vb": 200,
              "vref_fnl": 760, "V1": 270, "vref": 760}
 
 _a0 = CallableDict()
 
 
+def _get_sigma_mu_adjustment(sat_region, trt, imt, epi_adjs_table):
+    """
+    Get the sigma_mu_adjustment (epistemic uncertainty) factor to be applied
+    to the mean ground-motion. Values are provided by authors in the
+    electronic supplement for PGA and SA (not PGV) for each saturation regions.
+    """
+    # Map region to those within the adjustment table
+    if sat_region is None:
+        e_reg = 'Global'
+    elif sat_region in ['TW_N', 'TW_S']:
+        e_reg = 'Taiwan'
+    else:
+        e_reg = sat_region
+
+    # Get values for the saturation region
+    adjs = epi_adjs_table.loc[e_reg]
+    if trt == const.TRT.SUBDUCTION_INTERFACE:
+        add = 'interface'
+    else:
+        add = 'intraslab'
+
+    # Get period-dependent epistemic standard deviation (equation 27)
+    period = imt.period
+    if period < adjs[f'T1_{add}']:
+        eps_std = adjs[f'SigEp1_{add}']
+    elif period >= adjs[f'T1_{add}'] and period < adjs[f'T2_{add}']:
+        p1 = adjs[f'SigEp1_{add}'
+                  ] - (adjs[f'SigEp1_{add}'] - adjs[f'SigEp2_{add}'])
+        p2 = (np.log(period/adjs[f'T1_{add}']) / 
+              np.log(adjs[f'T2_{add}']/adjs[f'T1_{add}']))
+        eps_std = p1 * p2
+    else:  # Must be SA with a period larger than or equal to T2
+        eps_std = adjs[f'SigEp2_{add}']
+
+    return eps_std
+
+
 @_a0.add(const.TRT.SUBDUCTION_INTERFACE)
 def _a0_1(trt, region, basin, C, C_PGA):
     """
@@ -310,6 +351,17 @@ def get_stddevs(C, rrup, vs30):
 class ParkerEtAl2020SInter(GMPE):
     """
     Implements Parker et al. (2020) for subduction interface.
+
+    :param str region: Choice of sub region ("AK", "CAM", "SA", "TW",
+                       "Cascadia", "JP").
+    :param str saturation_region: Choice of saturation region ("Aleutian",
+                                  "AK", "Cascadia", "CAM_S", "CAM_N", "JP_Pac",
+                                  "JP_Phi", "SA_N", "SA_S", "TW_W", "TW_E")
+    :param str basin: Choice of basin region ("Out" or "Seattle")
+    :param float sigma_mu_epsilon: Number of standard deviations to multiply
+                                   sigma_mu (which is the standard deviations 
+                                   of the median) for the epistemic uncertainty
+                                   model
     """
 
     DEFINED_FOR_TECTONIC_REGION_TYPE = const.TRT.SUBDUCTION_INTERFACE
@@ -333,9 +385,11 @@ class ParkerEtAl2020SInter(GMPE):
     #: Required distance measure is closest distance to rupture, for
     #: interface events
     REQUIRES_DISTANCES = {'rrup'}
-    REQUIRES_ATTRIBUTES = {'region', 'saturation_region', 'basin'}
+    REQUIRES_ATTRIBUTES = {'region', 'saturation_region', 'basin', 
+                           'sigma_mu_epsilon'}
 
-    def __init__(self, region=None, saturation_region=None, basin=None):
+    def __init__(self, region=None, saturation_region=None, basin=None, 
+                 sigma_mu_epsilon=0.0):
         """
         Enable setting regions to prevent messy overriding
         and code duplication.
@@ -346,6 +400,9 @@ def __init__(self, region=None, saturation_region=None, basin=None):
         else:
             self.saturation_region = saturation_region
         self.basin = basin
+        self.sigma_mu_epsilon = sigma_mu_epsilon
+        with open(EPI_ADJS) as f:
+            self.epi_adjs_table = pd.read_csv(f.name).set_index('Region')
 
     def compute(self, ctx: np.recarray, imts, mean, sig, tau, phi):
         """
@@ -381,6 +438,12 @@ def compute(self, ctx: np.recarray, imts, mean, sig, tau, phi):
             # Take the exponential to get PGA, PSA in g or the PGV in cm/s
             mean[m] = fp + fnl + fb + flin + fm + c0 + fd
 
+            if self.sigma_mu_epsilon and imt != PGV: # Assume don't apply to PGV
+                # Apply epistemic uncertainty scaling
+                sigma_mu_adjust = _get_sigma_mu_adjustment(
+                    self.saturation_region, trt, imt, self.epi_adjs_table)
+                mean[m] += sigma_mu_adjust * self.sigma_mu_epsilon
+
             sig[m], tau[m], phi[m] = get_stddevs(C, ctx.rrup, ctx.vs30)
 
     COEFFS = CoeffsTable(sa_damping=5, table="""\
@@ -454,6 +517,7 @@ class ParkerEtAl2020SSlabB(ParkerEtAl2020SSlab):
     REQUIRES_SITES_PARAMETERS = {'vs30', 'z2pt5'}
 
 
+
 add_alias('ParkerEtAl2020SInterAleutian', ParkerEtAl2020SInter,
           region="AK", saturation_region="Aleutian")
 add_alias('ParkerEtAl2020SInterAlaska', ParkerEtAl2020SInter,

openquake/hazardlib/gsim/parker_2020_epi_adj_table.csv

@@ -0,0 +1,12 @@
+Region,SigEp1_interface,SigEp2_interface,T1_interface,T2_interface,SigEp1_intraslab,SigEp2_intraslab,T1_intraslab,T2_intraslab
+Global,0.4,0.4,0.2,0.4,0.35,0.22,0.15,2
+AK,0.15,0.1,1,4,0.15,0.12,0.5,1
+Aleutian,0.15,0.1,1,4,0.25,0.18,0.3,0.8
+Cascadia,0.43,0.33,0.2,0.5,0.35,0.16,0.2,3
+CAM_N,0.2,0.1,0.7,4,0.25,0.19,0.3,0.8
+CAM_S,0.2,0.1,0.7,4,0.25,0.19,0.3,0.8
+JP_Pac,0.13,0.1,0.7,4,0.13,0.11,0.3,0.8
+JP_Phi,0.13,0.1,0.7,4,0.21,0.15,0.3,0.7
+SA_N,0.13,0.1,0.7,4,0.21,0.15,0.3,0.7
+SA_S,0.13,0.1,0.7,4,0.21,0.15,0.3,0.7
+Taiwan,0.23,0.14,0.6,3,0.13,0.11,0.3,0.8