# Seismological Research Letters

- © 2015 by the Seismological Society of America

## ABSTRACT

Faults have been increasingly integrated into seismic‐hazard assessments. We have developed a package of MATLAB (http://www.mathworks.com/products/matlab, last accessed January 2016) tools (called *FiSH*), designed to help seismic‐hazard modelers analyze fault data. These tools enable the derivation of expected earthquake rates, given common fault data, and allow researchers to test the consistency between the magnitude–frequency distributions (MFDs) assigned to a fault and some available observations. The basic assumption of *FiSH* is that the geometric and kinematic features of a fault are the expression of its seismogenic potential. Three tools have been designed to integrate the variable levels of information available: (1) the first tool allows users to convert fault geometry and slip rates into a global budget of the seismic moment released in a given time frame, taking uncertainties into account; (2) the second tool computes the recurrence parameters and associated uncertainties from historical and/or paleoseismological data; and (3) the third tool outputs time‐independent or time‐dependent earthquake rates for different MFD models. We present a test case to illustrate the capabilities of *FiSH*, on the Paganica normal fault in central Italy that ruptured during the 2009 L’Aquila earthquake sequence (mainshock *M*_{w} 6.3). The source codes are open, and we encourage users to handle the scripts, communicate with us regarding bugs, and/or suggest further improvements. Our intent is to distribute these tools in order to help researchers to pinpoint potential inconsistencies and obtain reliable fault‐based seismic‐hazard evaluations.

## INTRODUCTION

Probabilistic seismic hazard assessment (PSHA) has been used extensively by the scientific community over the last two decades and has led to the distribution of many seismic‐hazard analysis tools (e.g., Bender and Perkins, 1987; Ordaz *et al.*, 2013; Pagani *et al.*, 2014). The majority of these tools have been conceived to customize the representation of the results of a given PSHA model in terms of maps, hazard curves and spectra, and disaggregation analyses (e.g., for Italy, http://esse1-gis.mi.ingv.it/s1_en.php, last accessed December 2015; Europe, http://www.efehr.org:8080/jetspeed/portal/hazard.psml, last accessed December 2015; and the United States of America, http://earthquake.usgs.gov/hazards/apps/, last accessed December 2015). Several applications have focused on the visualization and/or analysis of individual seismic‐hazard components, such as the ground‐motion prediction equation or the magnitude–frequency distribution (MFD) of seismic sources (e.g., Wiemer, 2001; Field *et al.*, 2003; http://www.opensha.org/apps, last accessed December 2015). Tools that are solely dedicated to converting geological field data on faults into seismic‐hazard model components (such as source geometries and earthquake rupture forecasts) are rare. In general, fault source modeling uses simplistic rules based on the size and slip rate of mapped faults to estimate seismicity rates by adopting *a priori* MFDs (of the Gutenberg–Richter type, following the characteristic earthquake model, or by hybrid formulations).

The practice of using fault sources is growing in popularity, including in regions with moderate seismic activity, such as the European countries. In these areas, fault identification may be affected by similarly large uncertainties in the historical and instrumental seismic histories of more active areas that have not been inhabited for long period time. Certain studies have effectively applied a time‐dependent perspective to combine historical and instrumental seismic data with geological and paleoseismological information, partially compensating for a lack of information (e.g., in Italy, Pace *et al.*, 2006, and Peruzza *et al.*, 2011; in New Zealand, Stirling *et al.*, 2012; in California, Field *et al.*, 2014, 2015).

*FiSH* is a graphical user interface (GUI)‐supported tool that is designed to help seismic‐hazard modelers analyze fault data. *FiSH* is a package of MATLAB routines that is designed to evaluate expected earthquake rates, given certain fault data constraints, and to test the consistency between the MFD that is assigned to a fault and any available observations. The *FiSH* output can be used as input in other software packages that use faults to evaluate seismic hazards, particularly the software packages OpenQuake (Pagani *et al.*, 2014), which was developed for the Global Earthquake Model, and CRISIS2008 (Ordaz *et al.*, 2013).

The initial concept of *FiSH* originated in the early 2000s because of the encountered difficulties with using faults in PSHA for central Italy (Barchi *et al.*, 2000; Peruzza and Pace, 2002; Pace *et al.*, 2006). Since then, new theoretical formulations have been developed (Peruzza *et al.*, 2010), and the first MATLAB routines were released in 2013 through an Italian project (Italian Presidenza del Consiglio dei Ministri—Dipartimento della Protezione Civile [INGV‐DPC], Project S2‐2012, in Pace *et al.*, 2013). We tested and implemented these tools over the past two years, and this article represents the official release of *FiSH* v.1.02. We illustrate the capabilities of *FiSH* with a test case that involves the Paganica fault, the causative fault of the L’Aquila earthquake in 2009 that led to 309 casualties (Boncio *et al.*, 2010). The example shows that we can derive time‐independent and time‐dependent earthquake rates that are mutually consistent and hopefully more realistic than what would be obtained with simplistic rules by considering (1) field geological data on fault geometry and slip rates, (2) earthquake series from historical and paleoseismological studies, and (3) the year of the last known major earthquake to have occurred on the fault.

## GENERAL FEATURES

The basic assumption of *FiSH* is that the geometric and kinematic features of a fault are the expression of its seismogenic potential. Defining the seismogenic sources that underlie the use of *FiSH* is important. Taking advantage of the segmentation model defined by Boncio *et al.* (2004), a “seismogenic master fault” in this work is a major structure that can be considered substantially continuous at depth for several kilometers. Seismogenic master faults are separated from each other by first‐order structural or geometrical complexities (Fig. 1). For central Italy, Boncio *et al.* (2004) defined these complexities as the following barriers: (1) 3–4 km fault gaps among aligned structures; (2) sharp bends or intersections with cross structures that are 4–9 km long (often transfer faults) and oriented at nearly right angles to the intersected faults; and (3) overlapping or underlapping en echelon arrangements with separations between 2 and 5 km faults. More recently, Field *et al.* (2015) used a list of plausibility criteria to determine which subsections of a fault can interact to produce an earthquake, including distance and azimuthal changes between fault sections, cumulative azimuthal and rake changes, and the Coulomb criterion.

Based on this framework, a seismogenic source is a fault that can break entirely during an earthquake. As shown in Figure 1, the shape and size of the seismogenic source must be evaluated in three dimensions by defining the depth extent of the fault, the dip angle, and the length along strike. These parameters describe the geometry of a seismogenic source in *FiSH*. The fault length and seismogenic depth are assumed to be known parameters when using our tools; similarly, the crustal properties (e.g., shear modulus and strain drop) may be derived from specific studies or the general literature.

*FiSH* is written in MATLAB (v.R2011b), a Mathworks (http://www.mathworks.com, last accessed December 2015) commercial software package widely used by geosciences researchers. A MATLAB license is mandatory to run *FiSH*, but, due to *FiSH*’s GUI, no knowledge of the MATLAB language is required to use the tools. The source codes are available at http://fish-code.com (last accessed January 2016), and we welcome implementations or user feedback. *FiSH* has been tested using several operating systems (Macintosh, Unix, and PC), and it should run on all platforms that support MATLAB.

Three tools (workflow in Fig. 2), designed to integrate variable levels of information that is available for each fault, have been developed. The first tool allows users to convert common fault data, such as geometry and slip rates, into a global budget of the seismic moment that was released in a given time frame, taking uncertainties into account. The second tool computes the recurrence parameters (RP) and associated uncertainties from historical and/or paleoseismological data. The third tool outputs time‐independent or time‐dependent earthquake rates for different MFD models.

In particular, the codes that are currently provided with *FiSH* accomplish the following tasks:

Moment budget (MB) uses fault data, such as the style of faulting, geometry, and slip rate variability, to calculate the global budget of the seismic moment rate that is permitted by the structure based on predefined size–magnitude relationships. We calculate MB using two values: the maximum magnitude

*M*_{max}and the associated mean recurrence time*T*_{mean}, including their uncertainties in terms of the standard deviation of a Gaussian distribution for*M*_{max}and the coefficient of variation (CV, standard deviation of the recurrence times over their mean) for*T*_{mean}.RP uses the time series of past major earthquake occurrences (historical events assigned to a fault and/or paleoseismological data that include the uncertainty in dates). Using a Monte Carlo approach, RP calculates the distributions of recurrence times of major events on that fault and their variability.

Activity rates (AR) use the output from MB to calculate the Poisson and time‐dependent (if the elapsed time from the last earthquake is available) probabilities of earthquake occurrences and the expected annual rates based on different predefined MFD models. Currently,

*FiSH*supports four types of MFDs: a single magnitude value, a characteristic distribution with a Gaussian bell that is symmetric around a central magnitude, a classical Gutenberg–Richter type (Gutenberg and Richter, 1944) and a truncated Gutenberg–Richter type (Ordaz, 1999; Kagan, 2002).

## TOOLS BACKGROUND

In this section, we describe the *FiSH* routines and the formulas they use.

### MB Tool

MB combines standard data that were collected for faults or are assumed to be known, such as the size and rheological properties, with empirical size–magnitude relationships and allows for the treatment of their uncertainties. The input list that is defined by the user consists of the following: (1) the fault name, (2) a scale relationship code, (3) the length along strike (kilometers), (4) the dip angle (degrees), (5) the thickness of the seismogenic layer (kilometers), (6–7) the minimum and maximum slip rate (mm/yr), (8) the maximum observed magnitude (*M*_{w}, if any), (9) the standard deviation in magnitude of the observed event (if any), (10) the date of the last event (year, if any), (11) the shear modulus, and (12) the strain drop. The first 10 input data are mandatory, whereas the last two are optional and can be set homogeneously for all the faults in the GUI or entered into the input file (only the significand; exponents are set to 10^{−5} for strain drop and 10^{10} Pa for the shear modulus); fields 8–10 must be set as NaN (not a number) if they are unknown. The scale relationship code (see information in the MB interface in the “Scale Relationship” format button) allows the program to use the proper empirical relationship to transform the fault size into the maximum magnitude, that is, the magnitude that corresponds to the complete rupture of the entire fault plane. Currently, four relationships for tectonic and volcanic contexts have been implemented (Wells and Coppersmith, 1994, and Leonard, 2010, for tectonic contexts; Villamor *et al.*, 2001, and Azzaro *et al.*, 2015, for volcanic contexts); the scale relationship code assigns the chosen relationship(s) and the associated set of coefficients to the fault.

The user has to prepare an input text file (Fig. 3) with fields that are separated by blanks or commas and contain one row for each fault to be analyzed. Next, MB prompts the user for the output filename and for other choices; the shear modulus *μ* and strain drop *k* (defined as the displacement‐to‐length ratio, *D*/*L*) are values that are set as homogeneous by default because they are commonly assumed to be constant in a seismotectonic region (Scholz, 1990). However, these values can be also set differently for each fault (e.g., for the sensitivity testing of their impact on the hazard input). These values are present in the scalar seismic moment (*M*_{0}) formula as follows: (1)in which *L*, *W*, and *D* are, respectively, the along‐strike fault length, the down‐dip width, and the average displacement of a rectangular source (Fig. 1).

For each fault, depending on the input choice, up to five *M*_{max} values (and their errors) are computed:

A magnitude value

*MM*_{0}is directly computed by applying the standard formula*M*_{w}= 2/3 (log*M*_{0}−9.1), which was accepted by International Association of Seismology and Physics of the Earth’s Interior (IASPEI, 2005), to equation (1) to avoid the occasional rounding errors that can occur in the formula that was originally published by Hanks and Kanamori (1979).A value of magnitude (

*M*_{ASP}, ASP for aspect ratio) is computed by modifying the along‐strike dimension if it exceeds the length that is predicted by the aspect ratio relationships, which were derived by Peruzza and Pace (2002) from a slightly modified Wells and Coppersmith (1994) dataset, followed by the application of equation (1).Two magnitude values that depend on the choice of the scale relationship are calculated. For example, if the user chooses to use Wells and Coppersmith (1994), the two magnitudes are given by empirical relationships that were established for either the maximum subsurface fault rupture length (MRLD) or maximum rupture area (MRA).

A value that corresponds to the maximum observed magnitude (

*M*_{obs}), if available in the input file, is calculated.

Because all the empirical relationships and observations are affected by uncertainties, MB is designed to take these factors into account and return a maximum magnitude value and a standard deviation. The uncertainties in the empirical scaling relationship are taken from the studies of Wells and Coppersmith (1994), Villamor *et al.* (2001), Peruzza and Pace (2002), Leonard (2010), and Azzaro *et al.* (2015). Currently, the uncertainty in *MM*_{0} is fixed and set to 0.3, whereas the uncertainty in *M*_{obs} is assigned by the user. To combine the maximum magnitudes, MB draws a probability curve for each magnitude estimate by assuming a normal distribution. The user can define the number of standard deviations (*σ*) for truncating the normal distribution of magnitudes at both sides (left empty for unbounded distributions). MB successively sums the probability density curves and fits the summed curve to a normal distribution to obtain the mean of the maximum magnitude *M*_{max} and its standard deviation *σM*_{max}. MB also allows the introduction of weights to create the probability density curves of the magnitudes. Therefore, *M*_{max} represents an evaluation of the maximum rupture that is allowed by the fault geometry and the rheological properties. An output figure of this process is given for each fault (an example that is described later is given in Fig. 4), so the user can easily compare the distributions and the impact of different magnitudes on the final value.

To obtain the mean recurrence time of *M*_{max}, hereafter referred as *T*_{mean}, MB uses the criterion of “segment seismic moment conservation” that was proposed by Field *et al.* (1999), which divides the seismic moment that corresponds to *M*_{max} by the moment rate given a slip rate (2)in which *T*_{mean} is the mean recurrence time in years, Char_Rate is the annual mean rate of occurrence, *M*_{max} is the previously computed mean maximum magnitude, *μ* is the shear modulus, *V* is the average long‐term slip rate (obtained by averaging the minimum and maximum slip rates defined by the user in the input file), and *L* and *W* are the geometrical parameters of the fault. The coefficients 1.5 and 9.1 are the standard values that are adopted by IASPEI (2005) for the relationship between the moment magnitude and seismic moment (in N·m).

To take into account the uncertainties in *M*_{max} and the slip rates and to explore how these uncertainties affect the variability in *T*_{mean}, we introduced formal error propagation following Peruzza *et al.* (2010). The general formula for error propagation is given by a Taylor series: (3)and the variance is given by (4)

The variance is applied to the *T*_{mean} obtained from equation (2) using partial derivatives of the slip rate *V* and magnitude *M*_{max}. The variance of *T*_{mean} becomes (5)in which *dM*_{max} and *dV* are generic small variations in the magnitude and slip rate. Here, they are substituted with the standard deviation *σM* (previously described) and *σV* for the analysis of error propagation. The parameter *σV* is assumed to be half of the user‐defined minimum and maximum slip rates (given in the input file).

The output of MB consists of a figure with curves of the magnitude distributions, as well as the derived mean and standard deviation (e.g., Fig. 4). The output also includes a file (Fig. 3b) that contains the following information for each fault: the identification number, the mean value of *M*_{max} and its standard deviation, *T*_{mean} (yr), the CV for *T*_{mean}, the time elapsed (yr) since the last earthquake (if given), the seismic moment rate (N·m·yr^{−1}, defined as the ratio between the seismic moment that corresponds to *M*_{max} and *T*_{mean}), and the name of the fault. This output can be used as the input for the AR tool.

### RP Tool

Recent work has focused on deriving statistical parameters, such as the mean and standard deviation of the recurrence times, from paleoseismological or historical earthquake time series by considering the uncertainties in dating (e.g., Parsons, 2008). The user should thus check them with the parameters that are derived from the moment budget to enlighten the potential limitations that are associated with paleoseismological and historical datasets, for example, incomplete catalogs or uncertain magnitudes of the causative events.

The RP tool is designed to take these uncertainties into account and show how they can influence the desired parameters. Starting from an input file that contains the youngest and oldest years of occurrence for each event in the series, RP produces *n* simulations of the earthquake catalog (hereafter referred to as “paleoeqs‐simulations”) using a uniform distribution for the occurrences within its window of uncertainty. This simplistic view does not consider the most recent developments in age modeling by paleoseismologists (e.g., the OxCal program, Bronk Ramsey and Lee [2013], available at http://c14.arch.ox.ac.uk/embed.php?File=oxcal.html#program, last accessed December 2015) but, conversely, enables the users to handle the literature data, which are usually provided as rough time intervals. Future implementations of *FiSH* could introduce arbitrary probability density distributions. For each paleoeqs‐simulation, RP estimates the arithmetic mean of the recurrence times and its standard deviation (we call this case “arithmetic”). RP also fits each paleoeqs‐simulation using three probability functions, assuming that the events follow a Brownian passage time (BPT), a Weibull, and a Poisson distribution.

For the BPT distribution (Matthews *et al.*, 2002), the probability density function is (6)in which *T*_{mean} is the mean recurrence time, *α* is the aperiodicity of the distribution, and *t* is the time in years. RP returns two values, *T*_{mean} and *α*.

For the Weibull distribution (Patel *et al.*, 1976), the probability density function is (7)in which *b*>0 is the shape parameter and *a*>0 is the scale parameter of the distribution.

The *T*_{mean} and variance (*σ*) of the Weibull distribution are provided by (8)(9)RP returns the values *a*, *b*, *T*_{mean}, and CV (CV=*σ*/*T*_{mean}).

In a Poisson distribution of the earthquakes, the times between events are described by an exponential distribution, for which the probability density function is (10)in which *λ* is the inverse of the mean recurrence time (*T*_{mean}) obtained by RP.

Because this fit is performed for each of the paleoeqs‐simulations, *n* output parameters are obtained. These values are shown as contour graphs (for the arithmetic, BPT, and Weibull distributions) and histograms of *T*_{mean} (for the Poisson distribution), see an example in Figure 6.

The input required by the RP is an index file (with a .txt extension, see Fig. 5) that lists the names of the files that contain the earthquake time series for each fault. These files (with a .paleo extension) should be prepared with one row for each event and two columns for the earliest and latest years of earthquake occurrences, thus accounting for the dating uncertainties; for historical events, the same year has to be repeated. Once the input file is uploaded, RP prompts the user to set the number of simulations *n* and the output name. The user can define a *seed* number in the GUI to control the random generator so that the user can reproduce the same random simulation, and the calculations are replicable using the same *seed* number. This structure enables the users to treat, for example, only some of the faults that were previously characterized in MB or to check the changes from different age modeling/interpretations on the same fault trench.

The outputs consist of a figure (e.g., Fig. 6) that includes sample data, contour graphs of the paleoeqs‐simulations according to the probability distributions and histogram for each row of the index file, and a summary file (Fig. 5) with the statistical parameters of all the rows.

The pairs of *T*_{mean} and CV that are obtained by RP can be entered into the AR tool depending on the user’s need.

### AR Tool

AR is designed to return annual earthquake rates. The user first enters the standard output of the MB tool or a new file from the RP results that contains the following for each row/fault: (1) the identification number, (2) *M*_{max}, (3) its standard deviation *σM*_{max}, (4) *T*_{mean}, (5) CV, (6) the seismic moment rate, (7) the elapsed time from the last assigned earthquake (if given), and (8) the name of the fault. All the fields must be filled, but not all the inputs will necessarily be used; the basic information for modeling Gutenberg–Richter (GR) MFDs, that is, the minimum magnitude and *b*‐value, is handled separately in a second input file. We hope that this choice will allow a user to easily perform several tests using a common input. Once the user has added the input file(s), AR prompts the user to name the output file, choose the MFD, set an observation window and magnitude bin, and flag the figure options. The standard output consists of two text files with the annual rates and probabilities in the time window that was selected in the GUI (e.g., Fig. 7); additional graphics are plotted for each row/fault (e.g., Fig. 8).

AR balances the annual earthquake rates over the range of magnitudes for the MFD with the moment rate. The details of this process for the various choices of MFDs are explained below.

#### Single‐Value Model (Poisson)

In the single‐value Poisson model, the MFD collapses into a single value that is given by (*M*_{max}, 1/*T*_{mean}) in the input file. Similar to the MB tool, the seismic moment rate is given by the seismic moment of *M*_{max} (using the standard conversion formula) divided by the *T*_{mean}. If this value differs from the sixth value that is listed in the AR input file (e.g., the MB output has been modified according to the RP results or by imposing the seismic moment rate, for example, by geodetic considerations), a warning message appears. If the option that forces the use of the input seismic moment rate is flagged, a new *T*_{mean} is computed; a warning message informs the user if this new recurrence time differs with respect to the original second and third input values. The same flow chart applies to the other choices that are described hereafter; thus, the code serves to check inconsistencies between geometrical and energetic formulations and to perform sensitivity testing on the seismicity rates.

#### Single‐Value Model (BPT)

In the single‐value BPT model, AR also uses the CV and the elapsed time of the input file. As in the previous case, the MFD collapses into a single value that is given by the *M*_{max} of the input file, but its annual rate is given in terms of a time‐dependent probability that is assigned to the fault using the BPT formula (equation 6).

The conditional probability (*P*|elap) that an event occurs during the next Δ*T* yr, given an elapsed time *T*_{elap} since the last event, is defined as follows: (11)

In this case, the *T*_{elap} is a mandatory input, Δ*T* is the observation window that is set in the GUI, and the CV of the mean recurrence time is the *α* value in equation (6). MB returns a CV value, which is defined as the ratio between *σT*_{mean} and *T*_{mean}, but the user can assign different values, for example, after analyzing the date sequence with RP. Similar to the previous (Poisson) case, a check on the seismic moment rate is performed.

Thus, using the simplification that was proposed by Wu *et al.* (1995), a fictitious recurrence time *T*_{fict} for *M*_{max} is computed by solving the equivalence of the probabilities as follows: (12)in which *P*|elap, which is the conditional probability from the BPT model (equation 11), is set equal to the probability of a Poisson process given a selected observation period Δ*T*. The inverse of *T*_{fict} is used to obtain the rate of occurrence at *M*_{max}.

The output files for each fault are similar to the previous case; the annual rates and probabilities do not refer to any time window, but only to the Δ*T* yr after the starting time, which is used to compute the elapsed time; forward retrospective testing can be performed at any starting time by modifying only the *T*_{elap} parameter in the input file. If the figure option is set in the GUI, an additional plot that represents the Poisson and BPT probabilities versus the elapsed time is given (see the example in Fig. 8, which is described in more detail in the next section).

#### Single‐Value Model (User-Defined Probability)

In the most generic case, AR can use probabilities that are defined by the user: with this option, the input file must also contain an additional column with the probability of the full rupture of the fault in the selected time window. A fictitious recurrence time *T*_{fict} for *M*_{max} is now computed by solving the equivalence of the probabilities in equation (10), where *P*|elap is the user‐defined probability in the input file. Incorrect results may derive from the uncontrolled utilization of the Δ*T* window in the GUI. The inverse of *T*_{fict} is used to obtain the rate of occurrence at *M*_{max}. This choice can be used, for example, if different probability density functions from BPT (e.g., Weibull, which is characterized by the RP tool) are used to model time dependency.

#### Characteristic Gaussian Model (Poisson)

In addition to the Single‐Value Model (Poisson) section, AR can use *M*_{max} with its standard deviation (*σM*_{max}). The MFD is provided by a Gaussian curve that is centered on the maximum magnitude in the input file and ranges from a lower bound *M*_{maxlb}=*M*_{max}−*σM*_{max} to an upper bound *M*_{maxub}=*M*_{max}+*σM*_{max} with a step that must be specified by the user in the GUI. The seismic moment rate is thus partitioned and balanced over this range of magnitudes. One output file lists *M*_{maxlb} in the user‐defined bin width and the annual rates for all the magnitude bins; the other lists the Poisson probability that at least one earthquake will occur with a magnitude greater than or equal to *M*_{maxlb}. The particular format of the first output file is dictated by the potential need for an output that can be easily used as input for the OpenQuake (Pagani *et al.*, 2014) and CRISIS2008 (Ordaz *et al.*, 2013) software when fault sources are used.

#### Characteristic Gaussian Model (BPT)

Similar to the previous case and the Single Value Model (User‐Defined Probability) section, AR can consider the BPT model to be valid under the assumption that an earthquake with a magnitude in the range *M*_{maxlb}–*M*_{maxub} allows the renewal process to begin again. Thus, the probability of equation (6) must be representative of an earthquake with a magnitude in the range *M*_{maxlb}–*M*_{maxub}. Therefore, *T*_{mean} from equation (6) is the inverse of the cumulative annual rate of earthquakes with *M*_{w}≥*M*_{maxlb}, and the fictitious recurrence time *T*_{fict} by equation (12) is computed similarly to that in the Single Value Model (User‐Defined Probability) section. The output is the same as in the Characteristic Gaussian Model (Poisson) section.

#### Characteristic Gaussian Model (User‐Defined Probability)

If the input file contains an additional column with the probability that is assigned to the full rupture, AR can compute seismicity rates and probabilities for the magnitude bins of a Gaussian bell. The logical flowchart is the same as in the previous case.

#### Classical GR

The identification number is used by AR when modeling seismicity rates with a GR distribution. An additional input file that lists the identification number, the lower threshold magnitude of the distribution (*M*_{t}), and the *b*‐value of the Gutenberg–Richter distribution is required. In this case, the distribution is truncated at both ends by *M*_{t} and *M*_{max}, as (13)in which *N*_{(M)} is the number of earthquakes with a magnitude greater than or equal to *M*_{t} and less than or equal to *M*_{max}, and *a* is the logarithm of the number of earthquakes with *M*≥0. AR computes *a* to balance the total seismic moment rate (given by summing the rates of moment magnitudes for the bins of magnitudes with a user‐defined step) with the seismic moment rate that was obtained by the pair *M*_{max} and *T*_{mean}.

The output consists of one file with seismicity rates that lists the *M*_{t}, the user‐defined bin, and the annual rates for each magnitude bin. Given the well‐known problems of decumulation, the use of this formulation has been progressively abandoned in PSHA.

#### Truncated GR

This case is similar to the previous one, and the additional input file is the same as in the classical GR case. The distribution is truncated at both ends by *M*_{t} and *M*_{max}, but the truncated GR is characterized by a smooth transition to *M*_{max} instead of the abrupt drop of the classical GR model. Following Kagan (2002), the cumulative density function of the truncated GR model is (15)in which *M*_{u} is evaluated by shifting the *M*_{max} of a magnitude bin to ensure a nonzero rate of earthquakes at *M*_{max}, and *β* is (2/3)*b*; this equivalence is handled within the code, and the input file must define the *b*‐value. As in classical GR, the seismic moment rate is partitioned and balanced over this range of magnitudes with a step that must be specified by the user in the GUI. The output is the same as in the previous case.

## PAGANICA FAULT CASE STUDY

The example shown below is intended to illustrate the capabilities of *FiSH*. The case study is taken from published articles, because presenting new data and/or interpretations is beyond the scope of this article.

The Paganica normal fault in central Italy ruptured during the 2009 L’Aquila earthquake sequence, causing an *M*_{w} 6.3 earthquake on 6 April. The geometrical and slip rate parameters that comprise the main input for *FiSH* are mainly taken from Pace *et al.* (2006, 2014), whereas the paleoseismological data are from Cinti *et al.* (2011). In the following, we describe in detail the results obtained using the three *FiSH* tools.

### MB Example

Following Peruzza *et al.* (2011), the last earthquake that was assigned to the fault is the 2009 L’Aquila event (*M*_{w} 6.3), which resulted in an elapsed time of 6 yr in 2015, and the uncertainty in the magnitude was fixed to 0.1. This variability is what can be assigned to the instrumental estimate of *M*_{w} by full seismic moment tensor inversion. The geometric parameters of the fault (a length along strike of 20 km, a dip angle of 50°, and a seismogenic thickness of 14 km) and slip rates (ranging from 0.6 to 0.8 mm/yr) were taken from previous geologic studies by Boncio *et al.* (2004) and Roberts *et al.* (2004), both of which were revised by Pace *et al.* (2010). The input file for MB is shown in Figure 3. The header in the input file is mandatory, and values must be given in the same order. To run MB, the user must define the path to the main folder in the MATLAB window and type “*mb*” into the command line of MATLAB to initiate an empty MB session. The MB GUI will open, allowing the user to load the input file, name the output file, and set the strain drop, the shear modulus, and eventually the number of standard deviations (*σ*) for truncating the normal distribution of magnitudes at both sides. The GUI window also contains information on the main features of the code.

MB accepts as many rows as the user needs, and the output will result in a text file that contains a row for each fault and as many figures (in .EPS format) as faults. These figures show the magnitude probability density functions, their sum, and the resulting mean and the standard deviation. In this example, the use of the relationship of Wells and Coppersmith (1994) for normal faults (scale relationship code WC94‐N), which uses the geometrical and kinematic parameters of the fault, suggests that the full rupture of the source is given *M*_{max}=6.5, which is associated with a *T*_{mean} of 922 yr. The text input and output files and output figure are shown in Figures 3 and 4, respectively. In this case, the use of 1‐*σ* truncation slightly moves *M*_{max} toward the maximum observed value (6.4), but a contemporaneous decrease in *T*_{mean} (653 yr) keeps the seismic moment rate almost fixed (from 7.678×10^{15} to 7.675×10^{15} N·m·yr^{−1}).

### RP Example

The RP tool accepts the dates of past events (historical or paleoearthquakes) as the input. An example input file, prepared using paleoseismological data from the Paganica fault (Cinti *et al.*, 2011), is shown in Figure 5. The series contains five earthquakes, of which the most recent is the A.D. 2009 earthquake. Three events occurred between A.D. 1461 and 760 B.C., and the oldest one occurred between 2900 B.C. and 760 B.C. Note that the header in the input file is mandatory, and values must be given in the same order. If no uncertainty is entered for an event, as is the case for the two most recent historical earthquakes, the youngest and oldest ages are the same. To run RP, set the path and type “*rp*” into the command line of MATLAB to start a new session. The RP GUI allows the user to load an input file and set the name of the output file and the number of simulations (in the Paganica example, we used 10,000 simulated paleoearthquake series).

The RP output figure, which is shown in Figure 6, contains the selected occurrence data (Figs. 5b and 6a), the number of simulations, the contour plots of the three two‐parameter methods (Fig. 6b,d–f), and a histogram of the average recurrence time from the Poisson distribution (Fig. 6c). With 10,000 paleoearthquake simulations, Figure 6b shows a quite large area of (*T*_{mean}, CV) pairs, with two peaks with the maximum frequency (more than 300 simulations) at values of ∼800−(0.55) and ∼1000 yrs (0.65), respectively. Figure 6c shows the histogram of *T*_{mean}, which has a range of 700–1200 yr, assuming a Poisson distribution. Figure 6d shows that the two peaks are more numerous (>600 pairs of *T*_{mean} and CV) and somewhat more clearly separated (with values of ∼800 [0.5] and 1200 yr [0.8], respectively) when using a BPT distribution to model the recurrence of the earthquake series, which suggests the two most likely recurrence behaviors of the Paganica fault. Finally, the results from the Weibull model (Fig. 6e and f) are given in terms of *a* and *b* coefficients (equation 7) and equivalent values of *T*_{mean} and CV (equations 8 and 9). The highest frequency of *a* and *b* pairs (>600 pairs) falls between 1200 and 1400 yr for *a* and between 1.3 and 1.5 for *b* and splits into two peaks whose highest frequency (>300 pairs) *T*_{mean}−CV results are between ∼800 and 1000 (0.4–0.6) and ∼1150 yr (0.7). The output also includes a text file that indicates the mode of pair(s) of parameters that are evaluated by each of the three methodsarithmetic, BPT, and Weibull and the mode of the mean recurrence time(s) that are obtained by assuming a Poisson distribution of events. These results are shown in Figure 8.

### AR Example

The structure of the input file for AR can be the same as the output file of MB (Fig. 7); conversely, the user can modify or edit a new, similar input file. To run AR, set the path and type “*ar*” into the command line of MATLAB to start a new AR session. The AR GUI prompts the user to load the input files, set the name of the output file, choose an MFD, define the probability window and the magnitude bin width, and select if the user wants the output to include a figure (one for each row of the input file) that shows hazard functions and the probability of occurrence from the BPT and Poisson models. The second input file is required only for MFDs that follow GR or truncated‐GR behavior; this second file contains the *b*‐value and *M*_{t} parameters, which were described in the Classical GR and Truncated GR sections. The GUI window also contains information on the main features of the code. The header in the input file is mandatory, and values must be given in the same order. AR accepts as many rows as the user needs and will return a file containing a row for each fault and as many figures (in .EPS format) as input rows. In the Paganica example, we use the output file from MB (Figs. 3b and 7a) to run AR with four possible MFD choices: *SingleValue*, *GaussCH*, *ClassicalGR*, and *TruncatedGR*.

Figure 8 shows the hazard functions for the time‐independent model (i.e., a Poisson distribution of earthquakes in time) and the time‐dependent model (i.e., a BPT distribution). The MFD is set to *GaussCH*, with a central value of *M*_{w} 6.5 and a standard deviation of 0.2, as listed in the input file (Fig. 7a). Because 6 yr have elapsed since the last earthquake occurred in 2009, the time‐dependent probability of the occurrence of an event with 6.3≤*M*_{w}≤6.7 over the next 50 yr is ∼0.00001% (*P*|elap red circle in Fig. 8). The Poisson probability of an earthquake with 6.3≤*M*_{w}≤6.7 in 50 yr is about 5% (Figs. 7c and 8). Together with the probability of occurrence, AR returns the annual earthquake rates for the given magnitude bins in the form of incremental rates, as shown in Figure 7b.

### Comparing Rates

In Figure 9, we compare the cumulative rates of occurrence that were computed by the four chosen MFDs with the rates that were calculated by RP, which only used paleoseismological data (Fig. 5). Because RP does not include the magnitude of the events, the rate of occurrence (i.e., the inverse of the mean recurrence time, which is computed using the arithmetic method in Figure 5c) is plotted as a horizontal line in Figure 9. The curves that were obtained from AR share the same seismic moment rate from the input file, 7.678×10^{15} N·m·yr^{−1} (Fig. 3b). Interestingly, the moment rate in an ideal situation, in which the fault rupture is dictated by single‐value behavior and the paloeseismological data are complete, should be the same as that from MB. In contrast, assuming complete paleoseismological data and that the fault ruptures in the greatest earthquakes have magnitudes that follow a Gaussian distribution, the horizontal line and the characteristic Gaussian (*GaussCH*) model should intersect at the lower bound of *M*_{max} (6.3 in this case), as should the GR and truncated‐GR cases. In Figure 9, the RP‐derived rate intersects three of the AR‐derived rates (the characteristic Gaussian, GR, and truncated GR) and lies below the single‐value model. We suggest two possible interpretations according to the limitations of the input data: (1) the paleoseismological data are complete for events with magnitudes in the range 6.3–6.7, if the fault has characteristic Gaussian behavior, and (2) the paleoseismological data are not complete for events with magnitudes lower than 6.1–6.2, if the fault has GR or truncated GR behavior, or the paleoseismological data are complete, and the slip rates overestimate the seismic moment that is released by the fault.

From a seismic‐hazard perspective, paleoseismologically inferred rates of occurrence should be cross‐checked with geologically derived long‐term slip rates and/or with rates of occurrence that are derived from other approaches (such as geodetic rates that are computed for seismogenic faults). Our intent is to distribute tools that can help researchers to pinpoint potential inconsistencies and obtain more reliable fault‐based seismic‐hazard evaluations.

## FINAL REMARKS

*FiSH* is available at http://fish-code.com (last accessed January 2016) and should be run in MATLAB R2011b. The compressed file also contains a collection of sample data files (inputs and outputs). Although tested and designed to work on any MATLAB‐supported platform, it is possible that some users will experience problems while running various functions. The source codes are open, and we encourage users to handle the scripts, communicate with us regarding bugs, and/or suggest further improvements. *FiSH* is intended to be useful for seismic‐hazard practitioners, but it must be used under the general principle “trash in–trash out.” The use of incomplete earthquake series and/or earthquake associations in RP or unrealistic geometries in MB can be only partially detected by the code. In our experience, the output figures help the user to detect input data problems. Although they have not been scheduled yet, new releases of *FiSH* are likely. These updates will largely be driven by research interests, and new features will be documented on the same website.

## DATA AND RESOURCES

*FiSH* is written in MATLAB (v.R2011b), a Mathworks software (http://www.mathworks.com, last accessed December 2015) and is available at http://fish-code.com (last accessed January 2016). All data used in the examples for moment budget (MB Example), recurrence parameter (RP Example), and activity rate (AR Example) came from published sources listed in the references.

## ACKNOWLEDGMENTS

This study has benefitted from partial funding that was provided by the Italian Presidenza del Consiglio dei Ministri Dipartimento della Protezione Civile (INGV‐DPC), Project S2‐2012 and V3‐2014, and the Ministero dell’Istruzione, dell’Università e della Ricerca (MIUR) funded project "Futuro in Ricerca" (FIRB) Abruzzo (Code RBAP10ZC8K_006). This article does not necessarily represent DPC’s official opinion and policies. We thank the reviewers for their accurate reviews and useful suggestions.