Source code for openquake.hazardlib.gsim.tusa_langer_2016

# -*- 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:`TusaLanger2016RepiBA08SE`,
               :class:`TusaLanger2016RepiBA08DE`,
               :class:`TusaLanger2016RepiSP87SE`,
               :class:`TusaLanger2016RepiSP87DE`,
               :class:`TusaLanger2016Rhypo`
"""
import numpy as np
from scipy.constants import g

from openquake.hazardlib.gsim.base import GMPE, CoeffsTable
from openquake.hazardlib import const
from openquake.hazardlib.imt import PGA, SA
from openquake.baselib.general import CallableDict

_compute_distance = CallableDict()


@_compute_distance.add("BA08SE", "BA08DE")
def _compute_distance1(kind, ctx, C):
    """
    Compute the distance function, equation (9):
    """
    mref = 3.6
    rref = 1.0
    rval = np.sqrt(ctx.repi ** 2 + C['h'] ** 2)
    return (C['c1'] + C['c2'] * (ctx.mag - mref)) * np.log10(
        rval / rref) + C['c3'] * (rval - rref)


@_compute_distance.add("SP87SE", "SP87DE")
def _compute_distance2(kind, ctx, C):
    """
    Compute the distance function, equation (5).
    """
    rval = np.sqrt(ctx.repi ** 2 + C['h'] ** 2)
    return C['c1'] * np.log10(rval)


@_compute_distance.add("Rhypo")
def _compute_distance3(kind, ctx, C):
    """
    Compute the distance function, equation (9):
    """
    mref = 3.6
    rref = 1.0
    rval = np.sqrt(ctx.rhypo ** 2 + C['h'] ** 2)
    return (C['c1'] + C['c2'] * (ctx.mag - mref)) * np.log10(
        rval / rref) + C['c3'] * (rval - rref)


_compute_magnitude = CallableDict()


@_compute_magnitude.add("BA08SE", "BA08DE", "Rhypo")
def _compute_magnitude1(kind, ctx, C):
    """
    Compute the magnitude function, equation (9):
    """
    return C['a'] + C['b1'] * ctx.mag + C['b2'] * ctx.mag ** 2


@_compute_magnitude.add("SP87SE", "SP87DE")
def _compute_magnitude2(kind, ctx, C):
    """
    Compute the magnitude function, equation (5).
    """
    return C['a'] + C['b1'] * ctx.mag


def _get_site_amplification(ctx, C):
    """
    Compute the site amplification function given by FS = eiSi, for
    i = 1,2,3 where Si are the coefficients determined through regression
    analysis, and ei are dummy variables (0 or 1) used to denote the
    different EC8 site classes.
    """
    ssa, ssb, ssd = _get_site_type_dummy_variables(ctx)

    return C['sA'] * ssa + C['sB'] * ssb + C['sD'] * ssd


def _get_site_type_dummy_variables(ctx):
    """
    Get site type dummy variables, which classified the ctx into
    different site classes based on the shear wave velocity in the
    upper 30 m (Vs30) according to the EC8 (CEN 2003):
    class A: Vs30 > 800 m/s
    class B: Vs30 = 360 - 800 m/s
    class C*: Vs30 = 180 - 360 m/s
    class D: Vs30 < 180 m/s
    *Not computed by this GMPE
    """
    ssa = np.zeros(len(ctx.vs30))
    ssb = np.zeros(len(ctx.vs30))
    ssd = np.zeros(len(ctx.vs30))
    # Class D; Vs30 < 180 m/s.
    ssd[ctx.vs30 < 180.0] = 1.0
    # Class B; 360 m/s <= Vs30 <= 800 m/s.
    ssb[(ctx.vs30 >= 360.0) & (ctx.vs30 < 800.0)] = 1.0
    # Class A; Vs30 > 800 m/s.
    ssa[ctx.vs30 >= 800.0] = 1.0
    if ((ctx.vs30 > 180) & (ctx.vs30 < 360)).any():
        raise Exception(
            'GMPE does not consider site class C (Vs30 = 180-360 m/s)')

    return ssa, ssb, ssd


[docs]class TusaLanger2016RepiBA08SE(GMPE):
    """
    Implements GMPE developed by Giuseppina Tusa and Horst Langer (2016) and
    published as "Prediction of ground motion parameters for the volcanic area
    of Mount Etna" Journal of Seismology, DOI 10.1007/s10950-015-9508-x.

    GMPE derives from earthquakes in the volcanic area of Mt. Etna in the
    magnitude range 3<ML<4.8 for epicentral distances <100 km, and for soil
    classes A, B, and D. Authors do NOT derive coefficients for site class C
    due to limited data. For implementation using hypocentral distance see
    :class:`TusaLanger2016Rhypo`.

    Two functional forms were considered by the authors: Sabetta and Pugliese,
    1987 (SP87) and a simplified version of Boore and Atkinson, 2008 (BA08).
    The GMPE distinguishes between shallow volcano-tectonic events related to
    flank movements (SE, focal depths <5km) and deeper events occurring due to
    regional tectonics (DE, focal depths >5km).

    Test tables are generated from a spreadsheet provided by the authors, and
    modified according to OQ format (e.g. conversion from cm/s2 to m/s2).

    Jan 2019: After noticing an anomalous-looking spike in the response spectra
    of the TusaLanger2016RepiBA08SE model at T=0.14s, we contacted the authors
    who found a mistake in one of the coefficients in the publication. It has
    been updated according to the authors suggestion.
    """
    kind = "BA08SE"

    #: Supported tectonic region type is 'volcanic' because the
    #: equations have been derived from data from Etna (Sicily, Italy)
    DEFINED_FOR_TECTONIC_REGION_TYPE = const.TRT.VOLCANIC

    #: Supported intensity measure types are PGA and SA
    DEFINED_FOR_INTENSITY_MEASURE_TYPES = {PGA, SA}

    #: Supported intensity measure component is the maximum of two horizontal
    #: components
    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}

    #: Required site parameter is Vs30
    REQUIRES_SITES_PARAMETERS = {'vs30'}

    #: Required rupture parameter is magnitude.
    REQUIRES_RUPTURE_PARAMETERS = {'mag'}

    #: Required distance measure is Repi
    REQUIRES_DISTANCES = {'repi'}

[docs]    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.
        """
        for m, imt in enumerate(imts):
            C = self.COEFFS[imt]

            imean = (_compute_magnitude(self.kind, ctx, C) +
                     _compute_distance(self.kind, ctx, C) +
                     _get_site_amplification(ctx, C))

            # convert from log10 to ln and from cm/s**2 to g
            mean[m] = np.log((10.0 ** (imean - 2.0)) / g)
            # Return stddevs in terms of natural log scaling
            sig[m] = np.log(10.0 ** C['SigmaTot'])

    # Coefficients from Table 9 (PGA) and Table 12 (SA); sigma values in log
    # Correction made to coeff b2 (0.14s) by Giuseppina Tusa in email (Jan 15, 2019)
    COEFFS = CoeffsTable(sa_damping=5, table="""
    IMT    a       b1      b2      c1      c2       h       c3      sA  sB      sD      SigmaIE SigmaIS SigmaTot
    pga   -0.568   0.475   0.037   -2.054  -0.015   2.317   0.006   0   0.452   0.477   0.222   0.221   0.390
    10.0  -4.794   0.827   0.054   -1.273   0.096   2.221  -0.002   0   0.363   0.340   0.164   0.243   0.351
    5.00  -4.010   0.646   0.083   -1.036   0.142   1.517  -0.003   0   0.296   0.241   0.163   0.263   0.363
    3.33  -5.200   1.456  -0.015   -0.953   0.092   1.728  -0.006   0   0.281   0.246   0.168   0.263   0.369
    2.50  -5.753   1.912  -0.069   -1.024   0.025   2.641  -0.006   0   0.343   0.353   0.158   0.255   0.359
    2.00  -5.531   2.016  -0.089   -1.139   0.037   3.044  -0.005   0   0.375   0.367   0.157   0.250   0.355
    1.67  -3.923   1.431  -0.023   -1.369   0.061   3.542  -0.003   0   0.403   0.384   0.150   0.238   0.342
    1.43  -3.023   1.184  -0.011   -1.461   0.114   3.668  -0.001   0   0.438   0.385   0.142   0.243   0.342
    1.25  -2.558   1.091  -0.014   -1.468   0.142   3.565  -0.001   0   0.433   0.379   0.131   0.246   0.341
    1.11  -1.921   0.847   0.011   -1.523   0.157   3.580   0.000   0   0.450   0.397   0.126   0.256   0.346
    1.00  -1.358   0.629   0.030   -1.522   0.191   3.356   0.000   0   0.465   0.416   0.124   0.263   0.353
    0.67  -1.348   0.824  -0.012   -1.631   0.152   3.229   0.003   0   0.532   0.510   0.128   0.247   0.339
    0.50  -2.278   1.574  -0.115   -1.897   0.075   3.761   0.005   0   0.486   0.523   0.144   0.234   0.337
    0.40  -2.601   1.771  -0.141   -1.916  -0.001   3.513   0.004   0   0.491   0.559   0.158   0.219   0.332
    0.33  -1.309   1.112  -0.051   -1.957  -0.020   3.442   0.004   0   0.460   0.510   0.178   0.212   0.339
    0.29   0.025   0.441   0.047   -2.134  -0.037   3.759   0.006   0   0.455   0.505   0.201   0.215   0.356
    0.25   0.127   0.413   0.051   -2.219  -0.051   3.695   0.007   0   0.468   0.528   0.220   0.221   0.374
    0.22   0.139   0.448   0.040   -2.259  -0.038   3.476   0.007   0   0.478   0.546   0.233   0.222   0.385
    0.20   0.136   0.455   0.038   -2.266  -0.044   3.363   0.007   0   0.468   0.576   0.240   0.223   0.392
    0.17   0.021   0.556   0.016   -2.240  -0.026   2.997   0.006   0   0.423   0.542   0.263   0.234   0.413
    0.14  -0.289   0.741  -0.013   -2.240  -0.023   2.676   0.006   0   0.403   0.516   0.277   0.236   0.428
    0.13  -0.119   0.593   0.011   -2.243  -0.035   2.402   0.007   0   0.412   0.490   0.285   0.236   0.437
    0.11  -0.015   0.476   0.031   -2.219  -0.044   2.316   0.006   0   0.411   0.456   0.288   0.233   0.439
    0.10  -0.221   0.522   0.030   -2.193  -0.062   2.300   0.006   0   0.412   0.435   0.283   0.231   0.437
    """)


[docs]class TusaLanger2016RepiBA08DE(TusaLanger2016RepiBA08SE):
    """
    Implements Tusa and Langer (2016) using the BA08 model and DE.

    Extends
    :class:`openquake.hazardlib.gsim.tusa_langer_2016.TusaLanger2016RepiBA08SE`
    because the same functional form is used, only the coefficients differ.
    """
    kind = "BA08DE"

    # Coefficients from Table 9 (PGA) and Table 13 (SA); sigma values in log
    COEFFS = CoeffsTable(sa_damping=5, table="""
    IMT    a       b1      b2       c1      c2      h       c3      sA  sB      sD      SigmaIE SigmaIS SigmaTot
    pga    2.527  -0.397   0.094   -1.998   0.315   9.6080  0.000   0  -0.200   0.002   0.161   0.276   0.399
    10.0  -3.019  -0.236   0.123   -0.716   0.144   2.0850 -0.004   0   0.115   0.006   0.170   0.203   0.342
    5.00  -2.383  -0.325   0.147   -0.695   0.062   2.0800 -0.002   0   0.120   0.003   0.167   0.213   0.347
    3.33  -2.170  -0.307   0.155   -0.633   0.002   2.2200 -0.002   0   0.107  -0.023   0.162   0.226   0.359
    2.50  -2.038  -0.210   0.136   -0.620   0.053   2.1350 -0.003   0   0.132  -0.026   0.169   0.215   0.348
    2.00  -1.977  -0.056   0.106   -0.625   0.120   2.3880 -0.004   0   0.139  -0.041   0.179   0.207   0.343
    1.67  -1.893   0.048   0.088   -0.635   0.158   2.8030 -0.005   0   0.128  -0.040   0.179   0.206   0.340
    1.43  -1.656   0.065   0.078   -0.667   0.199   3.2940 -0.005   0   0.136  -0.022   0.186   0.208   0.345
    1.25  -1.441   0.059   0.075   -0.702   0.216   4.0100 -0.005   0   0.145   0.004   0.189   0.206   0.344
    1.11  -1.054   0.025   0.076   -0.834   0.229   5.5460 -0.004   0   0.140  -0.003   0.193   0.203   0.346
    1.00  -0.945   0.041   0.074   -0.871   0.226   5.5910 -0.004   0   0.140  -0.007   0.191   0.201   0.345
    0.67   0.202  -0.354   0.121   -0.902   0.244   3.9200 -0.003   0   0.151   0.050   0.183   0.197   0.340
    0.50   1.322  -0.596   0.146   -1.135   0.256   5.9340 -0.002   0   0.062   0.037   0.181   0.201   0.342
    0.40   2.073  -0.688   0.152   -1.420   0.274   7.0510  0.000   0   0.021   0.024   0.186   0.207   0.351
    0.33   2.131  -0.635   0.144   -1.484   0.277   6.8740  0.000   0   0.037   0.035   0.191   0.223   0.362
    0.29   2.496  -0.600   0.138   -1.761   0.282   8.2450  0.002   0   0.046   0.050   0.195   0.238   0.371
    0.25   2.651  -0.482   0.116   -1.953   0.305   8.8410  0.003   0   0.021   0.051   0.189   0.252   0.380
    0.22   2.967  -0.358   0.097   -2.294   0.296   10.742  0.006   0  -0.014   0.058   0.183   0.260   0.383
    0.20   3.633  -0.334   0.091   -2.764   0.299   12.651  0.009   0  -0.061   0.666   0.177   0.269   0.389
    0.17   3.837  -0.311   0.085   -2.820   0.301   13.203  0.007   0  -0.145   0.079   0.178   0.282   0.403
    0.14   3.980  -0.328   0.081   -2.811   0.316   13.357  0.006   0  -0.197   0.051   0.177   0.292   0.414
    0.13   4.139  -0.335   0.077   -2.879   0.337   13.486  0.006   0  -0.232   0.015   0.183   0.306   0.427
    0.11   3.918  -0.339   0.077   -2.736   0.342   12.535  0.005   0  -0.241  -0.025   0.186   0.312   0.434
    0.10   3.746  -0.339   0.074   -2.612   0.354   12.019  0.004   0  -0.247  -0.054   0.192   0.317   0.443
    """)


[docs]class TusaLanger2016RepiSP87SE(TusaLanger2016RepiBA08SE):
    """
    Implements Tusa and Langer (2016) using the SP87 model and SE.

    Extends
    :class:`openquake.hazardlib.gsim.tusa_langer_2016.TusaLanger2016RepiBA08SE`
    with modification to the functional form and different coefficients.
    """
    kind = "SP87SE"

    # Coefficients from Table 8 (PGA) and Table 10 (SA); sigma values in log
    COEFFS = CoeffsTable(sa_damping=5, table="""
    IMT    a       b1      c1      h       sA  sB      sD      SigmaIE SigmaIS SigmaTot
    pga   -1.186   0.726  -1.719   1.551   0   0.357   0.376   0.223   0.229   0.393
    10.0  -5.799   1.335  -1.408   2.731   0   0.379   0.354   0.165   0.243   0.351
    5.00  -5.528   1.424  -1.274   2.510   0   0.336   0.280   0.165   0.263   0.364
    3.33  -4.975   1.464  -1.415   4.285   0   0.312   0.268   0.167   0.265   0.371
    2.50  -4.410   1.454  -1.568   5.747   0   0.354   0.348   0.156   0.258   0.360
    2.00  -4.116   1.435  -1.589   5.364   0   0.387   0.368   0.155   0.252   0.357
    1.67  -3.676   1.352  -1.611   4.623   0   0.410   0.386   0.149   0.238   0.343
    1.43  -3.328   1.265  -1.573   4.078   0   0.438   0.383   0.143   0.243   0.342
    1.25  -2.978   1.189  -1.568   3.866   0   0.432   0.378   0.132   0.247   0.341
    1.11  -2.808   1.143  -1.536   3.520   0   0.441   0.390   0.127   0.257   0.347
    1.00  -2.655   1.108  -1.530   3.248   0   0.454   0.406   0.124   0.264   0.354
    0.67  -2.077   0.954  -1.468   2.577   0   0.491   0.471   0.128   0.249   0.340
    0.50  -1.474   0.866  -1.595   2.795   0   0.426   0.468   0.145   0.238   0.339
    0.40  -1.113   0.778  -1.613   2.585   0   0.433   0.506   0.160   0.223   0.334
    0.33  -0.839   0.729  -1.662   2.575   0   0.410   0.463   0.179   0.216   0.340
    0.29  -0.739   0.725  -1.723   2.615   0   0.395   0.448   0.202   0.221   0.358
    0.25  -0.661   0.710  -1.757   2.484   0   0.396   0.458   0.220   0.229   0.377
    0.22  -0.574   0.688  -1.777   2.297   0   0.389   0.458   0.233   0.232   0.388
    0.20  -0.501   0.669  -1.795   2.247   0   0.376   0.484   0.241   0.233   0.395
    0.17  -0.357   0.643  -1.844   2.125   0   0.333   0.450   0.263   0.242   0.415
    0.14  -0.275   0.623  -1.873   1.921   0   0.307   0.417   0.277   0.244   0.430
    0.13  -0.334   0.627  -1.863   1.669   0   0.304   0.440   0.286   0.245   0.439
    0.11  -0.427   0.640  -1.843   1.593   0   0.303   0.341   0.289   0.242   0.442
    0.10  -0.538   0.653  -1.814   1.556   0   0.306   0.322   0.285   0.241   0.440
    """)


[docs]class TusaLanger2016RepiSP87DE(TusaLanger2016RepiSP87SE):
    """
    Implements Tusa and Langer (2016) using the SP87 model and DE.

    Extends
    :class:`openquake.hazardlib.gsim.tusa_langer_2016.TusaLanger2016RepiSP87SE`
    because the same functional form is used, only the coefficients differ.
    """
    kind = 'SP87DE'

    # Coefficients from Table 8 (PGA) and Table 11 (SA); sigma values in log
    COEFFS = CoeffsTable(sa_damping=5, table="""
    IMT    a       b1      c1      h        sA  sB      sD      SigmaIE SigmaIS SigmaTot
    pga   -0.377   0.765  -1.824   9.5270   0  -0.202   0.004   0.162   0.277   0.402
    10.0  -4.942   0.877  -1.035   6.1470   0   0.122   0.013   0.172   0.204   0.344
    5.00  -4.444   0.853  -0.879   4.4510   0   0.133   0.017   0.170   0.213   0.349
    3.33  -4.039   0.841  -0.841   4.6740   0   0.120  -0.011   0.165   0.227   0.361
    2.50  -3.804   0.871  -0.886   5.6350   0   0.141  -0.019   0.172   0.216   0.350
    2.00  -3.470   0.899  -1.015   7.4860   0   0.138  -0.044   0.181   0.208   0.345
    1.67  -3.164   0.923  -1.134   8.8600   0   0.123  -0.048   0.181   0.207   0.342
    1.43  -3.007   0.926  -1.140   9.0160   0   0.132  -0.029   0.187   0.209   0.348
    1.25  -2.838   0.928  -1.169   9.5550   0   0.138  -0.006   0.191   0.207   0.348
    1.11  -2.642   0.920  -1.200   9.7200   0   0.132  -0.013   0.194   0.204   0.348
    1.00  -2.555   0.913  -1.182   9.2180   0   0.134  -0.015   0.192   0.201   0.348
    0.67  -2.242   0.895  -1.130   6.9400   0   0.152   0.050   0.185   0.197   0.344
    0.50  -1.567   0.858  -1.305   8.1620   0   0.056   0.031   0.185   0.201   0.346
    0.40  -1.268   0.840  -1.359   7.5770   0   0.018   0.023   0.190   0.207   0.355
    0.33  -1.201   0.839  -1.349   6.9760   0   0.032   0.034   0.194   0.223   0.366
    0.29  -1.072   0.833  -1.385   6.9500   0   0.045   0.055   0.197   0.239   0.375
    0.25  -0.908   0.825  -1.442   6.8940   0   0.023   0.061   0.190   0.253   0.384
    0.22  -0.697   0.806  -1.505   7.4290   0   0.001   0.083   0.183   0.262   0.387
    0.20  -0.535   0.794  -1.558   7.6410   0  -0.031   0.109   0.177   0.271   0.393
    0.17  -0.097   0.777  -1.763   9.0800   0  -0.116   0.119   0.177   0.284   0.406
    0.14   0.215   0.758  -1.907   10.038   0  -0.172   0.085   0.176   0.294   0.416
    0.13   0.215   0.758  -1.907   10.038   0  -0.172   0.085   0.182   0.307   0.430
    0.11   0.239   0.747  -1.930   9.6770   0  -0.221   0.004   0.184   0.313   0.437
    0.10   0.219   0.745  -1.929   9.6990   0  -0.233  -0.031   0.191   0.318   0.446
    """)


[docs]class TusaLanger2016Rhypo(TusaLanger2016RepiBA08SE):
    """
    Implements the GMPE using the BA08 model and hypocentral distance (not
    described in Tusa and Langer, 2016). This version has been developed in the
    frame of V3-2012 INGV-DPC Project in order to perform PSHA calculations
    when topography is taken into consideration (e.g. the flanks of Mt Etna),
    hence dependence on vertical distance is required.

    Extends
    :class:`openquake.hazardlib.gsim.tusa_langer_2016.TusaLanger2016RepiBA08SE`
    because the same functional form is used, only the distance type and
    coefficients differ.
    """
    kind = 'Rhypo'

    # Required distance measure is Rhypo
    REQUIRES_DISTANCES = {'rhypo'}

    # Coefficients provided by Giuseppina Tusa in excel file
    # 'SpectralAccXLaura.xlsx' (email dated May 29, 2016) with modification
    # to pga coefficients in Catania (personal communication, July, 2016);
    # sigma values in log
    COEFFS = CoeffsTable(sa_damping=5, table="""
    IMT    a        b1       b2       c1       c2       h        c3          sA  sB        sD        SigmaTot
    pga    0.3290   0.1050   0.0760  -2.1110   0.0390   1.5530   0.0060000   0   0.45000   0.45700   0.3940
    0.10   0.8594   0.0525   0.0790  -2.2258   0.0066   1.8689   0.0068670   0   0.41397   0.42120   0.4424
    0.20   1.0619   0.0418   0.0804  -2.2684   0.0194   2.6969   0.0072450   0   0.47176   0.56672   0.3960
    0.25   0.9928   0.0316   0.0904  -2.2318   0.0067   3.1724   0.0067460   0   0.47069   0.51751   0.3779
    0.40  -1.8031   1.4272  -0.1051  -1.9414   0.0496   3.0608   0.0047120   0   0.49105   0.54594   0.3353
    0.50  -1.4907   1.2390  -0.0806  -1.9292   0.1235   3.3918   0.0049110   0   0.48505   0.50924   0.3401
    1.00  -0.6283   0.3082   0.0638  -1.5327   0.2391   2.7319   0.0005733   0   0.46483   0.40585   0.3549
    1.25  -1.8557   0.7888   0.0174  -1.4873   0.1883   3.0522  -0.0006545   0   0.43149   0.36715   0.3426
    2.00  -4.8586   1.7503  -0.0610  -1.1999   0.0767   2.8470  -0.0043670   0   0.36832   0.34906   0.3567
    """)