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 9b24d1be6afe746d00d028d7db6a5b06e2dce49f..b96fd53e43ce2f2a430095ac9be3527e621f1306 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
@@ -135,13 +135,14 @@ public abstract class AbrahamsonGulerce_2020 implements GroundMotionModel {
     COEFFS_AK_ADJUSTED = new CoefficientContainer("nga-sub-ak-interface-adjustment.csv");
   }
 
+  static final double VS30_ROCK = 1000.0;
+  static final double PHI_AMP_SQ = 0.09;
+
   private static final double C1S = 7.5;
   private static final double C4 = 10.0;
   private static final double C = 1.88;
   private static final double N = 1.18;
-
   private static final double VSS_MAX = 1000.0;
-  private static final double VS30_ROCK = 1000.0;
 
   private static final double D0 = 0.47;
 
@@ -401,33 +402,41 @@ public abstract class AbrahamsonGulerce_2020 implements GroundMotionModel {
     double Ï„lin = D0;
     double Ï„Sq = Ï„lin * Ï„lin;
 
-    double φSq = φlinSq(c.d1, c.d2, rRup);
+    double φSq = calc_φlinSq(c.d1, c.d2, rRup);
 
     double vsS = calcVsStar(vs30);
 
     /* Add nonlinear effects. */
     if (vsS < c.vlin) {
 
-      /* φamp^2 = 0.3 * 0.3 = 0.09 */
-      double φbSq = φSq - 0.09;
+      double φbSq = φSq - PHI_AMP_SQ;
 
-      double φlinSqPga = φlinSq(cPga.d1, cPga.d2, rRup);
-      double φbSqPga = φlinSqPga - 0.09;
+      double φlinSqPga = calc_φlinSq(cPga.d1, cPga.d2, rRup);
+      double φbSqPga = φlinSqPga - PHI_AMP_SQ;
 
-      double dSite = c.b * pgaRock *
-          ((1.0 / (pgaRock + C * pow((vs30 / c.vlin), N))) - (1.0 / (pgaRock + C)));
+      double dSite = calc_dSite(c.b, c.vlin, vs30, pgaRock);
       double dSiteSq = dSite * dSite;
 
-      φSq += dSiteSq * φbSqPga + 2 * dSite * sqrt(φbSqPga) * sqrt(φbSq) * c.ρW;
+      φSq += dSiteSq * φbSqPga + 2.0 * dSite * sqrt(φbSqPga) * sqrt(φbSq) * c.ρW;
 
       /* Ï„linPga = Ï„lin = D0 */
-      τSq += dSiteSq * τSq + 2 * dSite * τSq * c.ρB;
+      τSq += dSiteSq * τSq + 2.0 * dSite * τSq * c.ρB;
     }
 
     return sqrt(τSq + φSq);
   }
 
-  private static double φlinSq(double d1, double d2, double rRup) {
+  /*
+   * Return the partial derivative of the natural log of soil amplification
+   * relative to a hard rock site condition ln(SA[Vs30=1000]). This method is
+   * exposed for use in other GMM aleatory variability models.
+   */
+  static double calc_dSite(double b, double vLin, double vs30, double pgaRock) {
+    return b * pgaRock *
+        ((1.0 / (pgaRock + C * pow((vs30 / vLin), N))) - (1.0 / (pgaRock + C)));
+  }
+
+  static double calc_φlinSq(double d1, double d2, double rRup) {
     /* Eq. 5.2 */
     double φ = d1;
     if (rRup > 450.0) {