Merge remote-tracking branch 'upstream/master'

shakyground2/earthquake.py

@@ -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

shakyground2/magnitude_scaling_relations.py

+"""
+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

shakyground2/valid.py

@@ -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

tests/earthquake_test.py

@@ -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)

tests/magnitude_scaling_relations_test.py

+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):