From d01e0d91e0ee0676431b6393bdf850e055df5b8a Mon Sep 17 00:00:00 2001
From: Peter Powers <pmpowers@usgs.gov>
Date: Wed, 15 Jun 2022 14:09:30 -0600
Subject: [PATCH] added m9 adjustments

---
 .../nshmp/gmm/AbrahamsonGulerce_2020.java     | 61 +++++++++++++------
 .../earthquake/nshmp/gmm/KuehnEtAl_2020.java  | 39 ++++++++++--
 .../earthquake/nshmp/gmm/ParkerEtAl_2020.java | 59 ++++++++++++++----
 3 files changed, 126 insertions(+), 33 deletions(-)

diff --git a/src/main/java/gov/usgs/earthquake/nshmp/gmm/AbrahamsonGulerce_2020.java b/src/main/java/gov/usgs/earthquake/nshmp/gmm/AbrahamsonGulerce_2020.java
index 23bd5c54..b2e632f0 100644
--- a/src/main/java/gov/usgs/earthquake/nshmp/gmm/AbrahamsonGulerce_2020.java
+++ b/src/main/java/gov/usgs/earthquake/nshmp/gmm/AbrahamsonGulerce_2020.java
@@ -91,6 +91,7 @@ public abstract class AbrahamsonGulerce_2020 implements GroundMotionModel {
 
   private static final double VSS_MAX = 1000.0; // Errata 2021/7/12
   private static final double VS30_ROCK = 1000.0;
+  private static final double Z2P5_REF_AG = 1998.0;
 
   private static final double D0 = 0.47;
 
@@ -177,10 +178,10 @@ public abstract class AbrahamsonGulerce_2020 implements GroundMotionModel {
 
     double pgaRock = 0.0;
     if (in.vs30 < coeffs.vlin) {
-      pgaRock = exp(calcMean(coeffsPGA, slab(), basin(), adjust(), 0.0,
+      pgaRock = exp(calcMean(coeffsPGA, slab(), basin(), seattle(), adjust(), 0.0,
           in.Mw, in.rRup, in.zTor, VS30_ROCK, in.z2p5));
     }
-    double μ = calcMean(coeffs, slab(), basin(), adjust(), pgaRock,
+    double μ = calcMean(coeffs, slab(), basin(), seattle(), adjust(), pgaRock,
         in.Mw, in.rRup, in.zTor, in.vs30, in.z2p5);
 
     double σ = calcSigma(coeffs, coeffsPGA, pgaRock, in.rRup, in.vs30);
@@ -192,12 +193,15 @@ public abstract class AbrahamsonGulerce_2020 implements GroundMotionModel {
 
   abstract boolean basin();
 
+  abstract boolean seattle();
+
   abstract boolean adjust();
 
   private static double calcMean(
       Coefficients c,
       boolean slab,
       boolean basin,
+      boolean seattle,
       boolean adjust,
       double pgaRock,
       double Mw,
@@ -229,7 +233,7 @@ public abstract class AbrahamsonGulerce_2020 implements GroundMotionModel {
     }
 
     /* Eq. 3.7, 3.8 - Site-response scaling */
-    double vsS = getVsStar(vs30);
+    double vsS = calcVsStar(vs30);
     double fSite = c.a12 * log(vsS / c.vlin);
     if (vsS < c.vlin) {
       fSite += -c.b * log(pgaRock + C) + c.b * log(pgaRock + C * pow((vsS / c.vlin), N));
@@ -241,20 +245,23 @@ public abstract class AbrahamsonGulerce_2020 implements GroundMotionModel {
     double fBasin = 0.0;
     if (basin && !Double.isNaN(z2p5)) {
       double z2p5m = z2p5 * 1000.0;
-      /*
-       * Note: PEER report (and xlsx) uses factor rounded to 0.88, fortran code
-       * uses the calculated value. Using 0.878 would improve the match of
-       * lnZ2p5ref to the fortran code (within 0.01%).
-       */
-      double lnZ2p5ref = (vs30 >= 570)
-          ? 7.6
-          : (vs30 > 200)
-              ? 8.52 - 0.88 * log(vs30 / 200.0)
-              : 8.52;
+      double lnZ2p5ref = 8.52;
+      if (vs30 >= 570) {
+        lnZ2p5ref = 7.6;
+      } else if (vs30 > 200) {
+        lnZ2p5ref = 8.52 - 0.88 * log(vs30 / 200.0);
+      }
 
       /* Eq. 3.9, 3.11 - Cascadia basin scaling (Errata 2021/7/12) */
-      double zPrime = (z2p5m + 50.0) / (exp(lnZ2p5ref) + 50.0);
-      fBasin = (log(zPrime) > 0) ? c.a39 * log(zPrime) : 0.0;
+      double lnzPrime = calcLnzPrime(z2p5m, exp(lnZ2p5ref));
+      if (lnzPrime > 0) {
+        fBasin = c.a39 * lnzPrime;
+        if (seattle && z2p5 > 6.0 & c.imt.period() > 1.9) {
+          double lnzPrimeAdj = calcLnzPrime(z2p5m, Z2P5_REF_AG);
+          double m9adj = log(2.0) - c.a39 * lnzPrimeAdj;
+          fBasin += m9adj;
+        }
+      }
     }
 
     /* Eq. 3.1 */
@@ -271,7 +278,12 @@ public abstract class AbrahamsonGulerce_2020 implements GroundMotionModel {
     return μ;
   }
 
-  private static double getVsStar(double vs30) {
+  /* Args in meters. */
+  private static double calcLnzPrime(double z2p5, double z2p5ref) {
+    return log((z2p5 + 50.0) / (z2p5ref + 50.0));
+  }
+
+  private static double calcVsStar(double vs30) {
     return min(vs30, VSS_MAX);
   }
 
@@ -287,7 +299,7 @@ public abstract class AbrahamsonGulerce_2020 implements GroundMotionModel {
 
     double φSq = φlinSq(c.d1, c.d2, rRup);
 
-    double vsS = getVsStar(vs30);
+    double vsS = calcVsStar(vs30);
 
     /* Add nonlinear effects. */
     if (vsS < c.vlin) {
@@ -341,6 +353,11 @@ public abstract class AbrahamsonGulerce_2020 implements GroundMotionModel {
       return false;
     }
 
+    @Override
+    boolean seattle() {
+      return false;
+    }
+
     @Override
     boolean adjust() {
       return false;
@@ -365,6 +382,11 @@ public abstract class AbrahamsonGulerce_2020 implements GroundMotionModel {
       return false;
     }
 
+    @Override
+    boolean seattle() {
+      return false;
+    }
+
     @Override
     boolean adjust() {
       return true;
@@ -389,6 +411,11 @@ public abstract class AbrahamsonGulerce_2020 implements GroundMotionModel {
       return false;
     }
 
+    @Override
+    boolean seattle() {
+      return false;
+    }
+
     @Override
     boolean adjust() {
       return true;
diff --git a/src/main/java/gov/usgs/earthquake/nshmp/gmm/KuehnEtAl_2020.java b/src/main/java/gov/usgs/earthquake/nshmp/gmm/KuehnEtAl_2020.java
index 352ec391..90fd24b2 100644
--- a/src/main/java/gov/usgs/earthquake/nshmp/gmm/KuehnEtAl_2020.java
+++ b/src/main/java/gov/usgs/earthquake/nshmp/gmm/KuehnEtAl_2020.java
@@ -117,6 +117,7 @@ public abstract class KuehnEtAl_2020 implements GroundMotionModel {
   private static final double C = 1.88; // p. 22 !!! Value from python code
   private static final double N = 1.18; // p. 22 !!! Value from python code
   private static final double VSROCK = 1100;
+  private static final double LN_Z2P5_REF_KUEHN = log(200.0);
 
   private static final class Coefficients {
     final Imt imt;
@@ -244,15 +245,19 @@ public abstract class KuehnEtAl_2020 implements GroundMotionModel {
 
   abstract boolean basin();
 
+  abstract boolean seattle();
+
   @Override
   public LogicTree<GroundMotion> calc(GmmInput in) {
 
     double pgaRef = exp(calcPgaRef(coeffsPGA, in.Mw, in.rRup, in.zTor));
-    double μ = calcMean(coeffs, basin(), in.Mw, in.rRup, in.zTor, in.vs30, in.z2p5, pgaRef);
+    double μ = calcMean(coeffs, basin(), seattle(),
+        in.Mw, in.rRup, in.zTor, in.vs30, in.z2p5, pgaRef);
 
     // short periods can't be lower than PGA
     if (coeffs.imt.isSA() && coeffs.imt.period() <= 0.1) {
-      double μPga = calcMean(coeffsPGA, basin(), in.Mw, in.rRup, in.zTor, in.vs30, in.z2p5, pgaRef);
+      double μPga = calcMean(coeffsPGA, basin(), seattle(),
+          in.Mw, in.rRup, in.zTor, in.vs30, in.z2p5, pgaRef);
       μ = max(μ, μPga);
     }
 
@@ -283,6 +288,7 @@ public abstract class KuehnEtAl_2020 implements GroundMotionModel {
   private static double calcMean(
       Coefficients c,
       boolean basin,
+      boolean seattle,
       double Mw,
       double rRup,
       double zTor,
@@ -295,7 +301,7 @@ public abstract class KuehnEtAl_2020 implements GroundMotionModel {
     double fDepth = depthTerm(c, zTor);
     double fAtten = arcCrossingTerm(c, rRup);
     double fSite = siteAmpTerm(c, vs30, pgaRef);
-    double fBasin = basin ? basinTerm(c, vs30, z2p5) : 0.0;
+    double fBasin = basin ? basinTerm(c, vs30, z2p5, seattle) : 0.0;
 
     return c.θ1 + fMag + fGeom + fDepth + fAtten + fSite + fBasin;
   }
@@ -367,7 +373,8 @@ public abstract class KuehnEtAl_2020 implements GroundMotionModel {
   private static double basinTerm(
       Coefficients c,
       double vs30,
-      double z2p5) {
+      double z2p5,
+      boolean seattle) {
 
     // TODO add Seattle flag
 
@@ -378,8 +385,15 @@ public abstract class KuehnEtAl_2020 implements GroundMotionModel {
     double z2p5m = z2p5 * 1000.0;
     double δlnZ = log(z2p5m) - basinZref(vs30);
     double fb = c.θ11 + c.θ12 * δlnZ;
+    fb = Math.min(fb, c.θ11S); // Seattle cap
+
+    if (seattle && z2p5 > 6.0 & c.imt.period() > 1.9) {
+      double δz2p5adj = log(z2p5m) - LN_Z2P5_REF_KUEHN;
+      double m9adj = log(2.0) - (c.θ11 + c.θ12 * δz2p5adj);
+      fb += m9adj;
+    }
 
-    return Math.min(fb, c.θ11S); // Seattle cap
+    return fb;
 
     // switch (region) {
     // case Region.CASCADIA:
@@ -431,6 +445,11 @@ public abstract class KuehnEtAl_2020 implements GroundMotionModel {
     boolean basin() {
       return false;
     }
+
+    @Override
+    boolean seattle() {
+      return false;
+    }
   }
 
   private static abstract class Cascadia extends KuehnEtAl_2020 {
@@ -450,6 +469,11 @@ public abstract class KuehnEtAl_2020 implements GroundMotionModel {
     boolean basin() {
       return false;
     }
+
+    @Override
+    boolean seattle() {
+      return false;
+    }
   }
 
   private static abstract class Alaska extends KuehnEtAl_2020 {
@@ -469,6 +493,11 @@ public abstract class KuehnEtAl_2020 implements GroundMotionModel {
     boolean basin() {
       return false;
     }
+
+    @Override
+    boolean seattle() {
+      return false;
+    }
   }
 
   static class GlobalInterface extends Global {
diff --git a/src/main/java/gov/usgs/earthquake/nshmp/gmm/ParkerEtAl_2020.java b/src/main/java/gov/usgs/earthquake/nshmp/gmm/ParkerEtAl_2020.java
index 8b289b6a..a3e5b737 100644
--- a/src/main/java/gov/usgs/earthquake/nshmp/gmm/ParkerEtAl_2020.java
+++ b/src/main/java/gov/usgs/earthquake/nshmp/gmm/ParkerEtAl_2020.java
@@ -38,6 +38,11 @@ import gov.usgs.earthquake.nshmp.tree.LogicTree;
  * authors recommend {@code zHyp}.
  *
  * <li>The basin depth scaling model is that developed for Cascadia.
+ * Implementations include USGS developed depth tapering and restriction to long
+ * periods ({@code T &geq; 1 s}), as well as a Seattle-specific basin
+ * implementation that makes further adjustments based on M9 simulations for
+ * sites over the deepest depths ({@code Z2.5 > 6 km}) and for long periods
+ * ({@code T &qeq; 2 s}).
  *
  * <li>Implementations are provided for the global model as well as Cascadia and
  * Alaska regionalized models.</ul>
@@ -113,6 +118,7 @@ public abstract class ParkerEtAl_2020 implements GroundMotionModel {
   private static final double B4 = 0.1;
   private static final double V1 = 270.0; // m/s
   private static final double VREF = 760.0; // m/s
+  private static final double LN_Z2P5_REF_PARKER = log(1279.0);
 
   private static final double DB1 = 20.0; // km
   private static final double DB2 = 67.0; // km
@@ -244,12 +250,12 @@ public abstract class ParkerEtAl_2020 implements GroundMotionModel {
 
     double pgaRef = exp(calcMean(coeffsPGA, slab(), in.Mw, in.rRup, in.zTor));
     double μ = calcMean(
-        coeffs, slab(), basin(),
+        coeffs, slab(), basin(), seattle(),
         in.Mw, in.rRup, in.zTor, pgaRef, in.vs30, in.z2p5);
 
     double Ï•Total = phiTotal(coeffs, in.rRup, in.vs30);
     double σ = Maths.hypot(coeffs.τ, ϕTotal);
-    double σε = getσε(coeffs);
+    double σε = calcσε(coeffs);
 
     double[] μs = new double[] { μ + σε, μ, μ - σε };
 
@@ -260,6 +266,8 @@ public abstract class ParkerEtAl_2020 implements GroundMotionModel {
 
   abstract boolean basin();
 
+  abstract boolean seattle();
+
   /* Calc median at 760 reference condition */
   private static double calcMean(
       Coefficients c,
@@ -280,6 +288,7 @@ public abstract class ParkerEtAl_2020 implements GroundMotionModel {
       Coefficients c,
       boolean slab,
       boolean basin,
+      boolean seattle,
       double Mw,
       double rRup,
       double zHyp,
@@ -295,7 +304,7 @@ public abstract class ParkerEtAl_2020 implements GroundMotionModel {
     /* Equation 7 */
     double Fslin = siteTermLinear(c, vs30);
     double Fsnl = siteTermNonLinear(c, pgaRef, vs30);
-    double Fsb = basin ? getSiteBasinTerm(c, vs30, z2p5) : 0.0;
+    double Fsb = basin ? siteBasinTerm(c, vs30, z2p5, seattle) : 0.0;
     double Fs = Fslin + Fsnl + Fsb;
 
     return c.c0 + Fp + Fm + Fd + Fs;
@@ -365,7 +374,11 @@ public abstract class ParkerEtAl_2020 implements GroundMotionModel {
   }
 
   /* Equation 11, 12, 13 */
-  private static double getSiteBasinTerm(Coefficients c, double vs30, double z2p5) {
+  private static double siteBasinTerm(
+      Coefficients c,
+      double vs30,
+      double z2p5,
+      boolean seattle) {
 
     if (Double.isNaN(z2p5)) {
       return 0.0;
@@ -383,15 +396,24 @@ public abstract class ParkerEtAl_2020 implements GroundMotionModel {
         (θ1 * (1 + Maths.erf((log10(vs30) - log10(νμ)) / νσ / sqrt(2))) + θ0);
     double δz2p5 = log(z2p5m) - lnμz2p5;
 
-    double Fb;
+    double Fb = calcFb(c, δz2p5);
+
+    if (seattle && z2p5 > 6.0 & c.imt.period() > 1.9) {
+      double δz2p5adj = log(z2p5m) - LN_Z2P5_REF_PARKER;
+      double m9adj = log(2.0) - calcFb(c, δz2p5adj);
+      Fb += m9adj;
+    }
+
+    return Fb;
+  }
+
+  private static double calcFb(Coefficients c, double δz2p5) {
     if (δz2p5 <= (c.c_e1 / c.c_e3)) {
-      Fb = c.c_e1;
+      return c.c_e1;
     } else if (δz2p5 >= (c.c_e2 / c.c_e3)) {
-      Fb = c.c_e2;
-    } else {
-      Fb = c.c_e3 * δz2p5;
+      return c.c_e2;
     }
-    return Fb;
+    return c.c_e3 * δz2p5;
   }
 
   /* Eqsuation 18, 19, 20 */
@@ -422,7 +444,7 @@ public abstract class ParkerEtAl_2020 implements GroundMotionModel {
   }
 
   /* Equation 27 */
-  private static double getσε(Coefficients c) {
+  private static double calcσε(Coefficients c) {
     double σε;
     if (c.imt == Imt.PGA || c.imt == Imt.PGV || c.imt.period() < c.t1) {
       σε = c.σε1;
@@ -451,6 +473,11 @@ public abstract class ParkerEtAl_2020 implements GroundMotionModel {
     boolean basin() {
       return false;
     }
+
+    @Override
+    boolean seattle() {
+      return false;
+    }
   }
 
   private static abstract class Cascadia extends ParkerEtAl_2020 {
@@ -470,6 +497,11 @@ public abstract class ParkerEtAl_2020 implements GroundMotionModel {
     boolean basin() {
       return false;
     }
+
+    @Override
+    boolean seattle() {
+      return false;
+    }
   }
 
   private static abstract class Alaska extends ParkerEtAl_2020 {
@@ -489,6 +521,11 @@ public abstract class ParkerEtAl_2020 implements GroundMotionModel {
     boolean basin() {
       return false;
     }
+
+    @Override
+    boolean seattle() {
+      return false;
+    }
   }
 
   static final class GlobalInterface extends Global {
-- 
GitLab