diff --git a/src/main/java/gov/usgs/earthquake/nshmp/gmm/Gmm.java b/src/main/java/gov/usgs/earthquake/nshmp/gmm/Gmm.java
index d749ab0cd0a8f4e89b732a56f9b59273704efede..51137d860a0f39c8ee1066f900caeacfb8ca6a51 100644
--- a/src/main/java/gov/usgs/earthquake/nshmp/gmm/Gmm.java
+++ b/src/main/java/gov/usgs/earthquake/nshmp/gmm/Gmm.java
@@ -250,12 +250,12 @@ public enum Gmm {
SadighEtAl_1997.COEFFS_BC_HI,
SadighEtAl_1997.CONSTRAINTS),
- // /** @see Grazier_2018 */
- // GRAIZER_2018(
- // Graizer_2018.class,
- // Graizer_2018.NAME,
- // Graizer_2018.COEFFS,
- // Graizer_2018.CONSTRAINTS),
+ /** @see Grazier_2018 */
+ GRAIZER_2018(
+ Graizer_2018.class,
+ Graizer_2018.NAME,
+ Graizer_2018.COEFFS,
+ Graizer_2018.CONSTRAINTS),
/* Subduction Interface and Slab WUS 2008 2014 2018, AK 2007 */
diff --git a/src/main/java/gov/usgs/earthquake/nshmp/gmm/Graizer_2018.java b/src/main/java/gov/usgs/earthquake/nshmp/gmm/Graizer_2018.java
index 481c1b371100ded629b133036c049f85683d5044..5f8e799b8a7e60028ae5cc06671da9a750964777 100644
--- a/src/main/java/gov/usgs/earthquake/nshmp/gmm/Graizer_2018.java
+++ b/src/main/java/gov/usgs/earthquake/nshmp/gmm/Graizer_2018.java
@@ -29,7 +29,7 @@ import gov.usgs.earthquake.nshmp.gmm.GmmInput.Constraints;
import gov.usgs.earthquake.nshmp.tree.LogicTree;
/**
- * TODO:
+ * Outstanding questions:
*
* * is Q_SA, for G3, defined by Eq. 8?
*
@@ -45,7 +45,11 @@ import gov.usgs.earthquake.nshmp.tree.LogicTree;
* * Attempt to compute values in Table S1 and/or S2 using the manual regression
* code blocks in Matlab (so far, the matlab code does not reproduce the tables)
*
- *
+ * * "Residuals" calculated in matlab are not documented in published paper
+ * (they are mentioned)
+ */
+
+/**
* 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
@@ -128,11 +132,11 @@ public final class Graizer_2018 implements GroundMotionModel {
* 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
+ * The following Matlab code was used to compute sigma values at the MPRS
* periods:
*
* <pre>
- * % From Table 3 in Graizer (2018)
+ * % From Table 1 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];
@@ -155,8 +159,6 @@ public final class Graizer_2018 implements GroundMotionModel {
@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);
@@ -164,25 +166,15 @@ public final class Graizer_2018 implements GroundMotionModel {
double G1 = filterG1(in.Mw, fault_style); // M-scaling
double G2 = filterG2(in.Mw, in.rRup); // R-scaling
- // PGA is apparently G1 * G2 ???
+ // Matlab code GK17PGA.m returns G1 * G2
double G3 = filterG3(period, in.rRup); // Apparent attenuation scaling
double G4 = filterG4(period, in.Mw, in.rRup, fault_style); // Style of
// faulting style
- double G5 = filterG5(period, in.vs30); // Site-scaling
-
- /**
- * <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 G5 = filterG5(period, in.vs30, in.z2p5); // Site-scaling
+ /* Eq. 1 */
+ // SA = G1 * G2 * G3 * G4 * G5 * SA_NORM * εγ
double μ = log(G1 * G2 * G3 * G4 * G5 * sa_norm);
double σ = coeffs.σ;
@@ -212,6 +204,10 @@ public final class Graizer_2018 implements GroundMotionModel {
/* 5% damped response spectral shape, SA_norm (Figure 3) */
private static final double spectralShape(double T, double Mw, double rRup) {
+ /*
+ * Matlab code GK17SA.m uses truncated M and R values when calculating the
+ * spectral shape, this is not mentioned in the paper
+ */
double μ_mr = m1 * rRup + m2 * Mw + m3;
double I_mr = (a1 * Mw + a2) * exp(a3 * rRup);
double S_mr = s1 * rRup - (s2 * Mw + s3);
@@ -223,7 +219,7 @@ public final class Graizer_2018 implements GroundMotionModel {
sqrt((1.0 - t_ratio) * (1.0 - t_ratio) + 4.0 * Dsp * Dsp * t_ratio);
}
- /* Model coefficients from Graizer & Kalkan (2016) Table 3 */
+ /* Model coefficients from Graizer & Kalkan (2016) Table 3 and GK17PGA.m */
private static final double c1 = 0.14;
private static final double c2 = -6.25;
private static final double c3 = 0.37;
@@ -240,17 +236,29 @@ public final class Graizer_2018 implements GroundMotionModel {
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
+ private static final double c21 = 0.000854; // inferred from matlab GK17PGA.m
+ private static final double c22 = 1.0513; // inferred from matlab GK17PGA.m
/* 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
+ * Magnitude scaling filter is from previous model (Graizer and Kalkan
+ * (2016) Eq. 3) for larger magnitudes with new small-magnitude scaling.
+ *
+ * Paper implies k_scale is style-of-faulting factor from Graizer and Kalkan
+ * (2016) and does not define new coefficients c21 and c22, which must be
+ * inferred from matlab code GK17PGA.m (F1 term in matlab).
+ *
+ * "k_scale" as used in matlab code GK17PGA.m differs from Graizer and
+ * Kalkan (2016), which was dependent on fault style.
*/
double k_scale =
(fault_style == REVERSE) ? 1.28 : (fault_style == REVERSE_OBLIQUE) ? 1.14 : 1.0;
+ double F_matlab = 1.0; // F = 1.0 for geomean and RotD50
+ double k_scale_matlab = F_matlab / 1.12;
+
+ /* matlab code GK17PGA.m uses M <= 5.5 and M > 5.5 */
return (Mw < 5.5) ? (c21 * exp(c22 * Mw) * k_scale) : (c1 * atan(Mw + c2) + c3) * k_scale;
}
@@ -262,11 +270,19 @@ public final class Graizer_2018 implements GroundMotionModel {
/* 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 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));
+ // /*
+ // * w from Matlab GK17PGA.m uses 1.96 instead of Mw-dependent "4 * D2 * D2"
+ // * and pulls -sqrt out into log-space
+ // */
+ // double ml_w2 = 1.96 * sqrt(50. / R2); // matlab
+ // double w_matlab = exp(-0.5 * log((w3 + w4) / (w1 + ml_w2)));
+ // System.out.printf("Paper w = %.6e\nMatlab w = %.6e\n", w, w_matlab);
+
/* Eq. 3 */
double g1 = (1 - sqrt(rRup / R2)) * (1 - sqrt(rRup / R2));
double g2 = 4 * D2 * D2 * (rRup / R2);
@@ -279,23 +295,37 @@ public final class Graizer_2018 implements GroundMotionModel {
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);
+ double rRef = 120.0;
+
+ double Qsa = rRef * pow(freq, p);
+
+ /*
+ * Q-term (Qsa factor) in matlab code GK17SA.m has only a slight resemblance
+ * to Eq. 8 (and the G3 filter too, for that matter) and includes an added
+ * dependence on magnitude (Eq. 10?)
+ *
+ * β, with a value of 3.5 km defined in the paper does not appear in the
+ * matlab code
+ */
+ // double p_matlab = (freq <= 1.0) ? 0.90 : (freq <= 10.0) ? 0.91 : 0.95;
+ // double rRef_matlab = (rRup <= 140.0) ? 120.0 : 120. + (140. - rRup) /
+ // 100.;
+ // double Qsa_matlab = (0.2759 + 0.1379 * Mw) * rRef * 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
+ /* Eqs. 13 and 14 */
if ((fault_style != FaultStyle.REVERSE) || (Mw < 4.0)) {
return 1.0;
}
double dist = (rRup <= 90.0) ? rRup : 90;
- /* Eq. 12 */
+ /* Eq. 12 - coefficients b1, b2, b3, and λ inferred from matlab code */
double CSoF0 = (b1 + b2 / (1 + (T / 7.0) * (T / 7.0))) *
(b3 + 1 / (1 + pow(dist / 100.0, λ)));
@@ -303,7 +333,7 @@ public final class Graizer_2018 implements GroundMotionModel {
}
/* Filter G5 - Site-response term */
- private double filterG5(double T, double vs30) {
+ private double filterG5(double T, double vs30, double z2p5) {
/**
* <pre>
* Questions:
@@ -312,11 +342,15 @@ public final class Graizer_2018 implements GroundMotionModel {
* 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)
*
+ * 3. Should basin effects (sediment thickness) be added here to this term? Sediment
+ * thickness effects appear be calculated based on residuals in matlab code GK17SA.m
* </pre>
*/
double freq = 1.0 / T;
/*
+ * Eq. 15
+ *
* amplification coefficient constant for vs30 below 180 m/s and above 1100
* m/s
*/
@@ -327,22 +361,91 @@ public final class Graizer_2018 implements GroundMotionModel {
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));
- }
+ double g5_vs30 = 1 + kvs / sqrt((1 - (fvs / freq)) * (1 - (fvs / freq)) + 1.96 * (fvs / freq));
- /* Sediment thickness (Z2.5) */
- private double sedimentThickness(double vs30, double z2p5) {
- // Eq. 16 is a Z2.5 - Vs30 relation
- // ln(z2p5) = 7.108 - 1.145 * ln(vs30)
+ /*
+ * Matlab code in GK17SA.m, "Reducing to hard rock", has a vague similarity
+ * to Eq. 15. Terms Amp_Va and Amp_Vs are calculated and, in the formulation
+ * for Y, Amp_vs is divided by Amp_Va and the vs30 correction (residuals)
+ */
+ double basinAmp = basinAmplification(vs30, z2p5);
- // Matlab script has this relation:
- // ln(z2p5) = 7.521 - 1.233 * ln(VS30)
+ return g5_vs30 * basinAmp;
+ }
- // No equations in the paper are a function of Z2.5.
- // Paper refers to Table S2 for basin amplification (as a function of Vs30
- // and Vs30), which are specified as relative to the site amplification due
- // Vs30
+ /* Basin amplification */
+ private double basinAmplification(double vs30, double z2p5) {
+ /*
+ * The paper does not define an equation for the sediment thickness term
+ *
+ * Eq. 16 is a Z2.5 - Vs30 relation: ln(z2p5) = 7.108 - 1.145 * ln(vs30)
+ *
+ * Matlab code GK17SA.m uses Z2.5 directly and has a different Z2.5 - Vs30
+ * relation in a comment: ln(Z25) = 7.521 - 1.233 * ln(VS30)
+ *
+ * Paper refers to Table S2 for basin amplification (as a function of Vs30
+ * and period, T), which is specified as relative to the site amplification
+ * due to Vs30.
+ */
+ if (Double.isNaN(z2p5)) {
+ return 1.0;
+ }
+ /*
+ * Look up value from table of ratios (vs30_amp / basin_amp?) for vs30 and
+ * z2p5? Or calculate using residuals calculations (BasAmp term) in matlab
+ * code GK17SA.m?
+ *
+ * From matlab, 5 period-dependent coefs: ba0, ba1, ba2, ba3, ba4:
+ *
+ * - BasAmp = ba0 for 0 <= z2p5 <= 0.18
+ *
+ * - BasAmp = ba1 * z2p5^ba2 for 0.18 < z2p5 <= 0.67
+ *
+ * - BasAmp = 1 for 0.67 < z2p5 < 3.15 *** Does this interval make sense?
+ *
+ * - BasAmp = ba3 * z2p5^ba4 for z2p5 >= 3.15
+ *
+ * where BasAmp = ba1 * z2p5^ba2 or BasAmp = ba3 * z2p5^ba4, depending on
+ * value of z2p5
+ */
return 1.0;
}
+
+ // private static final double getPGA_matlab(double Mw, double rRup) {
+ // double F = 1.0; // F = 1.0 for geomean and RotD50
+ // double F1 = 0;
+ // double F2 = 0;
+ //
+ // // This is Equation 2 [G1(M,kscale)], but with kscale replaced with F/1.12
+ // // (instead of style-of-faulting factor from GK16)
+ // // c21 = 0.000854
+ // // c22 = 1.0513
+ // if (Mw > 5.5) { // Eq. 2 (M >= 5.5)
+ // F1 = log((c1 * atan(Mw + c2) + c3) * F / 1.12);
+ // } else {
+ // F1 = log((0.000854 * exp(1.0513 * Mw) * F / 1.12));
+ // }
+ //
+ // // This is Eq. 4 in the paper, Ro is R2, and Do is D2
+ // double Ro = c4 * Mw + c5;
+ // double Do = c6 * cos(c7 * (Mw + c8)) + c9;
+ //
+ // // Bottom of Eq. 3 in the paper ??? in log space without sqrt
+ // // Conversion from body to surface waves W
+ // double W = log((pow((1.0 - 50.0 / Ro), 2) + 4 * (Do * Do) * 50.0 / Ro) /
+ // (pow(1 - sqrt(50 / Ro), 2) + 1.96 * sqrt(50.0 / Ro)));
+ //
+ // if (rRup <= 50.0) {
+ // // attenuation of the order R-1
+ // F2 = -0.5 * log(pow(1 - (rRup / Ro), 2) + 4.0 * Do * Do * (rRup / Ro));
+ // } else {
+ // // attenuation of the order R-0.5
+ // F2 = -0.5 * log(pow(1 - sqrt(rRup / Ro), 2) + 1.96 * sqrt(rRup / Ro)) - 0.5
+ // * W;
+ // }
+ //
+ // return F1 + F2;
+ // }
+
}