Adding M9 Seattle Basin Adjustment + Campbell and Bozorgnia (2014) basin term to ModifiableGMPE capabilities

debian/changelog

@@ -1,3 +1,9 @@
+  [Christopher Brooks]
+  * Added M9 Basin amplification term for ModifiableGMPE
+    and assoc. unit tests
+  * Added the CB14 basin amplification term for
+    ModifiableGMPE and assoc. unit tests
+
   [Michele Simionato]
   * Added parameter `with_betw_ratio`
   * Added command `oq info peril`

openquake/hazardlib/contexts.py

@@ -1686,7 +1686,7 @@ def __eq__(self, other):
         return False
 
 
-# mock of a site collection used in the tests and in the SMTK
+# mock of a site collection used in the tests and in the SMT
 class SitesContext(BaseContext):
     """
     Sites calculation context for ground shaking intensity models.
@@ -1699,7 +1699,7 @@ class SitesContext(BaseContext):
     Only those required parameters are made available in a result context
     object.
     """
-    # _slots_ is used in hazardlib check_gsim and in the SMTK
+    # _slots_ is used in hazardlib check_gsim and in the SMT
     def __init__(self, slots='vs30 vs30measured z1pt0 z2pt5'.split(),
                  sitecol=None):
         self._slots_ = slots
@@ -1708,7 +1708,7 @@ def __init__(self, slots='vs30 vs30measured z1pt0 z2pt5'.split(),
             for slot in slots:
                 setattr(self, slot, getattr(sitecol, slot))
 
-    # used in the SMTK
+    # used in the SMT
     def __len__(self):
         return len(self.sids)
 
@@ -1745,7 +1745,7 @@ def get_dists(ctx):
 
 
 # used to produce a RuptureContext suitable for legacy code, i.e. for calls
-# to .get_mean_and_stddevs, like for instance in the SMTK
+# to .get_mean_and_stddevs, like for instance in the SMT
 def full_context(sites, rup, dctx=None):
     """
     :returns: a full RuptureContext with all the relevant attributes
@@ -1785,7 +1785,7 @@ def get_mean_stds(gsim, ctx, imts, **kw):
     return out[:, 0] if single else out
 
 
-# mock of a rupture used in the tests and in the SMTK
+# mock of a rupture used in the tests and in the SMT
 class RuptureContext(BaseContext):
     """
     Rupture calculation context for ground shaking intensity models.

openquake/hazardlib/gsim/mgmpe/cb14_basin_term.py

@@ -0,0 +1,85 @@
+# The Hazard Library
+# Copyright (C) 2012-2023 GEM Foundation
+#
+# This program 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.
+#
+# This program 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 this program.  If not, see <http://www.gnu.org/licenses/>.
+
+"""
+Module :mod:`openquake.hazardlib.mgmpe.cb14_basin_term` implements
+:class:`~openquake.hazardlib.mgmpe.CB14BasinTerm`
+"""
+import numpy as np
+
+from openquake.hazardlib import const
+from openquake.hazardlib.gsim.base import GMPE, registry
+from openquake.hazardlib.gsim.campbell_bozorgnia_2014 import CampbellBozorgnia2014
+
+
+def _get_cb14_basin_term(ctx, C, jpn_flag=False):
+    """
+    Get the basin response term defined in equation 20 of the Campbell and 
+    Bozorgnia (2014) GMM paper.
+
+    Currently the global basin term is provided (i.e. the Japan-regionalised
+    basin term is for now turned off).
+    """
+    z2pt5 = ctx.z2pt5
+    fb = np.zeros(len(z2pt5))
+    idx = z2pt5 < 1.0
+    fb[idx] = (C["c14"] + C["c15"] * jpn_flag) * (z2pt5[idx] - 1.0)
+    idx = z2pt5 > 3.0
+    fb[idx] = C["c16"] * C["k3"] * np.exp(-0.75) * (
+        1. - np.exp(-0.25 * (z2pt5[idx] - 3.)))
+
+    return fb
+
+
+class CB14BasinTerm(GMPE):
+    """
+    Implements a modified GMPE class that can be used to implement the Campbell
+    and Bozorgnia (2014) GMM's basin term
+
+    :param gmpe_name:
+        The name of a GMPE class
+    """
+    # Req Params
+    REQUIRES_SITES_PARAMETERS = {'z2pt5'}
+
+    # Others are set from underlying GMM
+    REQUIRES_DISTANCES = set() 
+    REQUIRES_RUPTURE_PARAMETERS = set()
+    DEFINED_FOR_INTENSITY_MEASURE_COMPONENT = ""
+    DEFINED_FOR_INTENSITY_MEASURE_TYPES = set()
+    DEFINED_FOR_STANDARD_DEVIATION_TYPES = {const.StdDev.TOTAL}
+    DEFINED_FOR_TECTONIC_REGION_TYPE = ""
+    DEFINED_FOR_REFERENCE_VELOCITY = None
+
+    def __init__(self, gmpe_name, **kwargs):
+        self.gmpe = registry[gmpe_name]()
+        self.set_parameters()    
+
+        # Need z2pt5 in req site params to ensure in ctx site col 
+        if 'z2pt5' not in self.gmpe.REQUIRES_SITES_PARAMETERS:
+            self.REQUIRES_SITES_PARAMETERS = frozenset(
+                self.gmpe.REQUIRES_SITES_PARAMETERS | {'z2pt5'})
+
+    def compute(self, ctx: np.recarray, imts, mean, sig, tau, phi):      
+        """
+        See :meth:`superclass method
+        <.base.GroundShakingIntensityModel.compute>`
+        for spec of input and result values.
+        """
+        self.gmpe.compute(ctx, imts, mean, sig, tau, phi)
+        for m, imt in enumerate(imts):
+             C = CampbellBozorgnia2014.COEFFS[imt]
+             mean[m] += _get_cb14_basin_term(ctx, C)

openquake/hazardlib/gsim/mgmpe/cy14_site_term.py

@@ -26,7 +26,7 @@
 from openquake.hazardlib.gsim.chiou_youngs_2014 import ChiouYoungs2014
 
 
-def _get_site_term(C, vs30, ln_y_ref):
+def _get_cy14_site_term(C, vs30, ln_y_ref):
     """
     Applies the linear and nonlinear site amplification term of Chiou &
     Youngs (2014) (excluding the basin amplification term)
@@ -104,4 +104,4 @@ def compute(self, ctx: np.recarray, imts, mean, sig, tau, phi):
         vs30 = ctx.vs30.copy()
         for m, imt in enumerate(imts):
             C = ChiouYoungs2014.COEFFS[imt]
-            mean[m] += _get_site_term(C, vs30, mean[m])
+            mean[m] += _get_cy14_site_term(C, vs30, mean[m])

openquake/hazardlib/gsim/mgmpe/m9_basin_term.py

@@ -0,0 +1,77 @@
+# The Hazard Library
+# Copyright (C) 2012-2023 GEM Foundation
+#
+# This program 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.
+#
+# This program 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 this program.  If not, see <http://www.gnu.org/licenses/>.
+
+"""
+Module :mod:`openquake.hazardlib.mgmpe.m9_basin_term` implements
+:class:`~openquake.hazardlib.mgmpe.M9BasinTerm`
+"""
+import numpy as np
+
+from openquake.hazardlib import const
+from openquake.hazardlib.gsim.base import GMPE, registry
+
+
+def _apply_m9_basin_term(ctx, imt, mean):
+    if imt.period > 1.9: # Only apply to long-period SA
+        fb_m9 = np.log(2.0)
+        idx = ctx.z2pt5 >= 6.0 # Apply only to sites with z2pt5 >= 6
+        mean[idx] += fb_m9
+
+    return mean
+
+
+class M9BasinTerm(GMPE):
+    """
+    Implements a modified GMPE class that can be used to implement the "M9"
+    US 2023 NSHM basin amplification adjustment.
+
+    This implementation is based on the description of the M9 adjustment 
+    within the Moschetti et al. (2024) EQ Spectra article on the conterminous
+    US 2023 NSHM GMC (pp. 1178).
+
+    :param gmpe_name:
+        The name of a GMPE class
+    """
+    # Req Params
+    REQUIRES_SITES_PARAMETERS = {'z2pt5'}
+
+    # Others are set from underlying GMM
+    REQUIRES_DISTANCES = set() 
+    REQUIRES_RUPTURE_PARAMETERS = set()
+    DEFINED_FOR_INTENSITY_MEASURE_COMPONENT = ""
+    DEFINED_FOR_INTENSITY_MEASURE_TYPES = set()
+    DEFINED_FOR_STANDARD_DEVIATION_TYPES = {const.StdDev.TOTAL}
+    DEFINED_FOR_TECTONIC_REGION_TYPE = ""
+    DEFINED_FOR_REFERENCE_VELOCITY = None
+
+    def __init__(self, gmpe_name, **kwargs):
+        self.gmpe = registry[gmpe_name]()
+        self.set_parameters()    
+
+        # Need z2pt5 in req site params to ensure in ctx site col 
+        if 'z2pt5' not in self.gmpe.REQUIRES_SITES_PARAMETERS:
+            self.REQUIRES_SITES_PARAMETERS = frozenset(
+                self.gmpe.REQUIRES_SITES_PARAMETERS | {'z2pt5'})
+
+    def compute(self, ctx: np.recarray, imts, mean, sig, tau, phi):      
+        """
+        See :meth:`superclass method
+        <.base.GroundShakingIntensityModel.compute>`
+        for spec of input and result values.
+        """
+        self.gmpe.compute(ctx, imts, mean, sig, tau, phi)
+        for m, imt in enumerate(imts):
+             mean[m] = _apply_m9_basin_term(ctx, imt, mean[m])

openquake/hazardlib/gsim/mgmpe/modifiable_gmpe.py

@@ -26,8 +26,11 @@
 from openquake.hazardlib.imt import from_string
 from openquake.hazardlib.gsim.mgmpe.nrcan15_site_term import (
     NRCan15SiteTerm, BA08_AB06)
-from openquake.hazardlib.gsim.mgmpe.cy14_site_term import _get_site_term
+from openquake.hazardlib.gsim.mgmpe.cy14_site_term import _get_cy14_site_term
 from openquake.hazardlib.gsim.chiou_youngs_2014 import ChiouYoungs2014
+from openquake.hazardlib.gsim.mgmpe.cb14_basin_term import _get_cb14_basin_term
+from openquake.hazardlib.gsim.campbell_bozorgnia_2014 import CampbellBozorgnia2014
+from openquake.hazardlib.gsim.mgmpe.m9_basin_term import _apply_m9_basin_term
 
 from openquake.hazardlib.gsim.nga_east import (
     TAU_EXECUTION, get_phi_ss, TAU_SETUP, PHI_SETUP, get_tau_at_quantile,
@@ -73,10 +76,25 @@ def cy14_site_term(ctx, imt, me, si, ta, phi):
     This function adds the CY14 site term to GMMs requiring it
     """
     C = ChiouYoungs2014.COEFFS[imt]
-    fa = _get_site_term(C, ctx.vs30, me)  # ref mean must be in natural log
+    fa = _get_cy14_site_term(C, ctx.vs30, me) # Ref mean must be in natural log
     me[:] += fa
 
 
+def cb14_basin_term(ctx, imt, me, si, ta, phi):
+    """
+    This function adds the CB14 basin term to GMMs requiring it.
+    """
+    C = CampbellBozorgnia2014.COEFFS[imt]
+    me[:] += _get_cb14_basin_term(ctx, C)
+
+
+def m9_basin_term(ctx, imt, me, si, ta, phi):
+    """
+    This function applies the M9 basin adjustment
+    """
+    me = _apply_m9_basin_term(ctx, imt, me)
+
+
 def add_between_within_stds(ctx, imt, me, si, ta, ph, with_betw_ratio):
     """
     This adds the between and within standard deviations to a model which has
@@ -248,6 +266,12 @@ def __init__(self, **kwargs):
             setattr(self, 'DEFINED_FOR_STANDARD_DEVIATION_TYPES',
                     {StdDev.TOTAL, StdDev.INTRA_EVENT, StdDev.INTER_EVENT})
 
+        if ('cb14_basin_term' in self.params or 'm9_basin_term' in self.params
+            ) and ( 'z2pt5' not in self.gmpe.REQUIRES_SITES_PARAMETERS):
+            tmp = list(self.gmpe.REQUIRES_SITES_PARAMETERS)
+            tmp.append('z2pt5')
+            self.gmpe.REQUIRES_SITES_PARAMETERS = frozenset(tmp)
+
         # This is required by the `sigma_model_alatik2015` function
         key = 'sigma_model_alatik2015'
         if key in self.params:
@@ -301,14 +325,15 @@ def compute(self, ctx: np.recarray, imts, mean, sig, tau, phi):
         <.base.GroundShakingIntensityModel.compute>`
         for spec of input and result values.
         """
+        # Set reference Vs30 if required
         if ('nrcan15_site_term' in self.params or
                 'cy14_site_term' in self.params):
             ctx_copy = ctx.copy()
             if 'nrcan15_site_term' in self.params:
                 rock_vs30 = 760.
             elif 'cy14_site_term' in self.params:
                 rock_vs30 = 1130.
-            ctx_copy.vs30 = np.full_like(ctx.vs30, rock_vs30)  # rock
+            ctx_copy.vs30 = np.full_like(ctx.vs30, rock_vs30) # rock
         else:
             ctx_copy = ctx
         g = globals()

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

@@ -0,0 +1,89 @@
+# The Hazard Library
+# Copyright (C) 2012-2023 GEM Foundation
+#
+# This program 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.
+#
+# This program 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 this program.  If not, see <http://www.gnu.org/licenses/>.
+
+import numpy as np
+import unittest
+
+from openquake.hazardlib.tests.gsim.mgmpe.dummy import new_ctx
+from openquake.hazardlib.contexts import simple_cmaker
+from openquake.hazardlib import valid
+from openquake.hazardlib.imt import PGA, SA
+from openquake.hazardlib.const import TRT, IMC
+from openquake.hazardlib.gsim.mgmpe.cb14_basin_term import CB14BasinTerm
+
+ae = np.testing.assert_equal
+aae = np.testing.assert_almost_equal
+
+exp_gmm_origin = np.array([[-1.55067862, -3.8245714 , -5.48623229],
+                           [-2.25723424, -3.72056054, -5.28051045],
+                           [-3.45764676, -4.67077352, -5.84789307]])
+
+exp_with_basin = np.array([[-1.35815033, -3.78220564, -5.293704  ],
+                           [-1.86104115, -3.63337845, -4.88431736],
+                           [-3.04727152, -4.58047065, -5.43751783]])
+
+class M9BasinTermTestCase(unittest.TestCase):
+
+    def test_instantiation(self):
+        mgmpe = CB14BasinTerm(gmpe_name='AtkinsonMacias2009')
+
+        # Check the assigned IMTs
+        expected = set([PGA, SA])
+        ae(mgmpe.DEFINED_FOR_INTENSITY_MEASURE_TYPES, expected)
+        # Check the TR
+        expected = TRT.SUBDUCTION_INTERFACE
+        ae(mgmpe.DEFINED_FOR_TECTONIC_REGION_TYPE, expected)
+        # Check the IM component
+        expected = IMC.RANDOM_HORIZONTAL
+        ae(mgmpe.DEFINED_FOR_INTENSITY_MEASURE_COMPONENT, expected)
+        # Check the required distances
+        expected = {'rrup'}
+        ae(mgmpe.REQUIRES_DISTANCES, expected)
+
+    def test_all(self):
+        """
+        Test that the M9 basin term applied to Kuehn et al. (2020)
+        provides the expected values (using sites with z2pt5 above
+        and below 6 km threshold and considering SAs with periods
+        above and below 1.9 s)
+        """
+        am09 = valid.gsim('AtkinsonMacias2009')
+
+        # Make original GMM
+        gmpe = valid.modified_gsim(am09)
+
+        # Make GMM with basin term using the mgmpe class directly
+        mgmpe_cls = CB14BasinTerm(gmpe_name='AtkinsonMacias2009')
+
+        # Make GMM with basin term using ModifiableGMPE and kwargs
+        mgmpe_val = valid.modified_gsim(am09, cb14_basin_term={})
+        # Make the ctx
+        imts = ['PGA', 'SA(1.0)', 'SA(2.0)']
+        cmaker = simple_cmaker([gmpe, mgmpe_cls, mgmpe_val], imts)                       
+        ctx = new_ctx(cmaker, 3)
+        ctx.dip = 60.
+        ctx.rake = 90.
+        ctx.z1pt0 = np.array([522.32, 516.98, 522.32])
+        ctx.z2pt5 = np.array([6.32, 3.53, 6.32])
+        ctx.rrup = np.array([50., 200., 500.])
+        ctx.vs30 = 1100.
+        ctx.vs30measured = 1
+        mea, _, _, _ = cmaker.get_mean_stds([ctx])
+
+        # Check expected is observed (original vs modified with basin term)
+        aae(mea[0], exp_gmm_origin)
+        aae(mea[1], exp_with_basin)
+        aae(mea[2], exp_with_basin)