**2. Methods**

*2.1. Maximum Likelihood Estimation of the Parameters*

The Gutenberg-Richter distribution and all its derivations were originally developed using magnitudes. If we use seismic moments (*Mom*) instead of magnitudes, the Gutenberg-Richter distributions (unbounded/truncated/tapered) correspond to the Pareto distributions defined in Kagan [3], with slope parameter *β* equal to 2/3 the b-value.

The probability density function of the tapered distribution is [3]:

$$f(Mon) \quad = \left(\frac{\beta}{Mon} + \frac{1}{C\_{Mon}}\right) \left(\frac{Mon\_{\text{min}}}{Mon}\right)^{\beta} \exp\left(\frac{Mon\_{\text{min}} - Mon}{C\_{Mon}}\right) \tag{1}$$
 
$$\text{for } Mon\_{\text{min}} \le Mon < \infty$$

where *Mommin* is the seismic moment of completeness of the catalog, *β* is the parameter controlling the slope of the distribution, and *CMom* is the corner moment that controls the tail of the distribution. We stress that it is always possible to pass from the seismic moment to the magnitude distribution (here, we adopt the relationship defined in Kanamori [17]. In this case the corner moment *CMom* is called "corner magnitude" (*CM*), and the seismic moment of completeness corresponds to the magnitude of completeness.

If we have a seismic moment of completeness that varies with time (*Mom*(*i*) *min*, Figure 2), we can easily rewrite Equation (1) with:

$$f\_{\left(i\right)}\left(\text{Mom}\right) \\ \quad = \left(\frac{\beta}{\text{Mom}} + \frac{1}{\mathbb{C}\_{\text{Mom}}}\right) \left(\frac{\text{Mom}\_{\text{min}}^{\left(i\right)}}{\text{Mom}}\right)^{\beta} \exp\left(\frac{\text{Mom}\_{\text{min}}^{\left(i\right)} - \text{Mom}}{\mathbb{C}\_{\text{Mom}}}\right) \\ \tag{2}$$
 
$$\text{for } \text{Mom}\_{\text{min}}^{\left(i\right)} \le \text{Mom} < \infty$$

This relationship allows referring each observation to the completeness that holds at the time of its occurrence: in this time frame, indeed, Equation (2) describes the statistical distribution that holds.

Being both the parameters of the distribution (*CMom* and *β*) in common to all these distributions (*Mom*(*i*) *min* is a parameter related to the seismic catalog, estimated independently), their likelihood holds for all such distributions. Thus, if we have a seismic catalog with

*N* earthquakes with moments *x*1, ... , *xi*, ... , *xN*, the log-likelihood of the tapered Pareto distribution becomes:

$$LL(\mathbf{x}\_1, \dots, \mathbf{x}\_N | \boldsymbol{\beta}, \mathbf{C}\_{\text{Mon}}) = \sum\_{i=1}^N \ln \left[ f\_{(i)}(\mathbf{x}\_i) \right] \tag{3}$$

In Equation (3) the probability density function *f*(*i*) depends on the seismic moment of completeness relative to the i-th earthquake. In Figure 3 we summarize the scheme of our methodology, applied to completeness thresholds relative to Figure 2: in this case the log-likelihood of Equation (3) is obtained by summing up the log-likelihoods relative to the three different thresholds of completeness. Notably, if the seismic moment of completeness is the same for all the events, Equation (3) becomes the classical log-likelihood for the tapered Pareto distribution [3].

**Figure 3.** Graphical representation of the log-likelihood computation scheme proposed in this paper in the case of three different magnitudes of completeness thresholds.

To maximize the likelihood of observations, and evaluate the maximum likelihood estimation (MLE) of the parameters *β* and *CMom*, we adopt a brute-force approach, that is, we evaluated the log-likelihood for many potential combinations of the parameters, covering the entire parameter space [12], This allows obtain the complete description of the log-likelihood function: the maximum of the function (*LL*max) is, by definition, the MLE of the parameters. Moreover, this approach allows evaluating also the shape of the log-likelihood function, which is particularly useful to assess the uncertainty associated with the parameters' estimation, as it will be shown in the next section.

We stress the simplicity of our approach: to move from Equation (1) to Equation (2) we only need to substitute *Mommin* with *Mom*(*i*) *min*, i.e., using the time-variable seismic moment of completeness instead of the fixed one. As the parameters of the distribution that we want to evaluate are in common to all periods, we can simply stack their likelihoods, passing to Equation (3). Noteworthy, this is based on the same principle exploited in Vere-Jones et al. [18], who used a standard log-likelihood function for a tapered Pareto distribution with one or more parameters that change with times (Vere-Jones et al. [18], Equation (19)), as also in that case, this is possible because each distribution holds at the time of the observation and the different likelihoods can be stacked by summing them (Figure 3).

To check the robustness of this approach, we test its performance by estimating the distribution parameters from synthetic catalogs, for which such parameters are known. To this end, we simulate, using the Taroni and Selva [11] toolbox (based on the Vere-Jones et al. [18] method), thousands of synthetic catalogs with different input parameters (*β*, *CMom*, and magnitude of completeness), obtaining a good agreemen<sup>t</sup> between the MLE of the parameters and the input parameters, as expected. The results are shown in Table 1. The goodness of the agreemen<sup>t</sup> should be evaluated based on the estimated uncertainty on the parameters. Thus, the results of this comparison are discussed in the next section.

**Table 1.** Input and estimated *β* and *CM* relative confidence region for thousands of simulated synthetic catalogs.


#### *2.2. Estimation of the Uncertainties*

To evaluate the uncertainties relative to the parameters' estimation, we use a widely applied method [3,7,8,16,19] based on asymptotic theory [20], sometimes called profilelikelihood confidence region estimation [21]. It states that, if we want to estimate the confidence region of the parameters (in our case, of the tapered Gutenberg-Richter), we have to "cut" the log-likelihood function at a fixed threshold, and then look at the contour plot of this cut. Different thresholds correspond to different confidence intervals. For example, to obtain a 95% confidence region, we have to look at the *LL* = *LLmax* − 2.995 threshold, where *LLmax* is the maximum of the log-likelihood [8]. Hereinafter, the confidence region that describes the uncertainties on the parameters' estimation will be represented by the contour plot of the selected threshold (see Figure 4 for an illustrative example).

This procedure is adopted to evaluate the goodness of the agreemen<sup>t</sup> between input and estimated parameters for the thousands of synthetic catalogs discussed in the previous paragraph. In particular, we verify that the input parameters are enclosed in the 95% confidence region for the estimated parameters about 95% of the simulations, obtaining a very good agreement. The results are shown in Table 1.

**Figure 4.** Contour plot of the bivariate log-likelihood function for the parameters of the Tapered Gutenberg-Richter distribution; different colors represent the different log-likelihood values, according to the color bar on the left; the black curve represents the 95% confidence region, and the black dot represents the maximum likelihood estimation MLE.
