rseismNet

Investigating the properties of the earthquake frequency-magnitude distribution (FMD) is the first necessary step of any seismicity analysis. The R package rseismNet provides a suite of statistical functions to determine the completeness magnitude $m_c$ and the seismicity properties below and above this threshold. “Net” refers to network, (1) for the seismic network properties described by the FMD shape below $m_c$, and (2) for the earthquake fault network properties described by the Gutenberg-Richter law above $m_c$. rseismNet also implements the Bayesian Magnitude of Completeness (BMC) method (Mignan et al., 2011, doi: 10.1785/0120100223), which maps $m_c$ based on a seismic network prior.

Basic functions

The following example makes use of all the basic functions of rseismNet, which are listed in R/fmd.R. It shows how to estimate the completeness magnitude $m_c$ from frequency magnitude distribution (FMDs) of different shapes, and how to avoid biased estimates of $\beta$. Accurate estimates of $\beta$ are important in seismic hazard assessment, as the probability of large and potentially damaging earthquakes is extrapolated from the value $\beta$ takes. A correct estimation of $m_c$ is also important in general seismicity statistics in order to maximize the sample size.

Let us first generate two earthquake magnitude vectors drawn from two different FMD models.

    n_eq <- 1e4
    theta_curved <- list(beta = log(10), mu = 2, sigma = 0.5)
    theta_angular <- list(kappa = 3 * log(10), beta = log(10), mc = 2)
    m_curved <- rseismNet::bfmd.sim(n_eq, theta_curved)
    m_angular <- rseismNet::efmd.sim(n_eq, theta_angular)

The vector m_curved is drawn from the FMD model first proposed by Ringdal (1975) and further developed by Ogata & Katsura (1993; 2006) by using the function bfmd.sim (“b” for bulk). The vector m_angular is drawn from the FMD model proposed by Mignan (2012) by using the function efmd.sim (“e” for elemental). The difference between the two models becomes clear once their FMD (computed with the function fmd) is plotted:

mdistr_curved <- rseismNet::fmd(m_curved)
mdistr_angular <- rseismNet::fmd(m_angular)
plot(mdistr_curved$mi, mdistr_curved$Ni, log = "y", col = "grey", main = "curved FMD")
points(mdistr_curved$mi, mdistr_curved$ni)
plot(mdistr_angular$mi, mdistr_angular$Ni, log = "y", col = "grey", main = "angular FMD")
points(mdistr_angular$mi, mdistr_angular$ni)

The first model leads to a curved FMD shape in the log-lin space, while the second one leads to an angular shape. Both models are of the form $n(m) \propto exp(-\beta m)q(m)$, the product of the theoretical Gutenberg-Richter law with slope $\beta = b\log(10)$ (Gutenberg & Richter, 1944) and a detection function $q(m)$, which gives the probability of a given magnitude $m$ being observed. In the first case, $q(m) = \mathcal{N}(\mu,\,\sigma^2)$, the cumulative Normal distribution, and in the second, $q(m < m_c) = exp(\kappa(m-m_c))$, an exponential function with detection parameter $\kappa > \beta$ and $q(m \geq m_c) = 1$. To learn more about the FMD shape ontology and where a given model applies, read Mignan (2012) and Mignan & Chen (2016). For $q(m)$ defined as a Gamma distribution (not yet implemented in rseismNet), see Kijko & Smit (2017).

Then the completeness magnitude $m_c$ is defined as the magnitude $m$ at which $q(m)$ tends to 1 (i.e., when all earthquakes are detected). For an angular FMD, $m_c$ is simply the mode of the distribution. For a curved FMD, $m_c$ is ambiguous and represented as $\mu+n\sigma$, but in practice it is approximated as the minimum magnitude above which a straight line is observed in log-lin space (i.e., above which the Gutenberg-Richter law is valid). Let us now evaluate $m_c$ by using three different methods, as defined in the function mc.val:

mc_mode_curved <- rseismNet::mc.val(m_curved, "mode")
mc_mbass_curved <- rseismNet::mc.val(m_curved, "mbass")
mc_gft_curved <- rseismNet::mc.val(m_curved, "gft")
plot(mdistr_curved$mi, mdistr_curved$ni, log = "y", main = "mc of a curved FMD")
abline(v = c(mc_mode_curved, mc_mbass_curved, mc_gft_curved), col = c("orange", "red", "brown"), lty = c("solid", "dashed", "dotdash"))

mc_mode_angular <- rseismNet::mc.val(m_angular, "mode")
mc_mbass_angular <- rseismNet::mc.val(m_angular, "mbass")
mc_gft_angular <- rseismNet::mc.val(m_angular, "gft")
plot(mdistr_angular$mi, mdistr_angular$ni, log = "y", main = "mc of an angular FMD")
abline(v = c(mc_mode_angular, mc_mbass_angular, mc_gft_angular), col = c("orange", "red", "brown"), lty = c("solid", "dashed", "dotdash"))

As of now, three FMD-based $m_c$ estimation methods are available in the mc.val function:

mode calculates the mode of the $m$ vector;
mbass (“median-based analysis of the segment slope”) determines the main breakpoints of the earthquake FMD. $m_c$ is defined as the change point that corresponds to the smallest probability of making an error when rejecting the null-hypothesis in a Wilcoxon-Mann-Whitney test (Amorese, 2007).
gft estimates the goodness-of-fit between the cumulative number of earthquakes observed and predicted by the Gutenberg-Richter law. $m_c$ is defined as the lowest magnitude bin at which a fixed threshold $R$ is first met. $R$ is defined as a normalized absolute difference, fixed to 0.95. If the threshold is not reached, 0.90 is used. If again the threshold is not reached, the mode is used instead (Wiemer & Wyss, 2000).

Both the mode and mbass methods are non-parametric while gft depends on the fitting of the Gutenberg-Richter law. Let us now use the function beta.mle to calculate the Maximum Likelihood Estimate of $\beta$ (Aki, 1965) and check whether the obtained $m_c$ estimates are correct (knowing that theta$beta = log(10) was fixed for the simulations).

Since all methods here give the correct $m_c$ estimate for the angular FMD case (with theta$mc = 2), we retrieve a reasonable estimate of the $\beta$-value, with:

beta_mode <- rseismNet::beta.mle(m_angular, mc_mode_angular)
beta_mode / log(10)
#> [1] 0.9808904
plot(mdistr_angular$mi, mdistr_angular$ni, log = "y", col = "grey", main = "Gutenberg-Richter law fit")
abline(v = mc_mode_angular, lty = "dotted", col = "red")
abline(a = log10(mdistr_angular$ni[which(mdistr_angular$mi >= mc_mode_angular)[1]]) + 
       beta_mode / log(10) * mc_mode_angular, b = -beta_mode / log(10), col = "red")

For a curved FMD however, the methods are likely to yield different $m_c$ estimates. The mode systematically underestimates it, leading in turn to an incorrect $\beta$-value:

beta_mode_curved <- rseismNet::beta.mle(m_curved, mc_mode_curved)
beta_mbass_curved <- rseismNet::beta.mle(m_curved, mc_mbass_curved)
beta_gft_curved <- rseismNet::beta.mle(m_curved, mc_gft_curved)
c(beta_mode_curved, beta_mbass_curved, beta_gft_curved) / log(10)
#> [1] 0.7106344 0.7908075 0.8872886
plot(mdistr_curved$mi, mdistr_curved$ni, log = "y", col = "grey", main = "Gutenberg-Richter law biases")
abline(v = c(mc_mode_curved, mc_mbass_curved, mc_gft_curved), col = c("orange", "red", "brown"), lty = c("solid", "dashed", "dotdash"))
abline(a = log10(mdistr_curved$ni[which(mdistr_curved$mi >= mc_mode_curved)[1]]) + 
       beta_mode_curved / log(10) * mc_mode_curved, b = -beta_mode_curved / log(10), 
       col = "orange")
abline(a = log10(mdistr_curved$ni[which(mdistr_curved$mi >= mc_mbass_curved)[1]]) + 
       beta_mbass_curved / log(10) * mc_mbass_curved, b = -beta_mbass_curved / log(10), 
       col = "red", lty = "dashed")
abline(a = log10(mdistr_curved$ni[which(mdistr_curved$mi >= mc_gft_curved)[1]]) + 
       beta_gft_curved / log(10) * mc_gft_curved, b = -beta_gft_curved / log(10), 
       col = "brown", lty = "dotdash")

rseismNet

Arnaud Mignan

2019-07-06

Installation

Basic functions