From e1da043eb34ed7564aaa89de7d376ade31a1225b Mon Sep 17 00:00:00 2001
From: Peter Powers <pmpowers@usgs.gov>
Date: Wed, 15 Jun 2022 22:20:50 -0600
Subject: [PATCH] kuehn seattle

---
 .../gov/usgs/earthquake/nshmp/gmm/Gmm.java    |  7 ++
 .../earthquake/nshmp/gmm/KuehnEtAl_2020.java  | 68 ++++++++++---------
 2 files changed, 43 insertions(+), 32 deletions(-)

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 62b0ffdc..6cb01c5b 100644
--- a/src/main/java/gov/usgs/earthquake/nshmp/gmm/Gmm.java
+++ b/src/main/java/gov/usgs/earthquake/nshmp/gmm/Gmm.java
@@ -533,6 +533,13 @@ public enum Gmm {
       KuehnEtAl_2020.COEFFS,
       KuehnEtAl_2020.CONSTRAINTS_INTERFACE),
 
+  /** @see KuehnEtAl_2020 */
+  KBCG_20_SEATTLE_M9_INTERFACE_BASIN(
+      KuehnEtAl_2020.SeattleInterfaceBasin.class,
+      KuehnEtAl_2020.SeattleInterfaceBasin.NAME,
+      KuehnEtAl_2020.COEFFS,
+      KuehnEtAl_2020.CONSTRAINTS_INTERFACE),
+
   /** @see KuehnEtAl_2020 */
   KBCG_20_CASCADIA_SLAB_BASIN(
       KuehnEtAl_2020.CascadiaSlabBasin.class,
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 03d06ae9..7db6a4f9 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
@@ -87,6 +87,10 @@ public abstract class KuehnEtAl_2020 implements GroundMotionModel {
    *
    * Value of k2 (hard-coded in Python) for PGV is a typo, the value should be
    * -1.955 but -1.995 is used
+   *
+   * TODO cascadia automatically include kuehn seattle behavior, seattle flag
+   * only for M9 : means by default cascadia interface is looking for site term
+   *
    */
 
   static final String NAME = "Kuehn et al. (2020)";
@@ -180,34 +184,14 @@ public abstract class KuehnEtAl_2020 implements GroundMotionModel {
         θ1 = (slab) ? coeffs.get("theta_1_slab_reg_Ca") : coeffs.get("theta_1_if_reg_Ca");
         θ6_2 = coeffs.get("theta_6_2_reg_Ca");
         θ7 = coeffs.get("theta_7_reg_Ca");
-
-        // if (basin == Basin.SEATTLE) {
-        // θ11 = 0;
-        // θ11S = coeffs.get("mean_residual_Seattle_basin");
-        // θ12 = 0;
-        // } else {
-        // θ11 = coeffs.get("intercept_Ca");
-        // θ11S = coeffs.get("mean_residual_Seattle_basin");
-        // θ12 = coeffs.get("slope_Ca");
-        // }
-
       } else if (zone == ALASKA) {
-
         θ1 = (slab) ? coeffs.get("theta_1_slab_reg_Al") : coeffs.get("theta_1_if_reg_Al");
         θ6_2 = coeffs.get("theta_6_2_reg_Al");
         θ7 = coeffs.get("theta_7_reg_Al");
-        // θ11 = 0;
-        // θ11S = 0;
-        // θ12 = 0;
-
-      } else {
-        // global
+      } else { // global
         θ1 = (slab) ? coeffs.get("mu_theta_1_slab") : coeffs.get("mu_theta_1_if");
         θ6_2 = coeffs.get("mu_theta_6");
         θ7 = coeffs.get("mu_theta_7");
-        // θ11 = 0;
-        // θ11S = 0;
-        // θ12 = 0;
       }
 
       θ10 = 0.0; // Table 4.1: fixed to 0
@@ -252,12 +236,12 @@ public abstract class KuehnEtAl_2020 implements GroundMotionModel {
   public LogicTree<GroundMotion> calc(GmmInput in) {
 
     double pgaRef = exp(calcPgaRef(coeffsPGA, in.Mw, in.rRup, in.zTor));
-    double μ = calcMean(coeffs, basin(), seattle(),
+    double μ = calcMean(coeffs, slab(), 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(), seattle(),
+      double μPga = calcMean(coeffsPGA, slab(), basin(), seattle(),
           in.Mw, in.rRup, in.zTor, in.vs30, in.z2p5, pgaRef);
       μ = max(μ, μPga);
     }
@@ -288,6 +272,7 @@ public abstract class KuehnEtAl_2020 implements GroundMotionModel {
 
   private static double calcMean(
       Coefficients c,
+      boolean slab,
       boolean basin,
       boolean seattle,
       double Mw,
@@ -302,7 +287,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, seattle) : 0.0;
+    double fBasin = basin ? basinTerm(c, vs30, z2p5, slab, seattle) : 0.0;
 
     return c.θ1 + fMag + fGeom + fDepth + fAtten + fSite + fBasin;
   }
@@ -375,6 +360,7 @@ public abstract class KuehnEtAl_2020 implements GroundMotionModel {
       Coefficients c,
       double vs30,
       double z2p5,
+      boolean slab,
       boolean seattle) {
 
     if (Double.isNaN(z2p5)) {
@@ -383,15 +369,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;
+    double fb = 0.0; // Seattle cap
+    if (seattle) {
+      fb = c.θ11S;
+      if (!slab && z2p5 > 6.0 & c.imt.period() > 1.9) {
+        fb = log(2.0); // M9 scaling
+      }
+    } else {
+      fb = Math.min(c.θ11 + c.θ12 * δlnZ, c.θ11S); // Seattle cap
     }
-
     return fb;
   }
 
@@ -532,6 +518,24 @@ public abstract class KuehnEtAl_2020 implements GroundMotionModel {
     }
   }
 
+  static final class SeattleInterfaceBasin extends CascadiaInterface {
+    static final String NAME = CascadiaInterface.NAME + " (Seattle basin)";
+
+    SeattleInterfaceBasin(Imt imt) {
+      super(imt);
+    }
+
+    @Override
+    boolean basin() {
+      return true;
+    }
+
+    @Override
+    boolean seattle() {
+      return true;
+    }
+  }
+
   static final class CascadiaSlabBasin extends CascadiaSlab {
     static final String NAME = CascadiaSlab.NAME + " (basin)";
 
-- 
GitLab