Source code for openquake.hazardlib.gsim.si_midorikawa_1999

# -*- coding: utf-8 -*-
# vim: tabstop=4 shiftwidth=4 softtabstop=4
#
# Copyright (C) 2013-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:`SiMidorikawa1999Asc`, :class:`SiMidorikawa1999SInter`,
:class:`SiMidorikawaSSlab`, :class:`SiMidorikawa1999SInterNorthEastCorrection`,
:class:`SiMidorikawa1999SInterSouthWestCorrection`,
:class:`SiMidorikawa1999SSlabNorthEastCorrection` and
:class:`SiMidorikawa1999SSlabSouthWestCorrection`.
"""
import numpy as np
from scipy.constants import g
from openquake.baselib.general import CallableDict
from openquake.hazardlib.gsim.base import GMPE
from openquake.hazardlib import const
from openquake.hazardlib.imt import PGV, PGA
from openquake.hazardlib.geo import (
    Mesh, RectangularMesh, SimpleFaultSurface)


# Subduction trench coordinates (table 3.5.2-1) page 3-150
SUB_TRENCH_LONS = np.array([
    143.50, 143.00, 141.90, 142.40, 143.25, 143.80, 144.20, 144.30, 144.65,
    146.80, 153.00])

SUB_TRENCH_LATS = np.array([
    24.00, 29.00, 33.80, 35.80, 36.55, 37.70, 39.20, 40.10, 41.00, 42.00, 45.50
])

# Volcanic front coordinates (table 3.5.2-2, page 3-150)
VOLCANIC_FRONT_LONS = np.array([
    122.00, 124.00, 128.30, 129.70, 130.80, 131.60, 132.00, 133.70, 134.90,
    136.90])

VOLCANIC_FRONT_LATS = np.array([
    24.50, 24.50, 27.90, 29.50, 31.50, 33.40, 34.90, 35.30, 35.30, 36.20])


def _construct_surface(lons, lats, upper_depth, lower_depth):
    """
    Utility method that constructs and return a simple fault surface with top
    edge specified by `lons` and `lats` and extending vertically from
    `upper_depth` to `lower_depth`.

    The underlying mesh is built by repeating the same coordinates
    (`lons` and `lats`) at the two specified depth levels.
    """
    depths = np.array([
        np.zeros_like(lons) + upper_depth,
        np.zeros_like(lats) + lower_depth
    ])

    mesh = RectangularMesh(
        np.tile(lons, (2, 1)), np.tile(lats, (2, 1)), depths
    )
    return SimpleFaultSurface(mesh)


def _get_min_distance_to_sub_trench(lons, lats):
    """
    Compute and return minimum distance between subduction trench
    and points specified by 'lon' and 'lat'

    The method creates an instance of
    :class:`openquake.hazardlib.geo.SimpleFaultSurface` to model the subduction
    trench. The surface is assumed vertical and extending from 0 to 10 km
    depth.
    The 10 km depth value is arbitrary given that distance calculation depend
    only on top edge depth. The method calls then
    :meth:`openquake.hazardlib.geo.base.BaseSurface.get_rx_distance`
    and return its absolute value.
    """
    trench = _construct_surface(SUB_TRENCH_LONS, SUB_TRENCH_LATS, 0., 10.)
    ctx = Mesh(lons, lats, None)
    return np.abs(trench.get_rx_distance(ctx))


def _get_min_distance_to_volcanic_front(lons, lats):
    """
    Compute and return minimum distance between volcanic front and points
    specified by 'lon' and 'lat'.

    Distance is negative if point is located east of the volcanic front,
    positive otherwise.

    The method uses the same approach as
    :meth:`_get_min_distance_to_sub_trench`
    but final distance is returned without taking the absolute value.
    """
    vf = _construct_surface(VOLCANIC_FRONT_LONS, VOLCANIC_FRONT_LATS, 0., 10.)
    ctx = Mesh(lons, lats, None)
    return vf.get_rx_distance(ctx)


def _apply_subduction_trench_correction(mean, x_tr, H, rrup, imt):
    """
    Implement equation for subduction trench correction as described in
    equation 3.5.2-1, page 3-148 of "Technical Reports on National Seismic
    Hazard Maps for Japan"
    """
    if imt.string == 'PGV':
        corr = np.maximum(1., 10 ** -0.012 * (rrup / 300.) ** 2.064)
        corr[H > 30] *= 10 ** ((-4.021E-5 * x_tr[H > 30] + 9.905e-3) *
                               (H[H > 30] - 30))
    else:
        corr = np.maximum(1., 10 ** 0.13 * (rrup / 300.) ** 3.2)
        corr[H > 30] *= 10 ** ((-8.1E-5 * x_tr[H > 30] + 2.0e-2) *
                               (H[H > 30] - 30))
    return np.log(np.exp(mean) * corr)


def _apply_volcanic_front_correction(mean, x_vf, H, imt):
    """
    Implement equation for volcanic front correction as described in equation
    3.5.2.-2, page 3-149 of "Technical Reports on National Seismic
    Hazard Maps for Japan"
    """
    V1 = np.zeros_like(x_vf)
    if imt.string == 'PGV':
        idx = x_vf <= 75
        V1[idx] = 4.28e-5 * x_vf[idx] * (H[idx] - 30)
        idx = x_vf > 75
        V1[idx] = 3.21e-3 * (H[idx] - 30)
    else:
        idx = x_vf <= 75
        V1[idx] = 7.06e-5 * x_vf[idx] * (H[idx] - 30)
        idx = x_vf > 75
        V1[idx] = 5.30e-3 * (H[idx] - 30)
    return np.log(np.exp(mean) * 10 ** V1)


def _apply_amplification_factor(AMP_F, mean, vs30):
    """
    Apply amplification factor to scale PGV value from 600 to 400 m/s vs30
    and convert mean from base 10 to base e.

    The scaling factor from 600 m/s to 400 m/s was defined by NIED.

    The scaling factor from 600 m/s to 800m/s is valid just for the elastic
    case as no adjustment for kappa was considered.
    """
    assert np.all(vs30 == vs30[0])
    if abs(vs30[0]-600.) < 1e-10:
        return mean * np.log(10)
    elif abs(vs30[0]-400.) < 1e-10:
        return mean * np.log(10) + np.log(AMP_F)
    elif abs(vs30[0]-800.) < 1e-10:
        return mean * np.log(10) - np.log(1.25)
    else:
        raise ValueError(
            'Si and Midorikawa 1999 do not support this Vs30 value')


[docs]def set_mean(imt, mag, hypo_depth, rrup, d, mean):
    """
    Return mean value as defined in equation 3.5.1-1 page 148
    """
    # clip magnitude at 8.3 as per note at page 3-36 in table Table 3.3.2-6
    # in "Technical Reports on National Seismic Hazard Maps for Japan"
    mag = np.fmin(mag, 8.3)
    if imt.string == 'PGV':
        mean[:] = (
            0.58 * mag +
            0.0038 * hypo_depth +
            d -
            1.29 -
            np.log10(rrup + 0.0028 * 10 ** (0.5 * mag)) -
            0.002 * rrup)
    else:
        mean[:] = (
            0.50 * mag +
            0.0043 * hypo_depth +
            d +
            0.61 -
            np.log10(rrup + 0.0055 * 10 ** (0.5 * mag)) -
            0.003 * rrup)
        mean[:] = np.log10(10**mean / (g*100))


set_sig = CallableDict()


[docs]@set_sig.add(const.TRT.ACTIVE_SHALLOW_CRUST)
def set_sig_asc(trt, rrup, sig):
    """
    Set standard deviation as defined in equation 3.5.5-2 page 151
    """
    sig[rrup <= 20] = 0.23
    idx = (rrup > 20) & (rrup <= 30)
    sig[idx] = 0.23 - 0.03 * np.log10(rrup[idx] / 20) / np.log10(30. / 20.)
    sig[rrup > 30] = 0.20
    # convert from log10 to ln
    sig[:] = np.log(10 ** sig)


[docs]@set_sig.add(const.TRT.SUBDUCTION_INTERFACE,
             const.TRT.SUBDUCTION_INTRASLAB)
def set_sig_sub(trt, pgv, sig):
    """
    Set standard deviation as defined in equation 3.5.5-1 page 151
    """
    sig[pgv <= 25] = 0.20
    idx = (pgv > 25) & (pgv <= 50)
    sig[idx] = 0.20 - 0.05 * (pgv[idx] - 25) / 25
    sig[pgv > 50] = 0.15
    # convert from log10 to ln
    sig[:] = np.log(10 ** sig)


[docs]class SiMidorikawa1999Asc(GMPE):
    """
    Implements GMPE developed by Hongjun Si and Saburoh Midorikawa (1999) as
    described in "Technical Reports on National Seismic Hazard Maps for Japan"
    (2009, National Research Institute for Earth Science and Disaster
    Prevention, Japan, pages 148-151).
    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
    DEFINED_FOR_TECTONIC_REGION_TYPE = const.TRT.ACTIVE_SHALLOW_CRUST

    #: Supported intensity measure type is PGV
    DEFINED_FOR_INTENSITY_MEASURE_TYPES = {PGV, PGA}

    #: Supported intensity measure component is greater of
    #: of two horizontal components :
    #: attr:`~openquake.hazardlib.const.IMC.GREATER_OF_TWO_HORIZONTAL`
    DEFINED_FOR_INTENSITY_MEASURE_COMPONENT = \
        const.IMC.GREATER_OF_TWO_HORIZONTAL

    #: Supported standard deviation type is total
    DEFINED_FOR_STANDARD_DEVIATION_TYPES = {const.StdDev.TOTAL}

    #: No ctx parameters are required
    # REQUIRES_SITES_PARAMETERS = set()
    REQUIRES_SITES_PARAMETERS = {'vs30'}

    #: Required rupture parameters are magnitude, and hypocentral depth
    REQUIRES_RUPTURE_PARAMETERS = {'mag', 'hypo_depth'}

    #: Required distance measure is Rrup
    REQUIRES_DISTANCES = {'rrup'}

    #: Amplification factor to scale PGV from 600 to 400 m/s vs30,
    #: see equation 3.5.1-1 page 148
    AMP_F = 1.41

[docs]    def compute(self, ctx: np.recarray, imts, mean, sig, tau, phi):
        """
        Implements equation 3.5.1-1 page 148 for mean value and equation
        3.5.5-2 page 151 for total standard deviation.

        See :meth:`superclass method
        <.base.GroundShakingIntensityModel.compute>`
        for spec of input and result values.
        """
        trt = self.DEFINED_FOR_TECTONIC_REGION_TYPE
        for m, imt in enumerate(imts):
            set_mean(imt, ctx.mag, ctx.hypo_depth, ctx.rrup, 0, mean[m])
            set_sig(trt, ctx.rrup, sig[m])
            mean[m] = _apply_amplification_factor(
                self.AMP_F, mean[m], ctx.vs30)


[docs]class SiMidorikawa1999SInter(SiMidorikawa1999Asc):
    """
    Implements GMPE developed by Hongjun Si and Saburoh Midorikawa (1999) as
    described in "Technical Reports on National Seismic Hazard Maps for Japan"
    (2009, National Research Institute for Earth Science and Disaster
    Prevention, Japan, pages 148-151).
    This class implements the equations for 'Subduction Interface'
    (that's why the class name ends with 'SInter').
    """
    #: Supported tectonic region type is subduction interface
    DEFINED_FOR_TECTONIC_REGION_TYPE = const.TRT.SUBDUCTION_INTERFACE

[docs]    def compute(self, ctx: np.recarray, imts, mean, sig, tau, phi):
        """
        Implements equation 3.5.1-1 page 148 for mean value and equation
        3.5.5-1 page 151 for total standard deviation.

        See :meth:`superclass method
        <.base.GroundShakingIntensityModel.compute>`
        for spec of input and result values.
        """
        trt = self.DEFINED_FOR_TECTONIC_REGION_TYPE
        for m, imt in enumerate(imts):
            if imt.string == 'PGV':
                d = -0.02
            else:
                d = 0.01
            set_mean(imt, ctx.mag, ctx.hypo_depth, ctx.rrup, d, mean[m])
            set_sig(trt, 10 ** mean[m], sig[m])
            mean[m] = _apply_amplification_factor(
                self.AMP_F, mean[m], ctx.vs30)


[docs]class SiMidorikawa1999SInterNorthEastCorrection(SiMidorikawa1999SInter):
    """
    Extend :class:`SiMidorikawa1999SInter` and takes into account
    correction for northeast Japan (i.e. proximity to subduction trench)
    """
    REQUIRES_SITES_PARAMETERS = {'lon', 'lat', 'vs30'}

[docs]    def compute(self, ctx: np.recarray, imts, mean, sig, tau, phi):
        """
        Implements equation 3.5.1-1 page 148 for mean value and equation
        3.5.5-1 page 151 for total standard deviation.

        See :meth:`superclass method
        <.base.GroundShakingIntensityModel.compute>`
        for spec of input and result values.
        """
        super().compute(ctx, imts, mean, sig, tau, phi)
        x_tr = _get_min_distance_to_sub_trench(ctx.lon, ctx.lat)
        for m, imt in enumerate(imts):
            mean[m] = _apply_subduction_trench_correction(
                mean[m], x_tr, ctx.hypo_depth, ctx.rrup, imt)


[docs]class SiMidorikawa1999SInterSouthWestCorrection(SiMidorikawa1999SInter):
    """
    Extend :class:`SiMidorikawa1999SInter` and takes into account
    correction for southwest Japan (i.e. proximity with volcanic front)
    """
    REQUIRES_SITES_PARAMETERS = {'lon', 'lat', 'vs30'}

[docs]    def compute(self, ctx: np.recarray, imts, mean, sig, tau, phi):
        """
        Implements equation 3.5.1-1 page 148 for mean value and equation
        3.5.5-1 page 151 for total standard deviation.

        See :meth:`superclass method
        <.base.GroundShakingIntensityModel.compute>`
        for spec of input and result values.
        """
        super().compute(ctx, imts, mean, sig, tau, phi)
        x_vf = _get_min_distance_to_volcanic_front(ctx.lon, ctx.lat)
        for m, imt in enumerate(imts):
            mean[m] = _apply_volcanic_front_correction(
                mean[m], x_vf, ctx.hypo_depth, imt)


[docs]class SiMidorikawa1999SSlab(SiMidorikawa1999SInter):
    """
    Implements GMPE developed by Hongjun Si and Saburoh Midorikawa (1999) as
    described in "Technical Reports on National Seismic Hazard Maps for Japan"
    (2009, National Research Institute for Earth Science and Disaster
    Prevention, Japan, pages 148-151).
    This class implements the equations for 'Subduction IntraSlab'
    (that's why the class name ends with 'SSlab').
    """
    #: Supported tectonic region type is subduction intraslab
    DEFINED_FOR_TECTONIC_REGION_TYPE = const.TRT.SUBDUCTION_INTRASLAB

[docs]    def compute(self, ctx: np.recarray, imts, mean, sig, tau, phi):
        """
        Implements equation 3.5.1-1 page 148 for mean value and equation
        3.5.5-1 page 151 for total standard deviation.

        See :meth:`superclass method
        <.base.GroundShakingIntensityModel.compute>`
        for spec of input and result values.
        """
        trt = self.DEFINED_FOR_TECTONIC_REGION_TYPE
        for m, imt in enumerate(imts):
            if imt.string == 'PGV':
                d = 0.12
            else:
                d = 0.22
            set_mean(imt, ctx.mag, ctx.hypo_depth, ctx.rrup, d, mean[m])
            set_sig(trt, 10 ** mean[m], sig[m])
            mean[m] = _apply_amplification_factor(
                self.AMP_F, mean[m], ctx.vs30)


[docs]class SiMidorikawa1999SSlabNorthEastCorrection(SiMidorikawa1999SSlab):
    """
    Extend :class:`SiMidorikawa1999SSlab` and takes into account
    correction for northeast Japan (i.e. proximity to subduction trench)
    """
    REQUIRES_SITES_PARAMETERS = {'lon', 'lat', 'vs30'}

[docs]    def compute(self, ctx: np.recarray, imts, mean, sig, tau, phi):
        """
        Implements equation 3.5.1-1 page 148 for mean value and equation
        3.5.5-1 page 151 for total standard deviation.

        See :meth:`superclass method
        <.base.GroundShakingIntensityModel.compute>`
        for spec of input and result values.
        """
        super().compute(ctx, imts, mean, sig, tau, phi)
        x_tr = _get_min_distance_to_sub_trench(ctx.lon, ctx.lat)
        for m, imt in enumerate(imts):
            mean[m] = _apply_subduction_trench_correction(
                mean[m], x_tr, ctx.hypo_depth, ctx.rrup, imt)


[docs]class SiMidorikawa1999SSlabSouthWestCorrection(SiMidorikawa1999SSlab):
    """
    Extend :class:`SiMidorikawa1999SSlab` and takes into account
    correction for southwest Japan (i.e. proximity to volcanic front)
    """
    REQUIRES_SITES_PARAMETERS = {'lon', 'lat', 'vs30'}

[docs]    def compute(self, ctx: np.recarray, imts, mean, sig, tau, phi):
        """
        Implements equation 3.5.1-1 page 148 for mean value and equation
        3.5.5-1 page 151 for total standard deviation.

        See :meth:`superclass method
        <.base.GroundShakingIntensityModel.compute>`
        for spec of input and result values.
        """
        super().compute(ctx, imts, mean, sig, tau, phi)
        x_vf = _get_min_distance_to_volcanic_front(ctx.lon, ctx.lat)
        for m, imt in enumerate(imts):
            mean[m] = _apply_volcanic_front_correction(
                mean[m], x_vf, ctx.hypo_depth, imt)