Source code for openquake.hazardlib.mfd.youngs_coppersmith_1985

# coding: utf-8
# The Hazard Library
# Copyright (C) 2013-2023 GEM Foundation
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU Affero General Public License as
# published by the Free Software Foundation, either version 3 of the
# License, or (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU Affero General Public License for more details.
#
# You should have received a copy of the GNU Affero General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.
"""
Module :mod:`openquake.hazardlib.mfd.youngs_coppersmith_1985` defines the
Youngs and Coppersmith 1985 MFD.
"""
import numpy

from openquake.baselib.python3compat import round
from openquake.hazardlib.mfd.base import BaseMFD

# width of the boxcar function representing the characteristic
# distribution
DELTA_CHAR = 0.5


[docs]class YoungsCoppersmith1985MFD(BaseMFD): """ Class implementing the MFD for the 'Characteristic Earthquake Model' as described in: "Implications of fault slip rates and earthquake recurrence models to probabilistic seismic hazard estimates", by Robert R. Youngs and Kevin J. Coppersmith and published in Bulletin of the Seismological Society of America, Vol. 75, No. 4, pages 939-964, 1985. The class implements the MFD under the following assumptions as reported at page 954: 1) Δ_mc (width of the boxcar distribution representing characteristic rate) is equal to 0.5 magnitude unit 2) m' (maximum magnitude value for the Gutenberg-Richeter part of the distribution) is equal to the absolute maximum magnitude minus Δ_mc (that is there is no gap between the Gutenberg-Richter distribution and the boxcar distribution) 3) the rate of events at the characteristic magnitude is equal to the rate of events for magnitude equal to m' - 1 :param min_mag: The lowest possible magnitude for the MFD. The first bin in the :meth:`result histogram <get_annual_occurrence_rates>` is aligned to make its left border match this value. :param a_val: The Gutenberg-Richter ``a`` value -- the intercept of the loglinear cumulative G-R relationship. :param b_val: The Gutenberg-Richter ``b`` value -- the gradient of the loglinear G-R relationship. :param char_mag: The characteristic magnitude defining the middle point of the characteristic distribution. That is the boxcar function representing the characteristic distribution is defined in the range [char_mag - 0.25, char_mag + 0.25]. :param char_rate: The characteristic rate associated to the characteristic magnitude, to be distributed over the domain of the boxcar function representing the characteristic distribution (that is λ_char = char_rate / 0.5) :param bin_width: A positive float value -- the width of a single histogram bin. Values for ``min_mag`` and the maximum magnitude (char_mag + 0.25) don't have to be aligned with respect to ``bin_width``. They get rounded accordingly anyway so that both are divisible by ``bin_width`` just before converting a function to a histogram. See :meth:`_get_min_mag_and_num_bins`. """ MODIFICATIONS = set() def __init__(self, min_mag, b_val, char_mag, char_rate, bin_width, total_moment_rate=None): if total_moment_rate is not None: beta = b_val * numpy.log(10) mu = char_mag + DELTA_CHAR / 2 m0 = min_mag # seismic moment (in Nm) for the maximum magnitude c = 1.5 d = 9.05 mo_u = 10 ** (c * mu + d) # equations (16) and (17) solved for N(min_mag) and N(char_mag) c1 = numpy.exp(-beta * (mu - m0 - 0.5)) c2 = numpy.exp(-beta * (mu - m0 - 1.5)) c3 = beta * c2 / (2 * (1 - c1) + beta * c2) c4 = (b_val * (10 ** (-c / 2)) / (c - b_val)) + \ (b_val * numpy.exp(beta) * (1 - (10 ** (-c / 2))) / c) n_min_mag = (1 - c1) * total_moment_rate / ( (1 - c3) * c1 * mo_u * c4) char_rate = c3 * n_min_mag a_val = numpy.log10( (n_min_mag - char_rate) / (10 ** (- b_val * min_mag) - 10 ** (- b_val * (char_mag - 0.25)))) elif char_rate is not None: a_incr = b_val * (char_mag - 1.25) + numpy.log10( char_rate / DELTA_CHAR) a_val = a_incr - numpy.log10(b_val * numpy.log(10)) else: raise ValueError('Either `char_rate` or `total_moment_rate` must ' 'be not None') self.min_mag = min_mag self.a_val = a_val self.b_val = b_val self.char_mag = char_mag self.char_rate = char_rate self.bin_width = bin_width self.check_constraints()
[docs] def get_min_max_mag(self): "Return the minimum and maximum magnitudes" mag, num_bins = self._get_min_mag_and_num_bins() return mag, mag + self. bin_width * (num_bins - 1)
[docs] def check_constraints(self): """ Checks the following constraints: * minimum magnitude is positive. * ``b`` value is positive. * characteristic magnitude is positive * characteristic rate is positive * bin width is in the range (0, 0.5] to allow for at least one bin representing the characteristic distribution * characteristic magnitude minus 0.25 (that is the maximum magnitude of the G-R distribution) is greater than the minimum magnitude by at least one magnitude bin. * rate of events at the characteristic magnitude is equal to the rate of events for magnitude equal to m_prime - 1. This is done by asserting the equality (up to 7 digit precision) :: 10 ** (a_incr - b * (m' - 1)) == char_rate / 0.5 where ``a_incr`` is the incremental a value obtained from the cumulative a value using the following formula :: a_incr = a_val + log10(b_val * ln(10)) and ``m' - 1 = char_mag - 1.25`` """ if self.min_mag <= 0: raise ValueError('minimum magnitude must be positive') if self.b_val <= 0: raise ValueError('b value must be positive') if self.char_mag <= 0: raise ValueError('characteristic magnitude must be positive') if self.char_rate is not None and self.char_rate <= 0: raise ValueError('characteristic rate must be positive') if not 0 < self.bin_width <= DELTA_CHAR: err_msg = 'bin width must be in the range (0, %s] to allow for ' \ 'at least one magnitude bin representing the ' \ 'characteristic distribution' % DELTA_CHAR raise ValueError(err_msg) if not self.char_mag - DELTA_CHAR / 2 >= self.min_mag + self.bin_width: err_msg = 'Maximum magnitude of the G-R distribution (char_mag ' \ '- 0.25) must be greater than the minimum magnitude ' \ 'by at least one magnitude bin.' raise ValueError(err_msg) a_incr = self.a_val + numpy.log10(self.b_val * numpy.log(10)) actual = 10 ** (a_incr - self.b_val * (self.char_mag - 1.25)) desired = self.char_rate / DELTA_CHAR if not numpy.allclose(actual, desired, rtol=0.0, atol=1e-07): err_msg = 'Rate of events at the characteristic magnitude is ' \ 'not equal to the rate of events for magnitude equal ' \ 'to char_mag - 1.25' raise ValueError(err_msg)
[docs] @classmethod def from_total_moment_rate(cls, min_mag, b_val, char_mag, total_moment_rate, bin_width): """ Define Youngs and Coppersmith 1985 MFD by constraing cumulative a value and characteristic rate from total moment rate. The cumulative a value and characteristic rate are obtained by solving equations (16) and (17), page 954, for the cumulative rate of events with magnitude greater than the minimum magnitude - N(min_mag) - and the cumulative rate of characteristic earthquakes - N(char_mag). The difference ``N(min_mag) - N(char_mag)`` represents the rate of noncharacteristic, exponentially distributed earthquakes and is used to derive the cumulative a value by solving the following equation :: 10 ** (a_val - b_val * min_mag) - 10 ** (a_val - b_val * (char_mag - 0.25)) = N(min_mag) - N(char_mag) which can be written as :: a_val = log10(N(min_mag) - N(char_mag)) / (10 ** (- b_val * min_mag) - 10 ** (- b_val * (char_mag - 0.25)) In the calculation of N(min_mag) and N(char_mag), the Hanks and Kanamori (1979) formula :: M0 = 10 ** (1.5 * Mw + 9.05) is used to convert moment magnitude (Mw) to seismic moment (M0, Newton × m) :param min_mag: The lowest magnitude for the MFD. The first bin in the :meth:`result histogram <get_annual_occurrence_rates>` is aligned to make its left border match this value. :param b_val: The Gutenberg-Richter ``b`` value -- the gradient of the loglinear G-R relationship. :param char_mag: The characteristic magnitude defining the middle point of characteristic distribution. That is the boxcar function representing the characteristic distribution is defined in the range [char_mag - 0.25, char_mag + 0.25]. :param total_moment_rate: Total moment rate in N * m / year. :param bin_width: A positive float value -- the width of a single histogram bin. :returns: An instance of :class:`YoungsCoppersmith1985MFD`. Values for ``min_mag`` and the maximum magnitude (char_mag + 0.25) don't have to be aligned with respect to ``bin_width``. They get rounded accordingly anyway so that both are divisible by ``bin_width`` just before converting a function to a histogram. See :meth:`_get_min_mag_and_num_bins`. """ # there is some contorsion here to get the correct seismicity rates, # see https://github.com/gem/oq-engine/issues/4930 tmp = cls(min_mag, b_val, char_mag, None, bin_width, total_moment_rate) calculated_moment_rate = sum( [rate * 10. ** (1.5 * mag + 9.05) for (mag, rate) in tmp.get_annual_occurrence_rates()]) misfit = calculated_moment_rate / total_moment_rate total_moment_rate_adjusted = total_moment_rate / misfit return cls(min_mag, b_val, char_mag, None, bin_width, total_moment_rate_adjusted)
[docs] @classmethod def from_characteristic_rate(cls, min_mag, b_val, char_mag, char_rate, bin_width): """ Define Youngs and Coppersmith 1985 MFD by constraing cumulative a value from characteristic rate. The cumulative a value is obtained by making use of the property that the rate of events at m' - 1 must be equal to the rate at the characteristic magnitude, and therefore by first computing the incremental a value, using the following equation:: 10 ** (a_incr - b_val * (m_prime - 1)) == char_rate / 0.5 where ``m' - 1 = char_mag - 1.25``. The cumulative a value is then obtained as :: a_val = a_incr - log10(b_val * ln(10)) :param min_mag: The lowest magnitude for the MFD. The first bin in the :meth:`result histogram <get_annual_occurrence_rates>` is aligned to make its left border match this value. :param b_val: The Gutenberg-Richter ``b`` value -- the gradient of the loglinear G-R relationship. :param char_mag: The characteristic magnitude defining the middle point of characteristic distribution. That is the boxcar function representing the characteristic distribution is defined in the range [char_mag - 0.25, char_mag + 0.25]. :param char_rate: The characteristic rate associated to the characteristic magnitude, to be distributed over the domain of the boxcar function representing the characteristic distribution (that is λ_char = char_rate / 0.5) :param bin_width: A positive float value -- the width of a single histogram bin. :returns: An instance of :class:`YoungsCoppersmith1985MFD`. Values for ``min_mag`` and the maximum magnitude (char_mag + 0.25) don't have to be aligned with respect to ``bin_width``. They get rounded accordingly anyway so that both are divisible by ``bin_width`` just before converting a function to a histogram. See :meth:`_get_min_mag_and_num_bins`. """ return cls(min_mag, b_val, char_mag, char_rate, bin_width)
def _get_rate(self, mag): """ Calculate and return the annual occurrence rate for a specific bin. :param mag: Magnitude value corresponding to the center of the bin of interest. :returns: Float number, the annual occurrence rate for the :param mag value. """ mag_lo = mag - self.bin_width / 2.0 mag_hi = mag + self.bin_width / 2.0 if mag >= self.min_mag and mag < self.char_mag - DELTA_CHAR / 2: # return rate according to exponential distribution return (10 ** (self.a_val - self.b_val * mag_lo) - 10 ** (self.a_val - self.b_val * mag_hi)) else: # return characteristic rate (distributed over the characteristic # range) for the given bin width return (self.char_rate / DELTA_CHAR) * self.bin_width def _get_min_mag_and_num_bins(self): """ Estimate the number of bins in the histogram and return it along with the first bin center value. Rounds ``min_mag`` and ``max_mag`` with respect to ``bin_width`` to make the distance between them include integer number of bins. :returns: A tuple of 2 items: first bin center, and total number of bins. """ min_mag = round(self.min_mag / self.bin_width) * self.bin_width max_mag = (round((self.char_mag + DELTA_CHAR / 2) / self.bin_width) * self.bin_width) min_mag += self.bin_width / 2.0 max_mag -= self.bin_width / 2.0 # here we use math round on the result of division and not just # cast it to integer because for some magnitude values that can't # be represented as an IEEE 754 double precisely the result can # look like 7.999999999999 which would become 7 instead of 8 # being naively casted to int so we would lose the last bin. num_bins = int(round((max_mag - min_mag) / self.bin_width)) + 1 return min_mag, num_bins
[docs] def get_annual_occurrence_rates(self): """ Calculate and return the annual occurrence rates histogram. :returns: See :meth: `openquake.hazardlib.mfd.base.BaseMFD.get_annual_occurrence_rates`. """ mag, num_bins = self._get_min_mag_and_num_bins() rates = [] for i in range(num_bins): rate = self._get_rate(mag) rates.append((mag, rate)) mag += self.bin_width return rates