Commit 38c3c9fe authored by Graeme Weatherill's avatar Graeme Weatherill
Browse files

Implements a more complete set of magnitude scaling relations

parent a6556a9f
Pipeline #20059 failed with stage
in 1 minute and 8 seconds
......@@ -3,59 +3,13 @@ Contains the core classes to represent an earthquake
"""
import datetime
from typing import Tuple
from math import radians, sin, tan, sqrt
from math import radians, sin, tan
from openquake.hazardlib.geo import Point, PlanarSurface
from openquake.hazardlib.geo.surface.base import BaseSurface
from openquake.hazardlib.source.rupture import ParametricProbabilisticRupture
from shakyground2 import valid
class DummyScalingRelation(object):
"""
A simple magnitude-area scaling relation to use as a placeholder until
updated by the synthetic rupture generator
"""
@staticmethod
def get_rupture_dimensions(
magnitude: float, rake: float = 0.0, dip: float = 90.0, aspect: float = 1.0,
thickness: float = 20.0
) -> Tuple[float, float, float]:
"""
Returns the area (km^2), length (km) and width (km) of the rupture using the PEER
Magnitude-Area scaling relation and constrained by the crustal thickness:
A = 10.0 ^ (magnitude - 4.0)
Args:
magnitude:
Earthquake magnitude
rake:
Rake of the earthquake rupture in degrees (in the range -180 to 180)
dip:
Dip of the earthquake rupture in degrees (in the range 0.0 to 90.0)
aspect:
Aspect ratio of the rupture
thickness:
Seismogenic thickness of the crust (in km)
Returns:
Dimensions of the rupture as a tuple of rupture area (in km^2), rupture length
(in km) and downdip rupture width (in km)
"""
# PEER defined magnitude-area scaling relation returns area in km^2
area = 10.0 ** (magnitude - 4.0)
width = sqrt(area / aspect)
# dz is the vertical dimension of the rupture (in km)
dz = width * sin(radians(dip))
if dz > thickness:
# Rupture is wider than the crustal thickness, so rescale it to the
# maximum width and break the aspect ratio
width = thickness / sin(radians(dip))
length = area / width
return area, length, width
class Earthquake(object):
"""
Shakemap event class requires a minimum of an ID, longitude, latitude,
......@@ -86,14 +40,14 @@ class Earthquake(object):
aspect:
Aspect ratio (length / width) of the rupture
mag_scale_rel:
Magnitude scaling relation
Magnitude scaling relation as instance of the BaseScalingRelation class or None
lsd:
Lower seismogenic depth (in km)
usd:
Upper seismogenic depth (in km)
datetime.date date:
date:
Date of the earthquake (year, month, day) as datetime.date object
datetime.time time:
time:
Time of the earthquake (hour, minute, second) as datetime.time object
surface:
Rupture surface as instance of :class:`openquake.hazardlib.geo.BaseSurface`
......@@ -121,7 +75,48 @@ class Earthquake(object):
date=None,
time=None,
):
"""
Initialize a new Earthquake object
Args:
id:
Unique ID for the earthquake
lon:
Earthquake longitude in decimal degrees (-180 to 180)
lat:
Earthquake latitude in decimal degrees (-90, to 90)
hypo_depth:
Hypocentral depth (in km)
magnitude:
Earthquake magnitude
strike:
Strike of the rupture (in degrees from north)
dip:
Dip of the rupture (in degrees from horizontal)
rake:
Rake of the earthquake rupture in decimal degrees according to the
Aki & Richards (1980) convention (-180 to 180)
aspect:
Aspect ratio (length / width) of the rupture
mag_scaling_relation:
Magnitude scaling relation as instance of the BaseScalingRelation class or None.
Defaults to :class:`shakyground2.magnitude_scaling_relations.PEERScalingRelation`
when missing
upper_seismogenic_depth:
Upper seismogenic depth (in km)
lower_seismogenic_depth:
Lower seismogenic depth (in km)
surface:
Rupture surface as instance of :class:`openquake.hazardlib.geo.BaseSurface`, or None
hypocenter_position:
Hypocentre location within the rupture plane as a tuple of
(fraction of along strike length, fraction of down-dip width)
date:
Date of the earthquake (year, month, day) as datetime.date object
time:
Time of the earthquake (hour, minute, second) as datetime.time object
"""
self.id = id
self.lon = valid.longitude(lon)
self.lat = valid.latitude(lat)
......@@ -134,7 +129,7 @@ class Earthquake(object):
self.usd, self.lsd = valid.seismogenic_thickness(
upper_seismogenic_depth, lower_seismogenic_depth
)
self.mag_scale_rel = mag_scaling_relation or DummyScalingRelation()
self.mag_scale_rel = valid.scaling_relation(mag_scaling_relation)
# Date and time should be parsed as datetime.date and datetime.time
# objects if defined, otherwise none
assert isinstance(date, datetime.date) or date is None
......
"""
Implements magnitude scaling relations for calculating the finite rupture dimensions
"""
import abc
from typing import Tuple
from math import radians, sin, exp, sqrt, pi, log
from scipy.stats import norm, chi2
class BaseScalingRelation(metaclass=abc.ABCMeta):
"""
Abstract base class representing an implementation of a magnitude scaling relation
"""
@abc.abstractmethod
def get_rupture_dimensions(self, magnitude: float, rake: float = 0.0, dip: float = 90.0,
aspect: float = 1.0, thickness: float = 20.0,
epsilon: float = 0.0) -> Tuple[float, float, float]:
"""
Args:
magnitude:
Earthquake magnitude
rake:
Rake of the earthquake rupture in degrees (in the range -180 to 180)
dip:
Dip of the earthquake rupture in degrees (in the range 0.0 to 90.0)
aspect:
Aspect ratio of the rupture
thickness:
Seismogenic thickness of the crust (in km)
Returns:
Dimensions of the rupture as a tuple of rupture area (in km^2), rupture length
(in km) and downdip rupture width (in km)
"""
class PEERScalingRelation(BaseScalingRelation):
"""
A simple magnitude-area scaling relation to use as a placeholder until
updated by the synthetic rupture generator
"""
def get_rupture_dimensions(self, magnitude: float, rake: float = 0.0, dip: float = 90.0,
aspect: float = 1.0, thickness: float = 20.0,
epsilon: float = 0.0) -> Tuple[float, float, float]:
"""
Returns the area (km^2), length (km) and width (km) of the rupture using the PEER
Magnitude-Area scaling relation and constrained by the crustal thickness
"""
# PEER defined magnitude-area scaling relation returns area in km^2
area = self.get_area(magnitude)
width = sqrt(area / aspect)
# dz is the vertical dimension of the rupture (in km)
dz = width * sin(radians(dip))
if dz > thickness:
# Rupture is wider than the crustal thickness, so rescale it to the
# maximum width and break the aspect ratio
width = thickness / sin(radians(dip))
length = area / width
return area, length, width
@staticmethod
def get_area(magnitude: float) -> float:
"""
Returns the area from the PEER magnitude-area scaling relation
A = 10.0 ^ (magnitude - 4.0)
"""
return 10.0 ** (magnitude - 4.0)
class StrasserEtAl2010Interface(PEERScalingRelation):
"""
Implements the magnitude-area scaling relation for subduction interface earthquakes
from Strasser et al. (2010)
Strasser FO, Arango MC, Bommer, JJ (2010) "Scaling of the Source Dimensions of
Interface and Intraslab Subduction-zone Earthquakes with Moment Magnitude", Seismological
Research Letters, 81: 941 - 950, doi:10.1785/gssrl.81.6.941
"""
@staticmethod
def get_area(magnitude: float) -> float:
"""
Returns the median area from the magnitude-area scaling relation
"""
return 10.0 ** (-3.476 + 0.952 * magnitude)
class StrasserEtAl2010Inslab(PEERScalingRelation):
"""
Implements the magnitude-area scaling relation for subduction in-slab earthquakes
from Strasser et al. (2020)
Strasser FO, Arango MC, Bommer, JJ (2010) "Scaling of the Source Dimensions of
Interface and Intraslab Subduction-zone Earthquakes with Moment Magnitude", Seismological
Research Letters, 81: 941 - 950, doi:10.1785/gssrl.81.6.941
"""
@staticmethod
def get_area(magnitude: float) -> float:
"""
Returns the median area from the magnitude-area scaling relation
"""
return 10.0 ** (-3.225 + 0.89 * magnitude)
class Stafford2014(BaseScalingRelation):
"""
Implements the hazard-consistent magnitude scaling relation described by Stafford (2014)
Stafford PJ (2014) "Source-Scaling Relationships for the Simulation of Rupture
Geometry within Probabilistic Seismic-Hazard Analysis", Bulletin of the Seismological
Society of America, 104(4): 1620 - 1635, doi: 10.1785/012013024
"""
# Model coefficients for the specific style of faulting
COEFFS = {
# Coefficients from Table 1
1: {"U": {"a0": -27.4922, "a1": 4.6656, "a2": -0.2033},
"SS": {"a0": -30.8395, "a1": 5.4184, "a2": -0.3044},
"N": {"a0": -36.9770, "a1": 6.3070, "a2": -0.1696},
"R": {"a0": -35.8239, "a1": 5.0680, "a2": -0.0457}},
# Coefficients from Table 2
2: {"U": {"b0": -2.300, "b1": 0.7167, "sigma": 0.2337},
"SS": {"b0": -2.300, "b1": 0.7167, "sigma": 0.2337},
"N": {"b0": -4.1055, "b1": 1.0370, "sigma": 0.2509},
"R": {"b0": -3.8300, "b1": 0.9982, "sigma": 0.2285}},
# Coefficients from Table 3
3: {"U": {"gamma0": -9.3137, "sigma": 0.3138, "rho": 0.3104},
"SS": {"gamma0": -9.3137, "sigma": 0.3138, "rho": 0.3104},
"N": {"gamma0": -9.2483, "sigma": 0.3454, "rho": 0.4336},
"R": {"gamma0": -9.2749, "sigma": 0.2534, "rho": 0.1376}},
# Coefficients from Table 4
4: {"U": 0.7574, "SS": 0.7574, "N": 0.8490, "R": 0.7496}
}
def get_rupture_dimensions(self, magnitude: float, rake: float = 0.0, dip: float = 90.0,
aspect: float = 1.0, thickness: float = 20.0,
epsilon: float = 0.0) -> Tuple[float, float, float]:
"""
Gets the rupture dimensions from for the given magnitude subject to the physical
constaints.
Args:
epsilon:
Number of standard deviations above or below the median to determine the
rupture width
"""
sof = self._get_sof(rake)
# Get mean and standard deviation of rupture width
mu_rw, sigma_rw, max_width, p_i = self.get_rupture_width(magnitude, dip, sof, thickness)
# Get mean and standard deviation of rupture area
mu_ra, sigma_ra = self.get_rupture_area(magnitude, sof, max_width, sigma_rw)
# Apply censoring of the rupture width by the seismogenic thickness
F_rw_max_norm = norm.cdf(log(max_width), loc=mu_rw, scale=sigma_rw)
ncdf_epsilon = norm.cdf(epsilon)
target = ncdf_epsilon / F_rw_max_norm
if target > 1:
target = ncdf_epsilon
epsilon_rw = norm.ppf(target)
# Retrieve the rupture width and the rupture area conditioned on the rupture width
width = exp(mu_rw + epsilon_rw * sigma_rw)
if width > max_width:
width = max_width
epsilon_ra = self.COEFFS[4][sof] * epsilon
area = exp(mu_ra + epsilon_ra * sigma_ra)
length = area / width
return area, length, width
@staticmethod
def _get_sof(rake: float) -> str:
"""
Determine the style of faulting: strike slip (SS), reverse (R), normal (N) or
unknown (U)
Args:
rake: Rake of the rupture (in degrees from -180 to 180)
Returns:
Style of faulting
"""
if rake is None:
return "U"
if (-45 <= rake <= 45) or (rake >= 135) or (rake <= -135):
# strike slip
return "SS"
elif rake > 0:
# thrust/reverse
return "R"
else:
# normal
return "N"
def get_rupture_width(self, magnitude: float, dip: float, sof: str, thickness: float) ->\
Tuple[float, float, float, float]:
"""
Returns the parameters needed to define the censored rupture width
distribution defined in equations 8 to 14
Args:
magnitude: magnitude of earthquake
dip: Dip of earthquake in degrees from 0 to 90.0
sof: Style-of-faulting class (as string)
thickness: Seismogenic thickness
"""
# Gets the probability of a full width rupture
rw_max = thickness / sin(radians(dip))
z_i = self.COEFFS[1][sof]["a0"] + \
self.COEFFS[1][sof]["a1"] * magnitude +\
self.COEFFS[1][sof]["a2"] * rw_max
# Probability of rupturing full seismogenic thickness from logistic regression
p_i = 1.0 / (1.0 + exp(-z_i))
# Median width from equation 16
ln_rw = self.COEFFS[2][sof]["b0"] + \
self.COEFFS[2][sof]["b1"] * magnitude
# Equations 18 - 20
phi_rw = (log(rw_max) - ln_rw) / self.COEFFS[2][sof]["sigma"]
phi_rw_ncdf = norm.cdf(phi_rw)
ln_rw_trunc = ln_rw - self.COEFFS[2][sof]["sigma"] *\
(norm.pdf(phi_rw) / phi_rw_ncdf)
mean_rw = p_i * log(rw_max) + (1.0 - p_i) * ln_rw_trunc
# Equations 21 - 22
stddev_rw = self._get_rupture_width_sigma(self.COEFFS[2][sof]["sigma"],
phi_rw,
phi_rw_ncdf,
p_i)
return mean_rw, stddev_rw, rw_max, p_i
@staticmethod
def _get_rupture_width_sigma(sigma: float, phi_rw: float, phi_rw_ncdf: float,
p_i: float) -> float:
"""
Returns the variabiliy in the rupture width described by equations 21 and 22
"""
denom = sqrt(2.0 * pi) * phi_rw_ncdf
if phi_rw_ncdf >= 0.0:
elem1 = sqrt(pi / 2.) * (1.0 + chi2.cdf(phi_rw, 3))
else:
elem1 = sqrt(pi / 2.) * (1.0 - chi2.cdf(phi_rw, 3))
elem2 = exp(-(phi_rw ** 2)) / denom
sigma_truncated = sqrt(((sigma ** 2.) / denom) * (elem1 - elem2))
return (1.0 - p_i) * sigma_truncated
def get_rupture_area(self, magnitude: float, sof: str, rw_max: float,
sigma_lnrw: float) -> Tuple[float, float]:
"""
Returns the rupture area conditioned upon the maximum rupture width and style of
faulting
Args:
magnitude: magnitude of earthquake
sof: Style-of-faulting class (as string)
rw_max: Maximum rupture width (in km)
sigma_lnrw: Standard deviation of the logarithmic rupture width
Returns:
ln_ra: Natural logarithm of rupture area
sigma_ra: Standard deviation of the natural logarithm of rupture area
"""
mw_crit = (log(rw_max) - self.COEFFS[2][sof]["b0"]) /\
self.COEFFS[2][sof]["b1"]
ln_ra = self.COEFFS[3][sof]["gamma0"] + log(10.) * magnitude
if magnitude > mw_crit:
# Equation 23
ln_ra -= ((log(10.) / 4.) * (magnitude - mw_crit))
# Get the sigma log rupture area (equation 28)
sigma = self.COEFFS[3][sof]["sigma"]
sigma_ra = sqrt(
(sigma ** 2.) + (sigma_lnrw ** 2.) +
(2.0 * self.COEFFS[3][sof]["rho"] * sigma_lnrw * sigma))
return ln_ra, sigma_ra
......@@ -4,6 +4,8 @@ consistent quantities
"""
from typing import Tuple, Optional, Union
from openquake.hazardlib.geo.nodalplane import NodalPlane
from shakyground2.magnitude_scaling_relations import (BaseScalingRelation,
PEERScalingRelation)
def longitude(lon: float) -> float:
......@@ -108,3 +110,17 @@ def hypocenter_position(hypo_pos: Tuple[float, float]) -> Tuple[float, float]:
if down_dip < 0.0 or down_dip > 1.0:
raise ValueError("Down dip position %.3f should be in the range 0 to 1" % down_dip)
return along_strike, down_dip
def scaling_relation(msr: Optional[BaseScalingRelation]):
"""
Verifies that the magnitude scaling relation is one supported by the
software, or return a default is none is provided
"""
if not msr:
# If no relation is defined then use the default
return PEERScalingRelation()
if not isinstance(msr, BaseScalingRelation):
raise TypeError("Magnitude Scaling Relation %s not instance of BaseScalingRelation" %
str(msr))
return msr
......@@ -6,42 +6,8 @@ import numpy as np
from datetime import date, time
from math import radians, sin, pi
from openquake.hazardlib.geo import Point
from shakyground2.earthquake import DummyScalingRelation, Earthquake
class DummyScalingRelationTestCase(unittest.TestCase):
"""
Simple tests for the dummy magnitude scaling relation
"""
def setUp(self):
self.msr = DummyScalingRelation()
def test_simple_case_within_thickness(self):
area, length, width = self.msr.get_rupture_dimensions(6.0)
self.assertAlmostEqual(area, 100.0)
self.assertAlmostEqual(length, 10.0)
self.assertAlmostEqual(width, 10.0)
self.assertAlmostEqual(length * width, area)
def test_aspect_greater_than_one(self):
area, length, width = self.msr.get_rupture_dimensions(6.0, aspect=2.0)
self.assertAlmostEqual(area, 100.0)
self.assertAlmostEqual(length / width, 2.0)
self.assertAlmostEqual(length * width, area)
def test_rupture_rescale(self):
area, length, width = self.msr.get_rupture_dimensions(6.0, thickness=5.0)
self.assertAlmostEqual(area, 100.0)
self.assertAlmostEqual(length, 20.0)
self.assertAlmostEqual(width, 5.0)
self.assertAlmostEqual(length / width, 4.0)
self.assertAlmostEqual(length * width, area)
def test_dipping_rupture_rescaled(self):
area, length, width = self.msr.get_rupture_dimensions(7.0, dip=60.0, thickness=15.0)
self.assertAlmostEqual(area, 1000.0)
self.assertAlmostEqual(length * width, area)
self.assertAlmostEqual(width, 15.0 / sin(radians(60.0)))
from shakyground2.earthquake import Earthquake
from shakyground2.magnitude_scaling_relations import PEERScalingRelation
class EarthquakeTestCase(unittest.TestCase):
......@@ -75,8 +41,8 @@ class EarthquakeTestCase(unittest.TestCase):
dip = 90.0
usd = 0.0
lsd = 20.0
area, length, width = DummyScalingRelation.get_rupture_dimensions(mag, dip=dip,
thickness=lsd - usd)
area, length, width = PEERScalingRelation().get_rupture_dimensions(mag, dip=dip,
thickness=lsd - usd)
plane0 = Earthquake.build_planar_surface(centroid, strike, dip,
length, width, lsd, usd)
# Length should be approximately 10 km, but allow a slight precision difference
......@@ -100,8 +66,8 @@ class EarthquakeTestCase(unittest.TestCase):
dip = 90.0
usd = 0.0
lsd = 20.0
area, length, width = DummyScalingRelation.get_rupture_dimensions(mag, dip=dip,
thickness=lsd - usd)
area, length, width = PEERScalingRelation().get_rupture_dimensions(mag, dip=dip,
thickness=lsd - usd)
plane0 = Earthquake.build_planar_surface(centroid, strike, dip, length, width,
lsd, usd)
......@@ -122,8 +88,8 @@ class EarthquakeTestCase(unittest.TestCase):
dip = 90.0
usd = 0.0
lsd = 10.0
area, length, width = DummyScalingRelation.get_rupture_dimensions(mag, dip=dip,
thickness=lsd - usd)
area, length, width = PEERScalingRelation().get_rupture_dimensions(mag, dip=dip,
thickness=lsd - usd)
plane0 = Earthquake.build_planar_surface(centroid, strike, dip, length, width,
lsd, usd)
......@@ -144,8 +110,8 @@ class EarthquakeTestCase(unittest.TestCase):
dip = 90.0
usd = 0.0
lsd = 10.0
area, length, width = DummyScalingRelation.get_rupture_dimensions(mag, dip=dip,
thickness=lsd - usd)
area, length, width = PEERScalingRelation().get_rupture_dimensions(mag, dip=dip,
thickness=lsd - usd)
plane0 = Earthquake.build_planar_surface(centroid, strike, dip, length, width,
lsd, usd)
......@@ -167,8 +133,8 @@ class EarthquakeTestCase(unittest.TestCase):
dip = 60.0
usd = 0.0
lsd = 10.0
area, length, width = DummyScalingRelation.get_rupture_dimensions(mag, dip=dip,
thickness=lsd - usd)
area, length, width = PEERScalingRelation().get_rupture_dimensions(mag, dip=dip,
thickness=lsd - usd)
plane0 = Earthquake.build_planar_surface(centroid, strike, dip, length, width,
lsd, usd)
self.assertAlmostEqual(width, 10.0 / sin(radians(dip)))
......@@ -188,8 +154,8 @@ class EarthquakeTestCase(unittest.TestCase):
strike = 0.0
dip = 90.0
# Create a simple surface
area, length, width = DummyScalingRelation.get_rupture_dimensions(self.mag, dip=dip,
thickness=lsd - usd)
area, length, width = PEERScalingRelation().get_rupture_dimensions(self.mag, dip=dip,
thickness=lsd - usd)
plane0 = Earthquake.build_planar_surface(centroid, strike, dip, length, width,
lsd, usd)
......
import unittest
from math import radians, sin
from shakyground2.magnitude_scaling_relations import (PEERScalingRelation,
StrasserEtAl2010Interface,
StrasserEtAl2010Inslab,
Stafford2014)
class PEERScalingRelationTestCase(unittest.TestCase):
"""
Simple tests for the PEER scaling relation, including rescaling according to the
seismogenic thickness
"""
def setUp(self):
self.msr = PEERScalingRelation()
def test_simple_case_within_thickness(self):
area, length, width = self.msr.get_rupture_dimensions(6.0)
self.assertAlmostEqual(area, 100.0)
self.assertAlmostEqual(length, 10.0)
self.assertAlmostEqual(width, 10.0)
self.assertAlmostEqual(length * width, area)
def test_aspect_greater_than_one(self):
area, length, width = self.msr.get_rupture_dimensions(6.0, aspect=2.0)
self.assertAlmostEqual(area, 100.0)
self.assertAlmostEqual(length / width, 2.0)
self.assertAlmostEqual(length * width, area)
def test_rupture_rescale(self):
area, length, width = self.msr.get_rupture_dimensions(6.0, thickness=5.0)
self.assertAlmostEqual(area, 100.0)
self.assertAlmostEqual(length, 20.0)
self.assertAlmostEqual(width, 5.0)
self.assertAlmostEqual(length / width, 4.0)
self.assertAlmostEqual(length * width, area)
def test_dipping_rupture_rescaled(self):
area, length, width = self.msr.get_rupture_dimensions(7.0, dip=60.0, thickness=15.0)
self.assertAlmostEqual(area, 1000.0)
self.assertAlmostEqual(length * width, area)
self.assertAlmostEqual(width, 15.0 / sin(radians(60.0)))
class StrasserEtAl2010InterfaceTestCase(unittest.TestCase):
def test_strasser_interface_msr(self):
magnitudes = [4.0, 5.0, 6.0, 7.0, 8.0]
msr = StrasserEtAl2010Interface()
# Target values from OpenQuake implementation
targets = [2.147830474, 19.2309173, 172.1868575, 1541.70045295, 13803.8426460]
for mag, target in zip(magnitudes, targets):
self.assertAlmostEqual(msr.get_area(mag), target, 4)
class StrasserEtAl2010InslabTestCase(unittest.TestCase):