Visual analytics of aftershock point cloud data in complex fault systems

Aftershock point cloud data provide direct evidence for the characteristics of underground faults. However, there has been a dearth of studies using state-of-the-art visual analytics methods to explore the data. In this paper, we present a novel interactive visual analysis approach for visualizing the aftershock point cloud. Our method employs a variety of interactive operations, rapid visual computing functions, flexible display modes, and various filtering approaches to present and explore the desired information for the fault geometry and aftershock dynamics. The case study conducted for the 2016 Central Italy earthquake sequence shows that the proposed approach can facilitate the discovery of the geometry of the four main fault segments and three secondary fault segments. It can also clearly reveal the spatiotemporal evolution of the aftershocks, helping to find the fluid-driven mechanism of this sequence. An open-source prototype system based on the approach is also developed and is freely available.

1 Introduction

Fault complexity of earthquakes is universal in all types of faults (strike-slip, normal, and reverse). Usually, the faults are divided into a number of subparallel segments along the strike by geometrical discontinuities (Manighetti et al., 2015). In some cases, the depth segmentation of a fault is also observed (Elliott et al., 2011) where the cascading earthquake occurs successively along the dip. The geometrical discontinuities in the complex fault structure can sometimes impede and stop the rupture propagation and limit the earthquake magnitude. In some cases these act as a high-permeability conduit to allow fluid migration and control the spatiotemporal sequence of the subsequent earthquakes (Walters et al., 2018). The elastic-rebound theory based on the earthquake cycle concept is commonly applied for long-term earthquake prediction. However, the complex fault structure complicates the seismic hazard assessment because theoretical models based on simplistic assumptions (e.g., single dynamic rupture) were found to be inapplicable to natural faults (Barbot et al., 2012). The time intervals for cascading fault rupture due to fault segmentation range from seconds to years on a case-by-case basis, while the control factors of the faults are still not clearly understood. Therefore, it is highly important to understand the complex fault geometry, and such understanding may provide useful insight into fault mechanics and help to constrain the existing theoretical models.

Due to the past lack of ability to conduct high-resolution observations, it is always difficult to define the geometry of the active fault system. Recent developments of both geodetic and seismic observation techniques make it possible to reveal the fault geometry at the surface and depth. Advanced synthetic aperture radar (SAR) geodesy provides near-field deformation and can be used to directly map the fault surface and constrain the deep fault through deformation inversion. The construction of an increasing number of seismic stations allows for the identification of the locations of a large number of aftershocks following the main shocks. The point cloud of the located aftershocks contains the information regarding both the three-dimensional coordinates and precise time of occurrence that can directly reveal the fault geometry and temporal evolution of the earthquake sequence. Current interpretation methods for the aftershock cloud are mostly case-specific rather than general, and to date, no standard processing procedure for the visualization of these data has been available. While many studies plot the cross-section profile of the aftershock point cloud to reveal the dip of the fault plane (Irmak et al., 2012; Hengesh and Whitney, 2016; Li et al., 2010; Elliott et al., 2013), these simple visualization methods usually have the following limitations:

1. Lack of interactive operations: the visualization methods provide a static display without an interactive interface to assist the users in effectively performing analyses through various human–computer interactions.

2. Lack of visual analysis: the current methods simply present these data without computing the feature (e.g., plane fitting or local outlier factor value calculation) to enhance the structure.

3. Lack of noise filtering: while the aftershock locations normally contain many uncertainties, most studies do not consider either the identification or removal of these noisy points.

As was shown in many recent studies, aftershock data can reveal much more information if these data are presented well, and this is particularly true for a complex fault system (Walters et al., 2018). Therefore, it is necessary to develop an interactive visual system to effectively mine information from the aftershock data.

Visual analytics is a research direction of the field of data visualization (Wong and Thomas, 2004), which integrates new computational and theory-based tools with innovative interactive techniques and visual representations to enable human–information discourse. How to efficiently transform large earth observation data into knowledge is a tremendous challenge in the era of big data (Xia et al., 2018). Visual analytics is highly advantageous and has good potential in complex data analytics for intuitive information representation and dynamic data exploration (Ward et al., 2015). The importance of scientific visualization has been recognized for a long time (McCormick, 1988). Specifically, in geophysics, the visualization of seismic refraction data has already been widely applied to explore the distribution of petroleum and gas (Liu et al., 2018; Raposo et al., 2016). However, to date, studies using state-of-the-art visual analytics methods to explore other kinds of seismic data (e.g., aftershock point cloud) to gain information about earthquakes have been lacking.

In this study, we propose a novel interactive visualization method for exploring the 3-D aftershock point cloud data. This method can help researchers to better understand the complex fault geometry in three ways:

1. A set of interactive operations is provided to stimulate creative analysts and exploit the researchers' background knowledge.

2. On-the-fly computation of the local outlier factor (LOF), 2-D projection, and plane fitting are supported in the visual analysis so that the hidden features of the fault structure can be discovered dynamically.

3. Noisy or irrelevant aftershock points can be filtered to reduce the distractions from invalid information.

Figure 1 Flowchart of the proposed approach.

2 Visual analytics of aftershock data

Figure 1 shows the flowchart of the proposed approach. Initially, we imported the aftershock point data from an earthquake catalog or the literature into our system. At the same time, the LOF values for all of the points are computed automatically. Then, the points are presented in a 3-D view window, allowing free rotation, panning, and zooming to observe the point structure. The earthquake fault, with data either from an existing source or from inferred information, can be imported into the 3-D view as well. After defining the parameters for the projection (i.e., strike and dip of the plane and azimuth of the projection direction), the 2-D points are rapidly computed and presented. Then, some interactions, i.e., coloring and filtering, can be performed on the 3-D point clouds and 2-D plane projection to enhance the structure. From the 2-D plane projection, two interactions are further provided for exploring the fault geometry and aftershock dynamics. The first is the selection and plane fitting of a point cluster on a potential plane. After interactively selecting the clustered points, a plane will be automatically fitted by a robust estimator to reveal the potential fault structure. The second is the observation of aftershock migration. Animation of the aftershocks is provided so that a preliminary recognition of the migration direction can be obtained. Then, the user can define a propagation direction to further explore the aftershock evolution. After assigning the direction, the distance traveled by the aftershocks along the specified direction as a function of time will be plotted, which is helpful in discovering the driving factor of aftershock migration. We note that when the preliminary exploration of the fault geometry or earthquake dynamics is finished, the users can return to the previous interactions to refine their interpretation of the data.

2.1 Visual design

As shown in Fig. 2, the user interface includes four main components. The first component is the operation tools panel with all the buttons, sliders, text boxes, and labels. It includes most interactions in the visual analytics. The second component is the 3-D view panel for presenting cloud points and fault planes. The third component is the 2-D view panel for presenting the points in the 2-D projection plane. The fourth component is a temporal panel that switches between 3-D fault plane fitting plotting and propagation distance–time plotting, corresponding to the explorations of the fault geometry and aftershock dynamics, respectively.

Figure 2 Visual interface of our approach consisting of four interrelated components: (a) operation tools, (b) a 3-D view, (c) a 2-D projection view, and (d) a temporal view for plotting a fitted fault plane or distance–time function.

2.2 Interactive operations

The easy-to-use interactions represent a significant difference between the proposed method and conventional aftershock interpretation methods. The traditional static displays prevent the analyst from directly and rapidly exploring the aftershock data. However, the easy-to-use interactive displays can lead to deeper insights into 3-D aftershock cloud point data and the discovery of new fault structures or spatiotemporal evolution patterns. These interactive operations are described as follows.

Three-dimensional view manipulation. The point cloud data can be displayed in a 3-D view control (Fig. 2b). Some typical 3-D interactions are available including zooming, panning, rotating, data tips, and data brushing (Fig. 2a), allowing interactive exploration and editing of the plotted data to improve the visual display.

Interactive projection. The cross-section plot is the main visualization approach for exploring the complex fault geometry (Fig. 2c). The faults are normally represented by linear features when the aftershocks are projected onto a plane perpendicular to the fault. We offer an interface to input the strike and dip of the projection plane for rapid plotting of a desired plane projection profile.

Dynamic filtering. Noisy or irrelevant data can significantly affect the interpretation of the aftershock point cloud. The filtering or denoising interaction is an important step in removing noisy or irrelevant aftershock points while retaining the fault geometry details. The filtering is based on the four features of the aftershocks, namely the depth, LOF, time, and magnitude (Fig. 2a). LOF is an index for finding anomalous data points that will be introduced below (Breunig et al., 2000). Filtering based on LOF is used to denoise the data.

Color switching. The feature for coloring the aftershock points can be quickly switched between the depth, LOF, time, and magnitude, allowing the users to effectively understand the data from different dimensions (Fig. 2a).

Fault plane fitting. To test the probability of an identified fault, the users first brush the points that form a linear feature in the cross-section plot (Fig. 2c). A fault plane will be automatically fitted to these points based on a robust estimator (Fig. 2d). This interaction combines expert knowledge and visual computing, allowing the user to quantitatively locate the potential fault structure.

Fault data fusing. Users can import the faults to the 3-D view in order to assist in the visualization and analysis of the aftershock point cloud (Fig. 2b). These faults can be hypothetical seismogenic faults or background fault systems.

Animation. The aftershocks can be plotted gradually to form an animation in both the 3-D view and the 2-D cross-section plot (Fig. 2b and c). An animated visualization has unique advantages and aids the viewers in interpreting the evolution of the aftershocks, making comparisons between different points in time, chunking the aftershock points into smaller datasets, focusing attention, and anticipating and understanding the spatiotemporal changes in the aftershocks.

Propagation distance–time plotting. To further quantify the migration process of the aftershocks, users can define a propagation direction (Fig. 2c) and plot the distance along this direction as a function of time (Fig. 2d). The temporal trend may reveal the mechanism (e.g., fluid diffusion) underlying the aftershock migration.

2.3 LOF calculation

Because of the limitations due to the resolution of seismic data inversion, there is some uncertainty in the location of the aftershocks, so that some earthquake points show large deviations from the fault plane. Usually, the original aftershock data do not contain any information about this uncertainty, and we cannot eliminate these outliers. In this study, we proposed to use the LOF, which is commonly used in the identification of the outliers in the lidar point cloud (Wang et al., 2019) to detect the anomalous points that deviate distinctly from the fault plane. The LOF (Breunig et al., 2000) is based on the concept of local density, which is estimated by the reachability distance of the object to its k nearest neighbors. Points with a substantially lower density than their neighbors are considered outliers. By comparing the local reachability density of the object to those of its neighbors, its LOF value can be obtained. The LOF is calculated as follows.

Let kd(p) be the k distance of a point p, defined as the distance from point p to the kth nearest neighbor. The reachability distance between points p and o is then given by

$$\begin{array}{}\text{(1)}& r{d}_k\left(\mathrm{p},\mathrm{o}\right)=max\left\{kd\left(\mathrm{o}\right),d\left(\mathrm{p},\mathrm{o}\right)\right\},\end{array}$$

where d(p,o) is the actual distance between p and o.

The local reachability density of an object p is then defined by

$$\begin{array}{}\text{(2)}& lr{d}_k\left(\mathrm{p}\right)={\displaystyle \frac{\left|N_k\left(\mathrm{p}\right)\right|}{{\sum }_{\mathrm{o}\in N_k\left(\mathrm{p}\right)}r{d}_k\left(\mathrm{p},\mathrm{o}\right)}},\end{array}$$

where N_k(p) is the kth nearest neighborhood, defined as the set of points within the distance of the k distance (p).

The LOF for point p is then obtained by

$$\begin{array}{}\text{(3)}& {\mathrm{LOF}}_k\left(\mathrm{p}\right)={\displaystyle \frac{{\sum }_{\mathrm{o}\in N_k\left(\mathrm{p}\right)}lr{d}_k\left(\mathrm{o}\right)}{\left|N_k\left(\mathrm{p}\right)\right|lr{d}_k\left(\mathrm{p}\right)}}.\end{array}$$

Figure 3 shows an example of the calculated LOF values, where an LOF value significantly larger than 1 indicates a sparse region.

Figure 3 An example of LOF scores for 2-D points. Each point is colored by its LOF score. The radius of red circles on five selected points is in direct proportion to their LOF score (plotted in the center). It is observed that the central points with high density have low LOF scores, while the surrounding points have high ones.

2.4 Two-dimensional projection

Projecting the 3-D aftershock point cloud to a 2-D plane can significantly enhance the fault geometry features. As mentioned above, we offer an interactive way to rapidly compute the projection when a plane and projection direction are assigned. To make the interactive operation user-friendly for geologists, we use geological terms, i.e., strike and dip (to define the orientation of the projection plane), and azimuth angle (to define the projection angle). We firstly calculate 3-D geographic coordinates of the projected points on the projection plane for aftershocks according to the line–plane intersection equation. The 3-D geographic coordinates are then transformed to a 2-D local coordinate system by two element rotations: a z axis rotation by an angle of strike − 90^∘ and an x axis rotation by an angle of − dip, as shown in Fig. 4a.

Figure 4 (a) Illumination of a projection to a plane defined by strike and dip. (b) Strike and dip parameter calculation from a normal vector of a fitted plane.

The rotation matrices (R_z and R_x) are given by

$$\begin{array}{}\text{(4)}& {\mathbf{R}}_z\left(\text{strike}-\mathrm{90}{}^{\circ }\right)=\left[\begin{array}{ccc}\cos (\text{strike}-\mathrm{90}{}^{\circ })& -\sin (\text{strike}-\mathrm{90}{}^{\circ })& \mathrm{0}\\ \sin (\text{strike}-\mathrm{90}{}^{\circ })& \cos (\text{strike}-\mathrm{90}{}^{\circ })& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right]\end{array}$$

and

$$\begin{array}{}\text{(5)}& {\mathbf{R}}_x\left(\mathrm{90}-\text{dip}\right)=\left[\begin{array}{ccc}\mathrm{1}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \cos \left(-\text{dip}\right)& -\sin \left(-\text{dip}\right)\\ \mathrm{0}& \sin \left(-\text{dip}\right)& \cos \left(-\text{dip}\right)\end{array}\right].\end{array}$$

The projected coordinates (P) from the aftershock point (u) is given by

$$\begin{array}{}\text{(6)}& \mathbf{P}={\mathbf{R}}_x{\mathbf{R}}_z\mathit{u},\end{array}$$

where the x and y coordinates of P correspond to the coordinates in the projection plane, and the z coordinate represents the distance to the plane.

2.5 Plane fitting

When the analyst identifies a potential linear feature from the cross-section plot, the plane fitting can rapidly evaluate the possible feature and provides the precise geometric parameters. In this study, we used an algorithm based on the singular value decomposition (SVD) to fit a plane to the given 3-D aftershock points. The calculation of the parameters of the plane ($ax+by+cz=d$) based on the singular value decomposition method is carried out as follows.

Given a set of points $\left(x_i,y_i,z_i,i=\mathrm{1}\mathrm{\dots }n\right)$ and their center point ($\overline{x},\overline{y},\overline{z}$), we form a matrix A as given by

$$\begin{array}{}\text{(7)}& \mathbf{A}=\left[\begin{array}{ccc}{x}_{\mathrm{1}}-\overline{x}& y_{\mathrm{1}}-\overline{y}& z_{\mathrm{1}}-\overline{z}\\ x_{\mathrm{2}}-\overline{x}& y_{\mathrm{2}}-\overline{y}& z_{\mathrm{2}}-\overline{z}\\ x_{\mathrm{3}}-\overline{x}& y_{\mathrm{3}}-\overline{y}& z_{\mathrm{3}}-\overline{z}\\ .& .& .\\ .& .& .\\ x_n-\overline{x}& y_n-\overline{y}& z_n-\overline{z}\end{array}\right].\end{array}$$

We apply an SVD to matrix A:

$$\begin{array}{}\text{(8)}& \mathbf{A}=UD{V}^{\mathrm{T}}.\end{array}$$

The parameter vector $[a,b,c]$ is just the smallest singular vector corresponding to the least singular value, and the last parameter d can be obtained by

$$\begin{array}{}\text{(9)}& d=a\overline{x}+b\overline{y}+c\overline{z}.\end{array}$$

To avoid disturbance from noisy points, we include an iterative step to identify outliers and exclude them from plane fitting. This step is carried out using the following algorithm:

1. Fit a plane using the available points.

2. Calculate the distance (D) of the points to the plane.

3. Calculate the root mean square error (RMSE) of all of the distances (D).

4. Find points with distances 3 times larger than RMSE and define them as outliers.

5. Remove these points, and then return to step 1.

6. Stop when there is no outlier found, and return the plane parameter.

Given the normal vector of the fitted plane ($a,b,c$), the strike and dip of the fault plane are given by (Fig. 4b)

$$\begin{array}{}\text{(10)}& \text{strike}=\text{atan}\mathrm{2}\left(\mathrm{sign}\left(c\right)\times \mathrm{a},\text{sign}\left(c\right)\times b\right)\end{array}$$

and

$$\begin{array}{}\text{(11)}& \text{dip}=\text{acos}(c/\sqrt{a^{\mathrm{2}}+b^{\mathrm{2}}+c^{\mathrm{2}}}),\end{array}$$

where atan2(Y,X) is the four quadrant arctangent of the elements X and Y, $-\mathit{\pi }<=\mathrm{atan}\mathrm{2}(Y,X)<=\mathit{\pi }$, and sign(.) is the signum function that returns 1 for a positive element and −1 for a negative element.

2.6 Prototype system

The proposed visualization approach has been implemented in a software package (Aftershock Visualization, AFV1.0; https://github.com/caigenszu, last access: 22 August 2019) and is freely accessible to the scientific community. This prototype system was developed in MATLAB with a graphical user interface (GUI), providing point-and-click control for all operations. An analyst with no programming background does not need to learn a new language or type commands to run the application. Meanwhile, the source code of the package is also provided so that advanced users can easily read the data and results throughout the visualization processing, and they can modify the package to include their custom-built modules. This package requires an input file consisting of the aftershock point cloud table in plain text with each row representing a point and the columns representing UTM-x, UTM-y, depth, magnitude, time (days from reference time), LOF, and an event flag (whether or not the shock is the main shock) in that order. An alternative input file is the fault plane used to assist data interpretation. Each fault plane is represented by one row, and the columns represent strike, dip, center UTM-x, center UTM-y, length, upper depth, and bottom depth.

3 Case study: the 2016 Central Italy earthquake sequence

On 24 August 2016 (01:36 UTC, local time 03:36), a destructive earthquake (M_w 6.2) occurred in Central Italy (the Amatrice earthquake). The epicenter was located at 42.70^∘ N, 13.23^∘ E, between the towns of Norcia and Amatrice. Two months later, two major earthquakes were triggered to the northwest of the Amatrice earthquake on 26 October (19:18 UTC) and 30 October (06:40 UTC) (Papadopoulos et al., 2017), with individual moment magnitudes of 6.1 (Visso earthquake) and 6.6 (Norcia earthquake). The earthquake sequence occurred across several fault systems (Fig. 5), including the Mt. Vettore–Mt. Bove Fault system (VBFS) and the Mt. Gorzano Fault system (GFS). The VBFS is a SW-dipping fault system composed of different splays and segments. The fault system scarps are exposed along the SW foothills of Mt. Vettore, Mt. Porche, and Mt. Bove, with a length of ∼27 km. The GFS is a SW-dipping fault with a ∼26 km long fault scarp developing along the foothills of Mt. Gorzano. These two faults are segmented by existing tectonic structures inherited from pre-Quaternary compressional tectonics (Pizzi and Galadini, 2009). A ∼3 km thick layer in which small events and some large extensional aftershocks occur is found below the seismogenic fault, limiting the seismicity to the first 8 km of the crust (Chiaraluce et al., 2017).

Figure 5 Seismotectonic setting of the Central Italy earthquake sequence. The three major earthquakes in the sequence are denoted by the red stars and beach ball symbols, and the two recent large historical events are colored in black. The dashed black rectangles represent the surface projection of the fault planes adopted in this study. The bold black lines are the two seismogenic normal fault systems, namely, the Mt. Vettore–Mt. Bove Fault system (VBFS) and the Mt. Gorzano Fault system (GFS). The pink line shows the simplified trace of the preexisting compressional front named the Olevano–Antrodoco–Sibillini (OAS) Thrust. Aftershocks are marked by the small dots with different colors. The blue, purple, and yellow points represent the events taking place on 24 August–24 October 2016, 24–30 October 2016, and 30 October 2016–8 October 2017, respectively.

Many recent studies reported in the literature have revealed a very complex fault geometry for this earthquake sequence. In addition to the main range-front normal faulting structures that were adopted in many early studies (Liu et al., 2017; Tinti et al., 2016), recent studies find increasing evidence for some additional minor antithetic and synthetic faults (Chiaraluce et al., 2017; Scognamiglio et al., 2018; Walters et al., 2018; Cheloni et al., 2017). The activation of this secondary fault structure has important implications for evaluating the seismic hazard in this section of the Central Apennines. The complex fault geometry also plays an important role in channeling fluid diffusion and affects the spatial order of the cascading earthquakes. The aftershock distribution provides direct evidence for the fault geometry and underground dynamics. By interactively visualizing these data, we can obtain useful information regarding both the fault geometry and earthquake dynamics. Additionally, a relocated aftershock dataset has been provided in the work of Chiaraluce et al. (2017), making this earthquake sequence an ideal case for applying our proposed interactive visualization method. Below, we demonstrate the benefits of our method in identifying fault geometry and aftershock migration.

3.1 Fault geometry identification

According to previous studies, several fault segments are involved in this earthquake sequence. Some of these are the main fault segments located in the GFS–VBFS, while others are secondary fault segments including a NE-dipping normal fault antithetic to the GFS–VBFS and a preexisting compressional structure that is likely related to a segment of the Olevano–Antrodoco–Sibillini Thrust. We will show how our interactive visualization method can facilitate the recognition of these structures.

Following the flowchart shown in Fig. 1, we first import the relocated aftershock data from the results obtained by Chiaraluce et al. (2017) to the 3-D view. Some of the basic interactions provided by this package, i.e., zooming, panning, rotating, data tips, and data brushing, can be used to explore the data. LOF values were then calculated for all points to identify the aftershock outliers. We can interactively select to show some of the points according to their location, depth, magnitude, LOF, and time, and we can project them to a desired plane. After we find a linear feature in the projection plane, we select the points comprising the linear feature and automatically fit a plane to the feature (Fig. 6).

Figure 6 Fault geometry identification and plane fitting. (a) A three-dimensional view of seven identified fault segments and an aftershock point cloud. Thick lines represent the surface traces of the faults, and the dashed lines are the underground fault boundaries. (b) Projections of the 3-D aftershocks. Red dashed lines denote the seven fitted fault planes.

First, we identify the four main fault planes by projecting the cloud point to a plane with a 70^∘ strike and a 90^∘ dip. This plane is approximately perpendicular to the GFS–VBFS, and it should reveal the cross sections of the fault plane well. To better illuminate the linear features of the faults, we filter the aftershock points according to the fault location and the spanning time of the aftershocks on the fault. To denoise the point cloud, the aftershocks with the LOF value larger than 2 are filtered. We then select the points around the linear feature for each fault plane, and we automatically obtain the suggested geometry parameters of the fitted plane as shown in Table 1, finding that the obtained parameters are in overall agreement with the main seismogenic fault planes in previous studies.

Table 1 Geometry parameters of the inferred fault from the aftershock point cloud. Strike and dip are given in degrees; all other parameters are measured in kilometers.

Second, we explore the point cloud to find other secondary faults. Since there is no hint of the existence of these faults, we use the basic interaction tools in the 3-D view and the filtering function in our package to find the possible fault geometry. An obvious cluster of aftershock points is found on an antithetic NE-dipping plane in the Norcia area. We therefore also project the cloud point to a plane with a 70^∘ strike and a 90^∘ dip. After filtering the points by location and LOF value, we obtain a very strong linear fault signal as shown in Fig. 6. Fault fitting gives a strike of −38^∘ and a dip of 59.8^∘ for the Norcia antithetic fault. However, other fault structures revealed in the literature are not very obvious in the aftershock cloud point. We then project the point cloud to a simplified seismogenic fault ($\text{strike}=\mathrm{160}{}^{\circ }$, $\text{dip}=\mathrm{50}{}^{\circ }$) along a horizontal vector with an azimuth of 15^∘, using the projection settings from Walters et al. (2018). We filter the points by the magnitude, location, timing, and LOF, and then, two linear structures are observed as shown in Fig. 6. The structure on the right is the NW-dipping Pian Piccole fault that generates a significant surface displacement signal in the interferometric synthetic aperture radar (InSAR) measurements. The structure on the left is an east-dipping fault that was only inferred in Walters et al. (2018). We pick the points around these linear structures for plane fitting and obtain a strike of −114^∘ and a dip of 46.5^∘ for the Pian Piccole fault and a strike of 7.8^∘ and a dip of 85.3^∘ for the inferred fault.

3.2 Migration of aftershocks

Besides the static analysis of aftershocks to find the fault geometry, the timing information of aftershocks has the potential to reveal earthquake dynamics related to fluid diffusion, afterslip, and/or pressure transients. Fluid diffusion is a main factor in triggering seismicity as observed in many tectonic regimes (Vidale and Shearer, 2006). It was also proposed to have been present in the nearby 2009 L'Aquila earthquake (Malagnini et al., 2012). The migration of the aftershocks could provide some evidence for such diffusive processes. Walters et al. (2018) found a significant fluid-driven aftershock migration between the Amatrice earthquake and Visso earthquake controlled by the fault intersection in this earthquake sequence.

Figure 7 Observation of aftershock migration. (a) Aftershock on the projection plane colored by time. The arrow represents north. It is observed that most aftershock points in the north happened more recently than in the south, implying a northward migration trend. (b) Propagation distance as a function of time. The red line represents the temporal trend that the propagation distances increase with time.

To observe the spatiotemporal process, we project the aftershocks of the seismogenic fault to be along a horizontal vector with an azimuth 15^∘ as described above. The points are colored by their rupture time. As shown in Fig. 7a, we can observe a trend by the color change from green to red in the northward direction, implying the clear directionality of the aftershock migration. We then use the animation operation to plot the aftershock points step by step as a function of time in both the 2-D projection plane and the 3-D view. The directional propagation of the aftershocks can be observed much more clearly.

Next, we can define a starting point and an assumed direction in the projection plane. After drawing the starting point and the propagation line, the distances of the aftershocks to the starting point along this direction can be calculated automatically. They are then plotted as a function of time to further reveal the time evolution of the aftershocks. Following the work of Walters et al. (2018), we defined a starting point in the peak slip of the Amatrice earthquake and a direction perpendicular to the Pian Piccolo fault. A clear temporal trend of the aftershocks was then observed in the plot (Fig. 7b). Therefore, the observed spatiotemporal pattern of aftershock migration can be interpreted by the seismologist to understand the earthquake mechanism and infer the possible controlling factor for the fault rupture. In this case, the diffusive-like aftershock temporal trend can be further interpreted using a pore fluid source model. Through the interactive visualization operations, we can observe that aftershock propagation is consistent with a diffusive process and can also observe the spatial concentration of the aftershocks along the fault intersection. This information can be used to derive the conclusion, as described in the work of Walters et al. (2018), that the intersecting structures act to channel the fluids and control the timing and order of the subsequent earthquakes.

4 Conclusion

We present a novel interactive approach and develop a prototype system to illustrate 3-D aftershock point clouds that can help the geophysicist and seismologist to better understand the geometry of a complex fault system and the spatiotemporal pattern of the aftershocks. Moreover, it can be applied in an educational scenario to facilitate student learning. In this approach, we design a set of interactive operations to facilitate the exploration of the data. Fast computation of LOF, 2-D projection, and plane fitting are implemented in addition to the visualization to reveal additional, hidden information, which is not obvious in the original data. Additionally, the point cloud can be visualized in many ways (e.g., 3-D view, 2-D plane projection, and animation), aiding viewers in interpreting the aftershock point cloud. A wide range of data filtering options (e.g., by magnitude, depth, location, time, and LOF) enable the users to avoid interference from invalid data in order to focus on the useful information.

It is widely accepted that visual analytics with interactive visual interfaces can amplify human cognitive capabilities, and it is suitable for solving some complex problems relying on closely coupled human and machine analysis. When applied to the aftershock catalog, visual analytics can therefore help to discover some new fault planes or weak signals of aftershock migrations, which might be hard to observe directly from the data. Specifically, the designed visual analytics procedure can assist the knowledge discovery from the aftershock catalog in several ways:

1. by accelerating the fault segment discovery processing through 3-D view manipulation (e.g., zooming, panning, and rotating);

2. by enhancing the recognition of fault segment patterns through rapid visual computing functions (e.g., interactive projection, plane fitting, and fault data fusing);

3. by reducing the influence from low-quality or irrelevant aftershock points through various filtering methods (e.g., depth, LOF, time, and magnitude filtering);

4. by enabling the aftershock migration exploration by providing more cognitive resources (e.g., animation and propagation distance–time plotting).

A case study of the 2016 Central Italy earthquake sequence shows that the proposed approach can complete a series of complex visualization operations in a user-friendly interactive manner. We can conveniently infer the fault geometry for both the main faults and the secondary faults, and we clearly observe the spatiotemporal pattern of the fluid-driven aftershock migration. In future work, we would like to examine the capability of this tool for more aftershock point data applications, such as analyzing tremor migration and tracking aseismic slip. Also, we plan to include more kinds of fault information (e.g., InSAR, GPS, source model, and fault scarp) in the visualization in order to further facilitate the analysis of earthquake fault behavior.