R/i_2014_nga.R

In wltcwpf/GMPE: A pacakge of commonly used Ground Motion Prediction Equations (GMPE)

#' The GMPE of CB 2014
#'
#' This function calculates the ground motion median values and standard deviations
#' @param M Moment magnitude, a numeric value
#' @param Vs30 Shear wave velocity averaged over top 30 m (in m/s).
#' @param T Period (sec); Use 1000 (by default) for output the array of Sa with original NGA West2 periods.
#' Note no PGV computation
#' @param Rrup Closest distance (km) to the ruptured plane
#' @param Fault_Type The indicator of fault type: 0 for strike slip & other non-reverse faulting;
#' 1 for reverse
#' @return A list of five elements is returned: med - median spectral acceleration prediction (in g);
#' sigma - logarithmic standard deviation of spectral acceleration prediction; phi - logarithmic
#' standard deviation of within event residuals; tau - logarithmic standard deviation of between
#' event residuals; period - the corresponding oscillator periods
#' @examples i_2014_nga(M = 5.5, Vs30 = 550, T = 1000, Rrup = 90, Fault_Type = 1)
#' @references Idriss, I. M. (2014). An NGA-West2 Empirical Model for Estimating
#' the Horizontal Spectral Values Generated by Shallow Crustal Earthquakes.
#' Earthquake Spectra, 30(3), 1155-1177.
#' @export
#' @importFrom stats approx
i_2014_nga <- function(M, Vs30, T = 1000, Rrup, Fault_Type) {

  period = c(0.01,	0.02,	0.03,	0.05,	0.075,	0.1,	0.15,	0.2,	0.25,	0.3,	0.4,	0.5,	0.75,
             1,	1.5,	2,	3,	4,	5,	7.5,	10, 0)

  if (Vs30 > 1200)
    Vs30 = 1200

  if (length(T) == 1 & T[1] == 1000) {

    # use pre-defined period
    Sa <- rep(0, length(period))

    Sigma <- rep(0, length(period))

    Phi <- rep(0, length(period))

    Tau <- rep(0, length(period))

    for (ip in 1:length(period)) {

      res <- i_2014_subroutine(M, ip, Rrup, Vs30, Fault_Type)

      Sa[ip] <- res[1]

      Sigma[ip] <- res[2]

      Phi[ip] <- res[3]

      Tau[ip] <- res[4]

    }

    T <- period

  } else {

    Sa <- rep(0, length(T))

    Sigma <- rep(0, length(T))

    Phi <- rep(0, length(T))

    Tau <- rep(0, length(T))

    for (i in 1:length(T)) {

      Ti <- T[i]

      if (min(abs(period - Ti)) > 0.0001) {

        # user defined period needs interpolation
        T_low <- max(period[period < Ti])

        T_high <- min(period[period > Ti])

        ip_low <- which(period == T_low)

        ip_high <- which(period == T_high)

        res_low <- i_2014_subroutine(M, ip_low, Rrup, Vs30, Fault_Type)

        res_high <- i_2014_subroutine(M, ip_high, Rrup, Vs30, Fault_Type)

        Sa[i] = approx(x = c(T_low, T_high), y = c(res_low[1], res_high[1]),
                           xout = Ti, rule = 2)$y

        Sigma[i] = approx(x = c(T_low, T_high), y = c(res_low[2], res_high[2]),
                              xout = Ti, rule = 2)$y

        Phi[i] = approx(x = c(T_low, T_high), y = c(res_low[3], res_high[3]),
                            xout = Ti, rule = 2)$y

        Tau[i] = approx(x = c(T_low, T_high), y = c(res_low[4], res_high[4]),
                            xout = Ti, rule = 2)$y

      } else {

        ip_T = which(abs(period- Ti) < 0.0001)

        res = i_2014_subroutine(M, ip_T, Rrup, Vs30, Fault_Type)

        Sa[i] <- res[1]

        Sigma[i] <- res[2]

        Phi[i] <- res[3]

        Tau[i] <- res[4]

      }
    }
  }

  res <- list()

  res$med <- Sa

  res$sigma <- Sigma

  res$phi <- Phi

  res$tau <- Tau

  res$period <- T

  return(res)
}

#' The subroutine of GMPE of Idriss 2014
#'
#' This is a subroutine of Idriss 2014
#' @param M Moment magnitude, a numeric value
#' @param ip The index of working period in the per-defined periods: (0.01,	0.02,	0.03,	0.05,	0.075,	0.1,
#' 0.15,	0.2,	0.25,	0.3,	0.4,	0.5,	0.75, 1,	1.5,	2,	3,	4,	5,	7.5,	10, 0)
#' @param R_RUP Closest distance (km) to the ruptured plane
#' @param Vs30 Shear wave velocity averaged over top 30 m (in m/s).
#' @param Fault_Type The indicator of fault type: 0 for strike slip & other non-reverse faulting;
#' 1 for reverse
#' @return An array of four values: (1) median spectral acceleration prediction (in g); (2)
#' logarithmic standard deviation of spectral acceleration prediction; (3) logarithmic
#' standard deviation of within event residuals; and (4) logarithmic standard deviation of between
#' event residuals at the period of ip.
#' @references Idriss, I. M. (2014). An NGA-West2 Empirical Model for Estimating
#' the Horizontal Spectral Values Generated by Shallow Crustal Earthquakes.
#' Earthquake Spectra, 30(3), 1155-1177.
#' @export
#' @importFrom stats approx
i_2014_subroutine <- function(M, ip, R_RUP, Vs30, Fault_Type) {

  period = c(0.01,	0.02,	0.03,	0.05,	0.075,	0.1,	0.15,	0.2,	0.25,	0.3,	0.4,	0.5,	0.75,
             1,	1.5,	2,	3,	4,	5,	7.5,	10,	0)

  if (M <= 6.75) {

    a1 = c(7.0887,	7.1157,	7.2087,	6.2638,	5.9051,	7.5791,	8.0190,	9.2812,	9.5804,	9.8912,	9.5342,
           9.2142,	8.3517,	7.0453,	5.1307,	3.3610,	0.1784,	-2.4301,	-4.3570,	-7.8275,	-9.2857,	7.0887)

    a2 = c(0.2058,	0.2058,	0.2058,	0.0625,	0.1128,	0.0848,	0.1713,	0.1041,	0.0875,	0.0003,	0.0027,	0.0399,
           0.0689,	0.1600,	0.2429,	0.3966,	0.7560,	0.9283,	1.1209,	1.4016,	1.5574,	0.2058)

    b1 = c(2.9935,	2.9935,	2.9935,	2.8664,	2.9406,	3.0190,	2.7871,	2.8611,	2.8289,	2.8423,	2.8300,	2.8560,
           2.7544,	2.7339,	2.6800,	2.6837,	2.6907,	2.5782,	2.5468,	2.4478,	2.3922,	2.9935)

    b2 = c(-0.2287,	-0.2287,	-0.2287,	-0.2418,	-0.2513,	-0.2516,	-0.2236,	-0.2229,	-0.2200,	-0.2284,
           -0.2318,	-0.2337,	-0.2392,	-0.2398,	-0.2417,	-0.2450,	-0.2389,	-0.2514,	-0.2541,	-0.2593,
           -0.2586,	-0.2287)

  } else {

    a1 = c(9.0138,	9.0408,	9.1338,	7.9837,	7.7560,	9.4252,	9.6242,	11.1300,	11.3629,	11.7818,	11.6097,
           11.4484,	10.9065,	9.8565,	8.3363,	6.8656,	4.1178,	1.8102,	0.0977,	-3.0563,	-4.4387,	9.0138)

    a2 = c(-0.0794,	-0.0794,	-0.0794,	-0.1923,	-0.1614,	-0.1887,	-0.0665,	-0.1698,	-0.1766,
           -0.2798,	-0.3048,	-0.2911,	-0.3097,	-0.2565,	-0.2320,	-0.1226,	0.1724,	0.3001,	0.4609,
           0.6948,	0.8393,	-0.0794)

    b1 = c(2.9935,	2.9935,	2.9935,	2.7995,	2.8143,	2.8131,	2.4091,	2.4938,	2.3773,	2.3772,	2.3413,	2.3477,
           2.2042,	2.1493,	2.0408,	2.0013,	1.9408,	1.7763,	1.7030,	1.5212,	1.4195,	2.9935)

    b2 = c(-0.2287,	-0.2287,	-0.2287,	-0.2319,	-0.2326,	-0.2211,	-0.1676,	-0.1685,	-0.1531,	-0.1595,
           -0.1594,	-0.1584,	-0.1577,	-0.1532,	-0.1470,	-0.1439,	-0.1278,	-0.1326,	-0.1291,	-0.1220,
           -0.1145,	-0.2287)

  }

  a3 = c(0.0589,	0.0589,	0.0589,	0.0417,	0.0527,	0.0442,	0.0329,	0.0188,	0.0095,	-0.0039,	-0.0133,
         -0.0224,	-0.0267,	-0.0198,	-0.0367,	-0.0291,	-0.0214,	-0.0240,	-0.0202,	-0.0219,
         -0.0035,	0.0589)

  zeta = c(-0.854,	-0.854,	-0.854,	-0.631,	-0.591,	-0.757,	-0.911,	-0.998,	-1.042,	-1.030,	-1.019,	-1.023,
           -1.056,	-1.009,	-0.898,	-0.851,	-0.761,	-0.675,	-0.629,	-0.531,	-0.586,	-0.854)

  gamma	= c(-0.0027,	-0.0027,	-0.0027,	-0.0061,	-0.0056,	-0.0042,	-0.0046,	-0.0030,	-0.0028,
            -0.0029,	-0.0028,	-0.0021,	-0.0029,	-0.0032,	-0.0033,	-0.0032,	-0.0031,	-0.0051,
            -0.0059,	-0.0057,	-0.0061,	-0.0027)

  phi	= c(0.08,	0.08,	0.08,	0.08,	0.08,	0.08,	0.08,	0.08,	0.08,	0.08,	0.08,	0.08,	0.08,	0.06,	0.04,	0.02,
          0.02,	0,	0,	0,	0,	0.08)

  ln_PSA = a1[ip] + a2[ip] * M + a3[ip] * (8.5 - M)^2 - (b1[ip] + b2[ip] * M) * log(R_RUP + 10) +
    zeta[ip] * log(Vs30) + gamma[ip] * R_RUP + phi[ip] * Fault_Type

  Sa = exp(ln_PSA)

  M_SE = max(min(M, 7.5), 5)

  if (period[ip] < 0.05) {

    Sigma = 1.18 + 0.035 * log(0.05) - 0.06 * M_SE

  } else if (period[ip] > 3) {

    Sigma = 1.18 + 0.035 * log(3) - 0.06 * M_SE

  } else {

    Sigma = 1.18 + 0.035 * log(period[ip]) - 0.06 * M_SE

  }

  tau <- sqrt(Sigma^2 - phi[ip]^2)

  return(c(Sa, Sigma, phi[ip], tau))

}