Based on Kuehn's paper at https://peer.berkeley.edu/sites/default/files/2020_04_kuehn_final.pdf
With GNZL projects' UHS data, we noticed a NaN issue with branches that include K_20_NZ in the Wellington region where there are enough faults data. (but not limited to Wellington).
Hence, we decided to do some investigation ourselves and see what is happening.
- Check the Hazard Curve in Wellington and see if we have any data at a certain exceedance rate
For instance, draw a line at a certain exceedance rate and see if we can get IM values(This is not from GNZL as we are currently re-generating GNZL data but these are in Wellington)
- Check which branch contains NaN at what exceedance rate level
For example, with Wellington_760_0.065_0p325 => Wellington, Vs30 = 760m/s, Z1.0 = 0.065(km) Z2.5 = 0.325(km)
We can see there are five branches are having NaN issue at exceedance rate = 0.0001 and 0.0002 and followed by the IM periods that is NaN
Exceedance rate = 0.0001, periods at 0.075, 0.1, 0.12, 0.15, 0.17, 0.2, 0.25 and 0.3 are NaN
Exceedance rate = 0.0002, periods at 0.1, 0.12, 0.15, 0.17, 0.2 and 0.25 are NaN
Also, we noticed that the K_20_NZ was the common model that gives us NaN
- Decide to use FLT IMDB to plots with Wellington station, every fault and only pSA's IM values. - Continue to use the same inputs, Wellington, Vs30 = 760, Z1.0 = 0.065, Z2.5 = 0.325
Here are plots, regardless of the faults, we can see the zig-zag shape for K_20NZ and for HikWgtnmin and HikWgtnmax faults, K_20_NZ stays way higher than every other model.
Finally, we decide to plot using Kuehn 2020 model that is in the openquake directly to plot and compare with the one on paper.
The following plot is from Kuehn's paper, Figure 7.19 and the following values are used as inputs during the manual calculation with OQ's Kuehn model(Without using any wrappers):Arg name Value mag 9.0 rrup 75 ztor 10 vs30 400 z1pt0 72.0386568757918 (in m) z2pt5 1.337495657487769 (in km) (Paper says Delta Z1.0 = Delta Z2.5 = 0 but we cannot pass the delta Z values, hence we manually figured out these by checking the reference Z value from Vs30)
Getting Delta Z values - Used the following references
(from Kuehn's paper)
(With Brendon)
def get_basin_response_term(C, region, vs30, z_value):
"""
Returns the basin response term, based on the region and the depth
to a given velocity layer
:param numpy.ndarray z_value:
Basin depth term (Z2.5 for JPN and CAS, Z1.0 for NZL and TWN)
"""
# Basin term only defined for the four regions: Cascadia, Japan,
# New Zealand and Taiwan
assert region in ("CAS", "JPN", "NZL", "TWN")
# Get c11, c12 and Z-model (same for both interface and inslab events)
c11 = C[REGION_TERMS_IF[region]["c11"]]
c12 = C[REGION_TERMS_IF[region]["c12"]]
CZ = Z_MODEL[region]
brt = np.zeros_like(z_value)
mask = z_value > 0.0
if not np.any(mask):
# No basin amplification to be applied
return 0.0
brt[mask] = c11 + c12 * (np.log(z_value[mask]) -
_get_ln_z_ref(CZ, vs30[mask]))
return brt
By looking at the first two reference images, I believe the second photo's Delta Z is equivalent to
then by having Delta Z1.0 = Delta Z2.5 = 0
f_basin(Z) becomes just a theta_11? which is equivalent to c11 from the function above.
So we can use something like this instead,
def get_basin_response_term(C, region, vs30, z_value):
"""
Returns the basin response term, based on the region and the depth
to a given velocity layer
:param numpy.ndarray z_value:
Basin depth term (Z2.5 for JPN and CAS, Z1.0 for NZL and TWN)
"""
# Basin term only defined for the four regions: Cascadia, Japan,
# New Zealand and Taiwan
assert region in ("CAS", "JPN", "NZL", "TWN")
# Get c11, c12 and Z-model (same for both interface and inslab events)
c11 = C[REGION_TERMS_IF[region]["c11"]]
c12 = C[REGION_TERMS_IF[region]["c12"]]
CZ = Z_MODEL[region]
brt = np.zeros_like(z_value)
mask = z_value > 0.0
if not np.any(mask):
# No basin amplification to be applied
return 0.0
# brt[mask] = c11 + c12 * (np.log(z_value[mask]) -
# _get_ln_z_ref(CZ, vs30[mask]))
# As if delta Z is 0, c12 * 0 means nothing
brt[mask] = c11
return brt
Outputs
- Comparison plots by overlapped images
- Seeing plots together
The plot from the Kuehn's paper
The plot using openquake's Kuehn model directly with the same inputs for NZL region
By looking at plots above, between the periods 0.01 ~ 1.0 seems reasonably align but not after period = 1.0
Bugs we found so far
1. Equation 4.8
Equation for Vs30 <= K1, paper says ln[PGA1100 + c * ln(Vs30 / k1) ^n]. However, both moel within openquake and the author's code does not include ln for (Vs30/k1)^n, so we think this is a typo/bug in the paper.
2. Equation 4.8
Two different equation based on whether Vs30 <= K1 or Vs30 > k1. However, author's code says Vs < k1 or Vs >= k1
OQ's Kuehn also compares Vs30 > k1 instead of Vs >= k1 in author's code.
The Paper also says Vs30 > k1. So we think this is another typo/bug in the paper.
Potential fix to match plots with paper
With getting a Magnitude Term using
The bilinear magnitude scaling has breakpoints at Mb,if and Mb,slab that are assigned based on the geometry of the subduction zone.
And for the long periods (T > 1 sec), some adjustments need to be done to Mb,if for interface events at T > 1 sec.
Where the following code block shows the adjustments by adding C["dm_b"] if its a SUBDUCTION_INTERFACE.
m_break = m_b + C["dm_b"] if trt == const.TRT.SUBDUCTION_INTERFACE else m_b
However, by looking at the author's source code.
mbreak <- (1 - fs) * (mb + dmb) + fs * mb (line at 385)
So when tectonic type is SUBDUCTION INTERFACE mbreak is (1 - 0) * (mb + dmb) + 0 * mb => mb + dmb
and this dmb is coming from
dmb <- interp_dmb(period_used, reg) (line at 327)
and the function interp_dmb looks like the following photo and as you can see in the comment, adjustment should only be done for Japan and South America instead of just when the Tectonic Type is SUBDUCTION INTERFACE.
Hence, we need to add extra conditions to the python code
From
m_break = m_b + C["dm_b"] if trt == const.TRT.SUBDUCTION_INTERFACE else m_b
To
m_break = m_b + C["dm_b"] if trt == const.TRT.SUBDUCTION_INTERFACE and self.region in ("JPN", "SAM") else m_b
which only makes adjustments to SUBDUCTION_INTERFAce and JAPAN or SOUTH AMERICA instead of every region.