magnitude_scaling_relations.py 10.9 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
"""
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