Source code for openquake.hazardlib.gsim.nz22.atkinson_2022

# -*- coding: utf-8 -*-
# vim: tabstop=4 shiftwidth=4 softtabstop=4
#
# Copyright (C) 2015-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:`Atkinson2022Crust`
               :class:`Atkinson2022SSlab`
               :class:`Atkinson2022SInter`
"""

from pathlib import Path

import numpy as np
from scipy.interpolate import interp1d

from openquake.hazardlib import const
from openquake.hazardlib.gsim.base import GMPE, CoeffsTable
from openquake.hazardlib.gsim.nz22.const import (
    periods_AG20,
    rho_Ws,
    rho_Bs,
    periods,
    theta7s,
    theta8s,
)
from openquake.hazardlib.imt import PGA, SA

Atk22_COEFFS = Path(
    Path(__file__).parent, "Atkinson22_coeffs_mod_v8b_sanjay_v2.csv"
)


def _fmag(suffix, C, mag):
    """
    ctx.magnitude factor
    """
    if suffix == "slab":
        # res = C['c0_' + suffix] + C['c1_' + suffix] * (mag - 6.0) + C['c2_' + suffix] * (mag - 6.0) ** 2
        # Modified as in RevisionsToBackbonev8 from Gail received on 21.06.2022.
        res = (
            C["c0_" + suffix]
            + C["c1_" + suffix] * mag
            + C["c2_" + suffix] * mag**2
        )
    else:
        res = C["c0_crust"] + C["c1_crust"] * mag + C["c2_crust"] * mag**2
    return res


def _fz_ha18(C, ctx):
    """
    Implements eq. 2,3,4,5 from page 5
    """
    # pseudo-depth

    # h = 10 ** (-0.1 + 0.2 * ctx.mag) The h term is modified after
    # receiving the modifications from Gail on Slack on 12.06.2022.  h
    # = 10 ** (0.3 + 0.11 * ctx.mag) Modified as in
    # RevisionsToBackbonev8 from Gail received on 21.06.2022. However,
    # there is a typo.
    h = 10 ** (-0.405 + 0.235 * ctx.mag)
    R = np.sqrt(ctx.rrup**2 + h**2)
    Rref = np.sqrt(1 + h**2)
    # The transition_distance
    Rt = 50
    # Geometrical spreading rates
    b1 = -1.3
    b2 = -0.5
    # Geometrical attenuation
    z = R**b1
    ratio = R / Rt
    z[R > Rt] = Rt**b1 * (ratio[R > Rt]) ** b2

    return np.log(z) + (C["b3"] + C["b4"] * ctx.mag) * np.log(R / Rref)


def _fgamma(suffix, C, ctx):
    if suffix == "crust":
        g1 = min(0.008, 0.005 + 0.0016 * np.log(C["f"]))
    elif suffix == "inter":
        g1 = min(0.006, 0.0045 + 0.0014 * np.log(C["f"]))
    else:
        g1 = min(0.005, 0.004 + 0.0012 * np.log(C["f"]))

    # a2 = max(0.002 + 0.0025 * np.log(max(C['f'], 35)), 0.0015)
    # a3 = 0.009 - 0.001 * np.log(max(C['f'], 35))
    # g2 = min(min(0.0065, a2), a3)

    gamma = np.zeros(ctx.rrup.shape)

    # gamma = -g1 * ctx.rrup + g2*(270.0 - np.clip(ctx.rrup, 270,
    # None)) Gail mentioned in personal communication (email
    # 13.06.2022) that now the modified F_gamma (see eq. 20-22 in
    # modifications posted on Salck 12.06.2022) does not include
    # gamma_2 term.
    gamma = -g1 * ctx.rrup
    return gamma


def _epistemic_adjustment_lower(C, ctx):
    # These are revised adjustments after Gail's post on slack 11th
    # May 2022 and in her revised report.  The lower branch adjustment
    # remains the same.  a = np.fmax(np.clip(0.5 - 0.1 *
    # np.log(ctx.rrup), 0.2, None), - 0.25 + 0.1 * np.log(ctx.rrup))
    # The following variable is after Gail's modifications received on
    # Slack 12.06.2022 The additional epistemic uncertainty for M>7
    # events was added in Gail's V8 modifications shared on 27.06.2022
    a = np.fmax(
        np.clip(0.6 - 0.13 * np.log(ctx.rrup), 0.3, None),
        -0.25 + 0.12 * np.log(ctx.rrup),
    )
    return np.clip(a, -np.inf, 0.5) + 0.15 * np.clip(ctx.mag - 7.0, 0, np.inf)


def _epistemic_adjustment_upper(C, ctx):
    # These are revised adjustments after Gail's post on slack 11th
    # May 2022 and in her revised report.  Only the upper brach is
    # modified.  a = np.fmax(np.clip(1.0 - 0.27 * np.log(ctx.rrup),
    # 0.2, None), - 0.25 + 0.1 * np.log(ctx.rrup)) The following
    # variable is after Modification from Gail recieved on Slack
    # 12.06.2022 The additional epistemic uncertainty for M>7 events
    # was added in Gail's V8 modifications shared on 27.06.2022
    a = np.fmax(
        np.clip(1.0 - 0.3 * np.log(ctx.rrup), 0.3, None),
        -0.25 + 0.12 * np.log(ctx.rrup),
    )
    return np.clip(a, -np.inf, 0.8) + 0.15 * np.clip(ctx.mag - 7.0, 0, np.inf)


[docs]def fs_SS14(C, pga_rock, ctx):
    # The site-term is implemnted from Seyhan and Stewart (2014; EQS).
    Vs_ref = 760.0
    flin = ctx.vs30 / Vs_ref
    flin[ctx.vs30 > C["Vc"]] = C["Vc"] / Vs_ref
    flin_func = C["c"] * np.log(flin)

    v_s = np.copy(ctx.vs30)
    v_s[ctx.vs30 > 760.0] = 760.0
    f_1 = 0.0
    f_3 = 0.1 * 981.0  # In Gail's model the GMM is in cm/s^2.
    # Nonlinear controlling parameter (equation 8)
    f_2 = C["f4"] * (np.exp(C["f5"] * (v_s - 360.0)) - np.exp(C["f5"] * 400.0))
    fnl = f_1 + f_2 * np.log((pga_rock + f_3) / f_3)
    return flin_func + fnl


def _get_pga_on_rock(suffix, C, ctx):
    pga_rock = np.exp(
        _fmag(suffix, C, ctx.mag) + _fz_ha18(C, ctx) + _fgamma(suffix, C, ctx)
    )
    return pga_rock


[docs]def get_stddevs(suffix, C):
    r"""
    Standard deviations given in COEFFS with suffix
    between event standard deviations as Be\_
    within event stdvs as We\_
    total as sigma\_
    """
    intra_e_sigma = np.sqrt(C["We_" + suffix] ** 2 + C["phiS2S"] ** 2)
    return [C["sigma_" + suffix], C["Be_" + suffix], intra_e_sigma]


[docs]def get_nonlinear_stddevs(suffix, C, C_PGA, imt, pga_rock, vs30):
    """
    Get the nonlinear tau and phi terms for Gail's model. For this
    implementation, I using the between-event and within-event
    correlation from AG20. Note that the Gail's nonlinear soil
    response term is identical to Seyhan and Stewart (2014) model.
    """
    period = imt.period
    pgar = pga_rock / 981.0
    # Linear Tau
    tau_lin = C["Be_" + suffix] * np.ones(vs30.shape)
    tau_lin_pga = C_PGA["Be_" + suffix] * np.ones(vs30.shape)

    # Linear phi
    intra_e_sigma = np.sqrt(C["We_" + suffix] ** 2 + C["phiS2S"] ** 2)
    intra_e_sigma_pga = np.sqrt(
        C_PGA["We_" + suffix] ** 2 + C_PGA["phiS2S"] ** 2
    )
    phi_lin = intra_e_sigma * np.ones(vs30.shape)
    phi_lin_pga = intra_e_sigma_pga * np.ones(vs30.shape)

    # Assume that the site response variability is constant with period.
    phi_amp = 0.3
    phi_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)
    if period < 0.01:
        rhoW = 1.0
        rhoB = 1.0
    else:
        rhoW = rho_W_itp(np.log(period))
        rhoB = rho_B_itp(np.log(period))

    f2 = C["f4"] * (
        np.exp(C["f5"] * (np.minimum(vs30, 760.0) - 360.0))
        - np.exp(C["f5"] * (760.0 - 360.0))
    )
    f3 = 0.1

    partial_f_pga = f2 * pgar / (pgar + f3)
    partial_f_pga = partial_f_pga * np.ones(vs30.shape)

    # nonlinear variance components
    phi2_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 [partial_f_pga, np.sqrt(tau2_NL), np.sqrt(phi2_NL)]
    return [np.sqrt(tau2_NL + phi2_NL), np.sqrt(tau2_NL), np.sqrt(phi2_NL)]


[docs]def get_backarc_term(trt, imt, ctx):
    """
    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 backarc sites as
    function of distance.
    """
    period = imt.period
    w_epi_factor = 1.008

    theta7_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.
    if period < 0:
        theta7 = 0.0
        theta8 = 0.0
    elif period >= 0 and period < 0.02:
        theta7 = 1.0988
        theta8 = -1.42
    else:
        theta7 = theta7_itp(np.log(period))
        theta8 = theta8_itp(np.log(period))

    dists = ctx.rrup

    if trt == const.TRT.SUBDUCTION_INTRASLAB:
        min_dist = 85.0
        backarc = np.bool_(ctx.backarc)
        f_faba = np.zeros_like(dists)
        fixed_dists = dists[backarc]
        fixed_dists[fixed_dists < min_dist] = min_dist
        f_faba[backarc] = theta7 + theta8 * np.log(fixed_dists / 40.0)
        return f_faba * w_epi_factor
    else:
        f_faba = np.zeros_like(dists)
        return f_faba


[docs]class Atkinson2022Crust(GMPE):
    """
    Implements Atkinson (2022) backbone model for New Zealand. For
    more info please refere to Gail Atkinson's NSHM report and linked
    revisions.
    """

    DEFINED_FOR_TECTONIC_REGION_TYPE = const.TRT.ACTIVE_SHALLOW_CRUST

    #: Supported intensity measure types are spectral acceleration,
    #: peak ground acceleration
    DEFINED_FOR_INTENSITY_MEASURE_TYPES = {PGA, SA}

    #: Supported intensity measure component is the RotD50
    DEFINED_FOR_INTENSITY_MEASURE_COMPONENT = const.IMC.RotD50

    DEFINED_FOR_STANDARD_DEVIATION_TYPES = {
        const.StdDev.TOTAL,
        const.StdDev.INTER_EVENT,
        const.StdDev.INTRA_EVENT,
    }

    REQUIRES_DISTANCES = {"rrup"}

    REQUIRES_RUPTURE_PARAMETERS = {"mag"}

    REQUIRES_SITES_PARAMETERS = {"vs30"}

    # define the epistemic uncertainities : Central/Lower/Upper

    REQUIRES_ATTRIBUTES = {"epistemic"}

    # define constant parameters
    suffix = "crust"

    def __init__(self, epistemic="Central", modified_sigma=False):
        """
        Aditional parameter for epistemic central,
        lower and upper bounds.
        """
        self.epistemic = epistemic
        self.modified_sigma = modified_sigma

[docs]    def compute(self, ctx: np.recarray, imts, mean, sig, tau, phi):
        trt = self.DEFINED_FOR_TECTONIC_REGION_TYPE
        C_PGA = self.COEFFS[PGA()]
        pga_rock = _get_pga_on_rock(self.suffix, C_PGA, ctx) * np.exp(
            get_backarc_term(trt, PGA(), ctx)
        )
        # Here the backarc term is applied as multiplication because the
        # pga_rock is in linear space not in log space

        for m, imt in enumerate(imts):
            C = self.COEFFS[imt]

            # compute mean
            mean[m] = (
                _fmag(self.suffix, C, ctx.mag)
                + _fz_ha18(C, ctx)
                + _fgamma(self.suffix, C, ctx)
                + fs_SS14(C, pga_rock, ctx)
                + get_backarc_term(trt, imt, ctx)
            )

            mean[m] = mean[m] - np.log(981.0)  # Convert the cm/s^2 to g.

            # In her email and slack post Gail mentioned that her
            # upper and lower branches are as 1.28 times of the delta.
            # So as to represent 10th and 90th percentile.
            # Consequently, the weights have also changed to 0.3, 0.4,
            # 0.3.  The scale factor of 0.9 is applied based upon the
            # discussion that it accounts for the reduction in
            # epistemic uncertainty when no perfect correlation is
            # assumed between rupture scenarios. See the note of Peter
            # and Brendon on slack.
            epistemic_scale_factor = 0.893
            if self.epistemic.lower() == "lower":
                mean[m] = (
                    mean[m]
                    - _epistemic_adjustment_lower(C, ctx)
                    * 1.28155
                    * epistemic_scale_factor
                )

            elif self.epistemic.lower() == "upper":
                mean[m] = (
                    mean[m]
                    + _epistemic_adjustment_upper(C, ctx)
                    * 1.28155
                    * epistemic_scale_factor
                )
            else:
                mean[m] = mean[m]
            # Aleatory Uncertainty terms.
            if self.modified_sigma:
                sig[m], tau[m], phi[m] = get_nonlinear_stddevs(
                    self.suffix, C, C_PGA, imt, pga_rock, ctx.vs30
                )
            else:
                sig[m], tau[m], phi[m] = get_stddevs(self.suffix, C)

    # periods given by 1 / 10 ** COEFFS['f']
    with Atk22_COEFFS.open() as coefs_file:
        COEFFS = CoeffsTable(sa_damping=5, table=coefs_file.read())


[docs]class Atkinson2022SInter(Atkinson2022Crust):
    """
    Atkinson 2022 for Subduction Interface in NZ.
    """

    DEFINED_FOR_TECTONIC_REGION_TYPE = const.TRT.SUBDUCTION_INTERFACE

    # constant table suffix
    suffix = "inter"
    # stress = 100


[docs]class Atkinson2022SSlab(Atkinson2022Crust):
    """
    Atkinson (2022) for Subduction IntraSlab in NZ.
    """

    DEFINED_FOR_TECTONIC_REGION_TYPE = const.TRT.SUBDUCTION_INTRASLAB

    # constant table suffix
    suffix = "slab"
    REQUIRES_SITES_PARAMETERS = {"vs30", "backarc"}