incomplete G18 implementation

src/main/java/gov/usgs/earthquake/nshmp/gmm/Grazier_2018.java

+/**
+ *
+ */
+package gov.usgs.earthquake.nshmp.gmm;
+
+import static gov.usgs.earthquake.nshmp.gmm.FaultStyle.REVERSE;
+import static gov.usgs.earthquake.nshmp.gmm.FaultStyle.REVERSE_OBLIQUE;
+import static gov.usgs.earthquake.nshmp.gmm.GmmInput.Field.MW;
+import static gov.usgs.earthquake.nshmp.gmm.GmmInput.Field.RAKE;
+import static gov.usgs.earthquake.nshmp.gmm.GmmInput.Field.RRUP;
+import static gov.usgs.earthquake.nshmp.gmm.GmmInput.Field.VS30;
+import static java.lang.Math.abs;
+import static java.lang.Math.atan;
+import static java.lang.Math.cos;
+import static java.lang.Math.exp;
+import static java.lang.Math.log;
+import static java.lang.Math.max;
+import static java.lang.Math.min;
+import static java.lang.Math.pow;
+import static java.lang.Math.sqrt;
+
+import java.util.Map;
+
+import com.google.common.collect.Range;
+
+import gov.usgs.earthquake.nshmp.Faults;
+import gov.usgs.earthquake.nshmp.gmm.GmmInput.Constraints;
+import gov.usgs.earthquake.nshmp.tree.LogicTree;
+
+/**
+ *
+ * TODO:
+ *
+ *
+ * * is Q_SA, for G3, defined by Eq. 8?
+ *
+ * * coefficients c21 and c22 are not defined in G18 or GK16 papers.
+ * coefficients b1, b2, b3, and λ are not defined in G18 or GK16 papers,
+ * inferred from matlab script
+ *
+ * * Do we want basin and non-basin flavors?
+ *
+ *
+ * Implementation of the GK17 ground motion model for shallow crustal
+ * continental earthquakes by Grazier (2018). This model computes the average
+ * (RotD50) horizontal components of PGA and 5% damped PSA at periods from 0.01
+ * to 10 s. The model includes a number of significant changes from the previous
+ * GK15 (Graizer and Kalkan, 2016).
+ *
+ * <p><b>References:</b><ul>
+ *
+ * <li>Graizer, V., 2018, GK17 Ground-Motion prediction equation for horizontal
+ * PGA and 5% damped PSA from shallow crustal continental earthquakes: Bulletin
+ * of the Seismological Society of America, v. 108, n. 1, p. 380-398.
+ * <b>doi:</b> <a href="http://dx.doi.org/10.1785/0120170158">
+ * 10.1785/0120170158</a></li>
+ *
+ * <li>Graizer, V., and E. Kalkan, 2016, Summary of the GK15 ground-motion
+ * prediction equation for horizontal PGA and 5% damped PSA from shallow crustal
+ * continental earthquakes: Bulletin of the Seismological Society of America, v.
+ * 106, n. 2, p. 687-707. <b>doi:</b> <a
+ * href="http://dx.doi.org/10.1785/0120150194"> 10.1785/0120150194</a></li>
+ *
+ * </ul>
+ *
+ * <p><b>Component:</b> RotD50 average horizontal components
+ *
+ * @author U.S. Geological Survey
+ * @see Gmm#GK_17
+ */
+public final class Grazier_2018 implements GroundMotionModel {
+
+  public static void main(String[] args) {
+    for (Imt imt : Imt.mprsImts()) {
+      double t = (imt.equals(Imt.PGA)) ? 0.01 : imt.period();
+      System.out.println(t);
+    }
+
+  }
+
+  static final String NAME = "Graizer 2018 (GK17)";
+
+  static final CoefficientContainer COEFFS = new CoefficientContainer("GK15.csv");
+
+  static final Constraints CONSTRAINTS = Constraints.builder()
+      .set(MW, Range.closed(4.0, 8.5))
+      .set(RRUP, Range.closed(0.0, 400.0))
+      .set(RAKE, Faults.RAKE_RANGE)
+      .set(VS30, Range.closed(150.0, 1500.0))
+      .build();
+
+  private static final class Coefficients {
+    final Imt imt;
+    final double σ;
+    final double τ;
+    final double φ;
+
+    Coefficients(Imt imt, CoefficientContainer cc) {
+      this.imt = imt;
+      Map<String, Double> coeffs = cc.get(imt);
+
+      this.σ = coeffs.get("sigma");
+      this.τ = coeffs.get("tau");
+      this.φ = coeffs.get("phi");
+    }
+  }
+
+  private final Imt imt;
+  private final double period;
+  private final Coefficients coeffs;
+
+  public Grazier_2018(Imt imt) {
+    this.imt = imt;
+    this.period = imt.equals(Imt.PGA) ? 0.01 : imt.period();
+
+    coeffs = new Coefficients(imt, COEFFS);
+  }
+
+  /**
+   * G18 provides table of sigmas at 16 periods that do not correspond exactly
+   * with USGS MPRS. The defined sigma values are interpolated in Matlab using
+   * the "shape-preserving piecewise cubic interpolation", and the resulting
+   * values are hard-coded here.
+   *
+   * The following Matlab code was used to compute sigma values at the MSRS
+   * periods:
+   *
+   * <pre>
+   * % From Table 3 in Graizer (2018)
+   * t     = [0.010 0.020 0.040 0.075 0.100 0.150 0.200 0.400 0.500 0.750 1.000 2.000 3.000 4.000 5.000 8.000 10.0 ];
+   * sigma = [0.631 0.633 0.653 0.713 0.734 0.720 0.693 0.673 0.693 0.732 0.738 0.783 0.774 0.733 0.746 0.789 0.736];
+   * tau   = [0.375 0.375 0.386 0.456 0.477 0.451 0.407 0.365 0.398 0.465 0.492 0.501 0.512 0.448 0.464 0.528 0.451];
+   * phi   = [0.507 0.511 0.527 0.549 0.557 0.562 0.560 0.565 0.568 0.565 0.549 0.601 0.580 0.581 0.584 0.587 0.582];
+   *
+   * % MSRS periods
+   * t_nga = [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.0 1.5 2.0 3.0 4.0 5.0 7.5 10.0];
+   *
+   * % interpolate at MSRS periods
+   * sigma_nga = interp1(t,sigma,t_nga,'pchip','extrap');
+   * tau_nga   = interp1(t,tau,t_nga,'pchip','extrap');
+   * phi_nga   = interp1(t,phi,t_nga,'pchip','extrap');
+   *
+   * fprintf('%6s,%7s,%7s,%7s\n','T','sigma','tau','phi');
+   * for i=1:length(t_nga)
+   *     fprintf('%6.3f,%.5f,%.5f,%.5f\n',t_nga(i),sigma_nga(i),tau_nga(i),phi_nga(i));
+   * end
+   * </pre>
+   */
+
+  @Override
+  public LogicTree<GroundMotion> calc(GmmInput in) {
+    /* Eq. 1 */
+    // SA = G1 * G2 * G3 * G4 * G5 * SA_NORM * εγ
+
+    var fault_style = GmmUtils.rakeToFaultStyle_NSHMP(in.rake);
+
+    double sa_norm = spectralShape(period, in.Mw, in.rRup);
+    double G1 = filterG1(in.Mw, fault_style);
+    double G2 = filterG2(in.Mw, in.rRup);
+
+    // PGA is apparently G1 * G2 ???
+
+    double G3 = filterG3(period, in.rRup);
+    double G4 = filterG4(period, in.Mw, in.rRup, fault_style);
+    double G5 = filterG5(period, in.vs30);
+
+    /**
+     * <pre>
+     * Questions:
+     * 1. Is a "G6" filter implied? To add the basin/deep sediment effects to the
+     * Vs30 amplification? No equation is provided but factors (relative to Vs30
+     * amplification) are provided in Table S2 (as a function of Vs30 and T)
+     *
+     * But then use Eq. 16 to convert Z2.5 to effective Vs30 to look up the basin effect???
+     * </pre>
+     */
+    // "G6" for basin effect? look up from table?
+
+    double μ = 1.0;
+    double σ = coeffs.σ;
+
+    return GroundMotions.createTree(μ, σ);
+  }
+
+  /* Spectral shape constants (Figure 3) */
+  private static final double m1 = -0.0012;
+  private static final double m2 = -0.38;
+  private static final double m3 = 4.08;
+
+  private static final double a1 = 0.017;
+  private static final double a2 = 1.27;
+  private static final double a3 = 0.0001;
+
+  private static final double Dsp = 0.75;
+
+  private static final double t1 = 0.001;
+  private static final double t2 = 0.59;
+  private static final double t3 = -2.45;
+
+  private static final double s1 = 0.0;
+  private static final double s2 = 0.077;
+  private static final double s3 = 0.3251;
+
+  private static final double ζ = 2;
+
+  /* 5% damped response spectral shape, SA_norm (Figure 3) */
+  private static final double spectralShape(double T, double Mw, double rRup) {
+    double μ_mr = m1 * rRup + m2 * Mw + m3;
+    double I_mr = (a1 * Mw + a2) * exp(a3 * rRup);
+    double S_mr = s1 * rRup - (s2 * Mw + s3);
+    double Tsp0 = max(0.3, abs(t1 * rRup + t2 * Mw + t3));
+
+    double t_ratio = pow(T / Tsp0, ζ);
+    double t_μ_S = (log(T) + μ_mr) / S_mr;
+    return I_mr * exp(-0.5 * t_μ_S * t_μ_S) +
+        sqrt((1.0 - t_ratio) * (1.0 - t_ratio) + 4.0 * Dsp * Dsp * t_ratio);
+  }
+
+  /* Model coefficients from Graizer & Kalkan (2016) Table 3 */
+  private static final double c1 = 0.14;
+  private static final double c2 = -6.25;
+  private static final double c3 = 0.37;
+  private static final double c4 = 2.237;
+  private static final double c5 = -7.542;
+  private static final double c6 = -0.125;
+  private static final double c7 = 1.19;
+  private static final double c8 = -6.15;
+  private static final double c9 = 0.6;
+
+  /* Constants from Graizer (2018) */
+  private static final double β = 3.5; // defined p. 385
+  private static final double b1 = 0.79; // inferred from matlab code
+  private static final double b2 = 0.27; // inferred from matlab code
+  private static final double b3 = 0.444; // inferred from matlab code
+  private static final double λ = 1.36; // inferred from matlab code
+  private static final double c21 = 0.000854; // inferred from matlab code
+  private static final double c22 = 1.0513; // inferred from matlab code
+
+  /* Filter G1 - magnitude-scaling */
+  private static final double filterG1(double Mw, FaultStyle fault_style) {
+    /*
+     * k_scale is style-of-faulting factor from Graizer and Kalkan (2016) Eq. 3
+     */
+    double k_scale =
+        (fault_style == REVERSE) ? 1.28 : (fault_style == REVERSE_OBLIQUE) ? 1.14 : 1.0;
+
+    return (Mw < 5.5) ? (c21 * exp(c22 * Mw) * k_scale) : (c1 * atan(Mw + c2) + c3) * k_scale;
+  }
+
+  /* Filter G2 - distance scaling */
+  private double filterG2(double Mw, double rRup) {
+    /* Eq. 4 */
+    double R2 = c4 * Mw + c5;
+    double D2 = c6 * cos(c7 * (Mw + c8)) + c9;
+
+    /* Eq. 3 */
+    double w1 = (1 - sqrt(50 / R2)) * (1 - sqrt(50 / R2));
+    double w2 = 4 * D2 * D2 * sqrt(50 / R2);
+    double w3 = (1 - (50 / R2)) * (1 - (50 / R2));
+    double w4 = 4 * D2 * D2 * (50 / R2);
+    double w = sqrt((w1 + w2) / (w3 + w4));
+
+    /* Eq. 3 */
+    double g1 = (1 - sqrt(rRup / R2)) * (1 - sqrt(rRup / R2));
+    double g2 = 4 * D2 * D2 * (rRup / R2);
+    double g3 = 4 * D2 * D2 * sqrt(rRup / R2);
+
+    return (rRup < 50.0) ? 1 / sqrt(g1 + g2) : w / sqrt(g1 + g3);
+  }
+
+  /* Filter G3 - Apparent Attenuation of Spectral Accelerations */
+  private double filterG3(double T, double rRup) {
+    double freq = 1.0 / T;
+
+    // TODO: Is this the right formulation of QSA to use?
+    /* Qsa from Eq. 8 (as a function of freq) */
+    double p = (freq < 1.0) ? 0.90 : (freq < 10.0) ? 0.91 : 0.95;
+    double Qsa = 120 * pow(freq, p);
+
+    return exp(-freq * rRup / β / Qsa);
+  }
+
+  /* Filter G4 - Fault-style scaling */
+  private double filterG4(double T, double Mw, double rRup, FaultStyle fault_style) {
+    // coefficients b1, b2, b3, and λ inferred from matlab code
+    if ((fault_style != FaultStyle.REVERSE) || (Mw < 4.0)) {
+      return 1.0;
+    }
+    double dist = (rRup <= 90.0) ? rRup : 90;
+
+    /* Eq. 12 */
+    double CSoF0 = (b1 + b2 / (1 + (T / 7.0) * (T / 7.0))) *
+        (b3 + 1 / (1 + pow(dist / 100.0, λ)));
+
+    return (Mw >= 5.5) ? CSoF0 : 1 + (CSoF0 - 1) * exp(Mw - 5.5);
+  }
+
+  /* Filter G5 - Site-response term */
+  private double filterG5(double T, double vs30) {
+    /**
+     * <pre>
+     * Questions:
+     * 1. is f in Eq. 15 frequency or something else that doesn't seem to be defined?
+     *
+     * 2. Should Eq. 15 reproduce Table S1: VS30 site amplification factors relative
+     * to the western hard rock of VS30 = 1100 m/s? (it doesn't seem to)
+     *
+     * </pre>
+     */
+    double freq = 1.0 / T;
+
+    /*
+     * amplification coefficient constant for vs30 below 180 m/s and above 1100
+     * m/s
+     */
+    double vs = max(180.0, min(1100.0, vs30));
+    double vss = max(275.0, vs);
+
+    /* Eq. 15 */
+    double kvs = -0.5 * log(vs / 1100.0);
+    double fvs = vss / 120.0 - 2.0;
+
+    return 1 + kvs / sqrt((1 - (fvs / freq)) * (1 - (fvs / freq)) + 1.96 * (fvs / freq));
+  }
+}

src/main/resources/gmm/coeffs/GK17.csv

+     T,  sigma,    tau,    phi
+ 0.010,0.63100,0.37500,0.50700
+ 0.020,0.63300,0.37500,0.51100
+ 0.030,0.64069,0.37845,0.51851
+ 0.050,0.66893,0.40170,0.53414
+ 0.075,0.71300,0.45600,0.54900
+ 0.100,0.73400,0.47700,0.55700
+ 0.150,0.72000,0.45100,0.56200
+ 0.200,0.69300,0.40700,0.56000
+ 0.250,0.68437,0.38956,0.56052
+ 0.300,0.67811,0.37633,0.56181
+ 0.400,0.67300,0.36500,0.56500
+ 0.500,0.69300,0.39800,0.56800
+ 0.750,0.73200,0.46500,0.56500
+ 1.000,0.73800,0.49200,0.54900
+ 1.500,0.76419,0.49776,0.57500
+ 2.000,0.78300,0.50100,0.60100
+ 3.000,0.77400,0.51200,0.58000
+ 4.000,0.73300,0.44800,0.58100
+ 5.000,0.74600,0.46400,0.58400
+ 7.500,0.78675,0.52450,0.58689
+10.000,0.73600,0.45100,0.58200