A self-organized critical model and multifractal analysis for earthquakes in Central Alborz, Iran.

This paper is devoted to a phenomenological study of the earthquakes in central Alborz, Iran. Using three observational quantities, namely the weight function, the quality factor, and the velocity model in this region, we develop a modified dissipative sandpile model which captures the main features of the system, especially the average activity field over the region of study. The model is based on external stimuli, the location of which is chosen (I) randomly, (II) on the faults, (III) on the low active points, (IV) on the moderately active points, and (V) on the highly active points in the region. We uncover some universal behaviors depending slightly on the method of external stimuli. A multi-fractal detrended fluctuation analysis is exploited to extract the spectrum of the Hurst exponent of the time series obtained by each of these schemes. Although the average Hurst exponent depends slightly on the method of stimuli, we numerically show that in all cases it is lower than 0.5, reflecting the anti-correlated nature of the system. The lowest average Hurst exponent is found to be associated with the case (V), in such a way that the more active the stimulated sites are, the lower the average Hurst exponent is obtained, i.e. the large earthquakes are more anticorrelated. Moreover, we find that the activity field achieved in this study provide information about the depth and topography of the basement, and also the area that can potentially be the location of the future large events. We successfully determine a high activity zone on the Mosha Fault, where the mainshock occurred on May 7th, 2020 (M[Formula: see text] 4.9).

Introduction

It is widely believed that the earthquake is a self-organized critical system. When an earthquake occurs and a fault slips, it causes an excess in the tension of the neighboring regions. The motion of the neighboring regions depends on their local strain, so that if it exceeds a threshold, then it slips. This dynamics is reminiscent of the sandpile dynamics, first invented by Bak-Tang-Wiesenfeld (BTW) model, where the energy spreads the system based on similar dynamic rules: when the local energy exceeds a threshold, the site topples, rising the energy of the neighboring sites by one unit. BTW and some other variants[1,2] are too unrealistic to give an acceptable description of the earthquake, e.g. it cannot explain noise that is seen in reality[3]. For more realistic situations, one needs more detailed model paying attention to the structure of the earth and also the information on how the seismic activities affect each other, which translates to how a perturbation propagates from region to region. The latter is crucial and pretty complicated since it depends on the material content of the Earth’s interior where the signal propagates. Recently, using a virtual seismometers, the correlations between seismic activities were incorporated in the model, based on which a complex network was designed on top of which sandpile dynamics were implemented[4]. The virtual seismometer can provide inter-event empirical Green’s function in the Earth’s interior. Therefore, using unconventional form of seismic interferometry, one earthquake beneath the Earth’s surface can turn to a receiver where recorded another event waveform[5]. One may think of this problem from another point of view: the cross-correlation between previously happened earthquakes gives us a set of valuable information about the structure of the earth and the signal propagation in the region. For example, suppose that we have a map (the place and the magnitude, M) of previous earthquakes in the region of interest, along with the corresponding time series. Cross-correlation between the events gives us a criteria of how events are related and to which extend they are correlated.

As the case study, we focus on the earthquakes in central Alborz in the present paper. Alborz range with seismic active east-west trending mountain belt extends across the north of Iran. The south Caspian block to the north, central Iran micro-plateau to the south surround the central Alborz with several folds and various faults (see Fig. 2). Three major tectonic events can invoke for Central Alborz including: I) shortening (led to thrusting and folding[6]), II) extension (led to Damavand Volcanism[7,8]), and III) collision-related compression (from middle Miocene to recent[9]). Studies of the crustal deformation by GPS measurements  [10] indicated compression between the Central Iranian micro-plateau and South Caspian blocks evince a partitioning  5 mm yr−1 range-perpendicular and  4 mm yr−1 along range-parallel. This region is very active and experiences large, and catastrophic earthquakes (e.g., the catastrophic Manjil-Rudbar in July th, 1990 with M 7.4), which are associated with major active faults (red line in Fig. 2). Historical[11-14] and instrumental recorded earthquakes in the central Alborz represent that many faults have the potential of an earthquake up to M 7.5.

Figure 2

The map of the study area. The know faults, circles and volcano depict by solid red lines, earthquakes and brown triangle, respectively. The pink square represent Tehran, capital of Iran. The historical events reported by[11–14] and instrumental earthquakes in the central Aborz located by Iranian Seismological Center, IrSC (http://irsc.ut.ac.ir). The white box shows a region with a fair resolution (Figs. 2–6). The black rectangular in the inset map indicates the study area. This figure has depicted by Generic Mapping Tools (GMT) version 6.2.0 ([31], https://www.generic-mapping-tools.org/).

The character of the seismic wave propagation effects directly depends on the nature of the distribution of the elastic parameters within the Earth’s interior. Besides the seismic source functions (focal mechanisms, the rupturing algorithms, and time duration, etc.), the elastic parameters of the seismic wave propagation (e.g., seismic velocity structures, attenuation models, etc.) are calculated using these recorded waveforms with a combination of the classical and new seismological processing methods. Today, different methods exist for calculating 1D, 2D, or 3D structures of Earth and various properties of seismic records are used, including amplitudes, travel times, full waveforms, etc. Inversion of the arrival times of seismic (body or surface) waves is one of the routines and popular techniques for imaging Earth’s interior[15,16]. Several studies about the crustal velocity structure of the Alborz range have been done using the traveltime tomography method and different kinds of seismic recording components. A 1D velocity model with two sedimentary layers ( km s−1; with a total thickness of 8 km) overlying two crystalline layers ( km s−1) was developed by Ashtari et al.[17] which employed the first arrival P- and S-waves inversion. This 1D model was updated by [18] with a similar method (first arrival time of P- and S-waves) by deploying dense temporary seismic stations on the southern edge of Central Alborz. Recently, the 1D model was updated in[19] by combining recorded data of all available run temporary and permanent seismic networks.

A 3D first arrival P-wave velocity model for Alborz Mountains was calculated in[20,21]. Although previous surface wave tomographic models (e.g., teleseismic model [22]; ambient seismic noise results [23,24]) can address relatively good resolution on a regional scale, they do not have required efficiency to provide an insight into the crustal structure (up to 30 km) in Central Alborz. Several P-wave tomographies (e.g. the study of Tehran City by Shirzad et al.[25]), surface wave tomography (e.g. Tehran City by [25], North Tehran and Mosha fault junction by[25]), and radial anisotropy (e.g. Tehran city by[26] and North Tehran and Mosha fault junctions by[27]) were studied providing a velocity model for a part of Central Alborz. Moreover, seismic attenuation, which generally leads the amplitude of seismic waves to decay, has been studied by several researchers[28-30] in central Alborz.

Although the method can be applied on-sites with anthropocentric earthquakes (e.g., geothermal, mine, etc.), in this paper, we use the already established data on the quality factor, velocity model, and cross-correlations of seismic activity of the region and calculate the weight field over the system under study. Then we apply a dissipative avalanche dynamics to the system by designing an activity-propagation algorithm based on the phenomenological parameters that were explored.

The paper has been organized as follows: in the next section, we present the previous analysis on the region under study and introduce the phenomenological quantities of interest. In Section "Weight Models", we introduce our model for the weight function, based on which our self-organized critical model is defined. We describe our model in Section The Dynamic Model, where the phenomenological parameters of the previous section are employed. Section "Measures and Results" is devoted to the numerical results and the activity field. The multi-fractal analysis is presented in Section "Multifractal detrended fluctuation analysis" where the generalized variance of the activity time series is analyzed, and the corresponding Hurst exponent is extracted. We close the paper with a conclusion.

Observational data, modified sandpile model

A network with evenly grid space was used in this study, which has been achieved by processing the waveforms of the occurred earthquakes in the Central Alborz region. Calculating the weight function between cells involves the interferometry of the recorded waveforms, attenuation, and seismic velocity model of the study area. The workflow of the paper is shown in Fig. 1. As is evident in this figure, we follow three objectives in this work:

Figure 1

The work flow of the paper. We start from the data set obtained by the velocity model, the Q-factor, and the weighted matrix between cells (calculated by signals-cross-correlation procedure). These quantities are employed as the input parameters of a modified sandpile model. At the end, the time series of the sandpile model is studied using the multifractal analysis of detrended time series.

1: Using the velocity model, Q-factor, and the weighted matrix between cells (calculated by signals-cross-correlation procedure), and the phenomenological parameters, we construct an effective correlated lattice.

2: We developed an effective dynamical (sandpile) model for the earthquake.

3: We apply the multifractal detrended fluctuation analysis (MDFA) to find the spectrum of the Hurst exponent of the system.

Earthquake dataset

We processed all tectonically earthquakes that occurred in the Central Alborz with M (upon which a higher concern in the community is), between 2006 and 2021 except an earthquake occurred on May 7th, 2020 (M 4.9). The reason for considering only this range is that the earthquakes with smaller magnitudes are often (but not necessarily) aftershocks, the location of which is almost a priori known. The other reason is that we do not have access to the catalog of these small events. These events were recorded by 48 seismic stations operated with three permanent seismic networks that included the following: (I) Iranian Seismological Center, IrSC, (II) Interational Institute of Earthquake Engineering and Seismology, IIEES, and (III) Tehran Disaster Management and Mitigation Organization, TDMMO. Seismic stations of IrSC and TDMMO have been equipped with SS1 (with 50 sps), CK1 (with 71.43 sps) short period, while the IIEES has been supported by Guralp CMG-3 broadband with 100 sps, respectively. From more than 3,000 events (see Fig. 2a), we selected 372 earthquakes based on following criteria: earthquakes magnitudes M 4, both horizontal and vertical location uncertainties less than 2 km, RMS 0.2 s, azimuthal gap °, and event recorded at least by 10 seismic stations.

Figure 6

The activity Field for 1000000 avalanches in the whole space, wherein each avalanche starts from a Fault site. The solid lines and black border represent the known faults and the region with good resolution, respectively.

Weight models

The size and strength of heterogeneities in the Earth’s interior can affect the coupling weight of neighborhood cells. The event interferometry approach[5] can provide the possibility of retrieving the coupling weight without a set of dense stations and/or expensive seismic imaging. For each event, we selected waveforms of the vertical component (Z-component) with signal-to-noise ratio, SNR 4.0, epicentral distance, 15 km dist 180 km, and then single station data preparation was done. This preparation includes removing mean, trend, correcting of instrument response, pre-filtering (with a 5-point zero-phase bandpass Butterworth filter) in a period range of 1–30 s, and then running time (one-bit) and frequency (whitening) domain normalizations to suppress the influence of instrument irregularities, human activities, and source time functions with nonuniform energy. The maximum amplitude of the envelope function within the expected signal window (1.5–3.5 kms−1) to the root-mean-square, rms, is the SNR definition in this study[32]. The amplitudes of the waveform parts beyond the expected signal window were transferred to zero. Next, we took the waveform part from the origin time, , to the end of the Rayleigh coda wave (1.5 kms−1) for the signals-cross-correlation procedure. Then, we applied the signal-cross-correlation (hereafter SCC) operator on the prepared waveforms of a pair-event which are recorded by a common station in the alignment of the inter-event line. The total number of inter-event raypaths depicts in Fig. 2b which is 47,864 paths. The SCC can be summarizedwhere , , , r, , and are the moment tensor solution, spatial gradient, inter-event Green’s function, coordinate vector, displacement and traction, respectively. Moreover, the indices 1 , 2, and S represent event1, event2, and Earth’s surface, respectively. Since we have taken into account only vertical components of signals, the full moment tensor solutions reduce from nine components (, , , , , , , , ) to one component (), which is a constant value. The effect of these constants ( and ) have been removed using a simple time-domain normalization operator in the data preparation step. Then, the maximum amplitude of the absolute SCC function, , was extracted and attributed to the corresponding cells along the inter-event raypath.

Given the limitations and non-uniform distribution in raypath coverage, we used a grid base tomography procedure to obtain SCC weight for all evenly grid space for the area under study. For tomography, the observed data, , is (the arrival time of ) and model, , could be the . Therefore, the relation between data and model is

The calculated data, , can be calculated using Fast Marching Method (FMM;[33,34]) grid base algorithm. When G, Green’s function, is a linear or near-linear function this formula can be

An iterative linearized damped-least squares inversion procedure can apply to minimize the observed and calculated data misfit[15,35]. Using Gauss-Newton gradient method, the relation (3) can bewhere and are the damping, and smoothing regularization parameters, and also and remark the data, and model errors, respectively, and Tr represent the transpose operator. To solve equation 4, the area under study was meshed with the grid size 14 km × 14 km, and the average of applied as initial input model, . Also, the regularization parameters, , and , were obtained by standard the L-curve by considering a trade-off between data misfit and model roughness. To stabilize the result, we used those observed data, , with residuals less than two standard deviations (2) in the inversion procedure. Figure 3a shows the obtained SCC weight map in the study area. The damping value, , was fixed to 255, and the smoothing parameter, , was 950.

Figure 3

(Color Online) (a) A two-dimensional (2D) projection of BTW lattice on the plane (). The red lines show known faults, and red circles are the events with the magnitude of M . The three-dimensional (3D) projection is shown for (b) Cross-Correlation Weight, (c) Q-factor[29], (d) P-wave velocity model[36] wherein The Cross-Correlation Weight and Q-factor are robust against .

The first arrival P-wave traveltimes can be used to obtain the 3D crustal velocity structure of the Central Alborz region. This 3D model can provide an insight into the detailed crustal velocity structure of the Central Alborz to better understand the fine-scale tectonics and seismic data transfer speed. Based on the event’s epicenter distribution and our study area, the 3D velocity model obtained by Afra et al.[36] can be an appropriate model. This model has been calculated using an iterative linearized, damped least-squares widely used inversion code SIMULPS[37]. Because the recovered anomalies of the model have been confirmed by the two different (including seismology and gravity) geophysical methods. Moreover, the reliability of this velocity model was performed by different resolution tests which are including the checkerboard resolution (both dense and sparse for checking lateral and smearing resolutions, respectively), input initial model uncertainties, events location uncertainties, and Resolution Diagonal Element[36]. In this study, the region with the fairly resolution is surrounded by a black thick border according to [36]. Fig. 3c represents this 3D velocity model from subsurface to depth of 30 km. It’s worthy to note that the 2D SCC weight was also parameterized as the depth of the velocity model as shown in Fig. 3b.

The quality factor (Q-factor) generally refers to the attenuation of the amplitude due to elastic and inelastic properties of the medium. The coda and body wave Q-factors are a powerful tool to study thermal, compositional, and deformational characteristics of Earth’s interior[38] in seismology. The body or surface waves are attenuated with rates greater than the calculated rates for geometrical spreading. The inverse of the Q-factor is composed of two factors, namely the scattering , and the attenuation along the propagation path of a seismic wave (due to the inhomogeneities within the earth and the geometrical spreading) , so that [39]. Naghavi et al.[29] obtained the Q-factor model for the Central Alborz. This model has been calculated by the Lg coda method using 1020 waveforms of the vertical component of 205 earthquakes with recorded by 35 short-period seismic stations between 2000–2009. In this study, we used the Q-factor model calculated by Naghavi et al.[29] for predicting source-receiver attenuation as depicted in Fig. 3d.

The dynamical model

Our dynamical model is a variant of dissipative continuous sandpile model[40,41], that is implemented on top of a three-dimensional cubic lattice with coordination number . The spatial extent of the region under study is the map presented in Figs. 2 and 3. We parameterized (meshed) the system so that the lattice points are fitted to existing data from which we picked the earth factors like the q-factor and the velocity filed. The resulting lattice consists of nodes in each plate parallel to the earth surface, and totally horizontally plates are considered in the perpendicular direction as shown if Fig. 3, so that the lattice includes sites. To each site of the lattice i, we attribute three intrinsic quantities obtained from the real data of the earth explained in the previous section [29,36]: the weight function (), the quality factor (), and the velocity Model (). The weight and other functions have designed so that their values for the connection between two neighboring sites is , wherein can be , , and . The resulting weight field, the Q-factor, and the velocity field model are shown in Fig. 3b, c and d. Using the weight field we construct the lattice as represented in Fig. 3a, where the position of faults and the points at which earthquakes have taken place are shown. The threshold field to be used in the dynamics of the system is proportional to the field given in Fig. 3b, i.e. . Once the network is constructed, we define the following dynamics, which is based on the local energy/stress in each site i, denoted by , taking values in the range (the summation j is over neighbors of i), so that a local status of the system is identified by the set . The initial state of the system is chosen to be random with uniform distribution.

The dynamics of the system are defined by local relaxations generated by local slipping of the fault, i.e. distributing the stress excess through the neighboring regions. The rate of stress transfer is related to the Q-factor, velocity model, and the weight of the connections. The local stimulation of a region is external and is implemented on a randomly chosen site i via , where r is a flat random number between 0 and 1), which favors a local slip of this site. This site is however static if its accumulated stress is lower than a threshold as a consequence of local static friction[42]. If exceeds , then site i is called unstable, leading to a toppling process (local relaxation), during which , where matrix is defined aswhere is the inelastic attenuation factor[38]where is an amplitude, f is frequency, and are the distance and the travel time between sites i and j. In our coarse grained model, we ignore the dependence on the frequency by setting . Moreover, to recover conservative dynamics in the limit , we set , so that the inelastic attenuation factor reduces to . The toppling rule is schematically shown in Fig. 4. Also we have km in the direction, km in the direction, and varies from 2 km to 4 km in the direction.

Figure 4

A schematic 2D set-up of the sandpile model. The rad site, i, is unstable since and sends stress to any neighboring sites, but the neighboring sites get energy.

Burst dynamics is a popular property of sandpiles, manifested by prominent avalanches which occur with a low frequency, called sometimes rare events[43-48]. Avalanches in our model are defined as a chain of local relaxations (topplings) occurring as a consequence of an external stimulus. In fact, when one site becomes unstable by an external stimulus, it topples, and as a consequence, the neighboring site might become unstable and topple in their turn, so that a chain of relaxations take place up to a time where no further unstable site is found. The duration D and the size S of avalanche are defined as the lifetime and the total number of topplings in the avalanche respectively. Then another random site for stimulation is chosen and so on. The average local stress grows almost linearly with time until reaching a stationary state after which the amount of stress that leaves the system through the boundaries is statistically equal to the number of input stress, for a good review see[49].

Measures and results

In this section, we present the results of the simulation of our modified sandpile model. The size of the lattice is fixed as explained in the previous section. We have tested five kinds of external stimuli, which are listed bellow:The case I represents the case were the local stress excess (resulting to an earthquake) takes place in a completely random region, while the case II realizes the situations where the earthquake starts on the fault. The cases III, IV, and V capture the cases where its starts from more active sites which is classified to M , M and M . For each case, we have generated over samples, i.e. the avalanches in the stationary states. The activity field for samples for the cases I, II and III are shown in Figs. 5, 6, and 7 respectively. By looking at Fig. 3a, we see a good correlation between the average activity field and the weight field. A much similarity is observed between the case I and the weight field but as the pattern of stimuli changes (cases II, III, IV, and V), the activity field show different patterns. Especially for the fault stimuli (case II) the activity in the vicinity of the fault positions is much higher, which is rather expected. A more realistic situation is the case where the stimulation takes place in the vicinity of the active points (cases III, IV, and V), e.g. Fig. 7 since these points are more active with respect to the rest regions. Apart from the observational data (the points where the earthquake has taken place), this activity field predicts the activity of the other regions, and explicitly shows important points which can be potentially the starting point of the upcoming earthquakes. Figure 8 shows a vertical cross-section of belt passing the Mosha fault, extended from along (51.50°E, 35.87°N) to (53.30°E, 35.60°N), where M.N.T.I refer to the Mosha-North Tehran fault intersection. Interestingly this shows an earthquake with M 4.9 occurred at (52.05°E, 35.78°N, 7km) on May 7th, 2020 (the blue region in the three lower graphs). This shows that our model supports the location of the next earthquake, while the input parameters (cross-correlation weight, velocity model, Q-factor, etc.) of our processing did not include the information of this earthquake.

Figure 5

The activity Field for 1000000 avalanches in the whole space, wherein each avalanche starts from a random site. The solid lines and black border indicate the known faults and the region with good resolution, respectively.

Figure 7

The activity Field for 1000000 avalanches in the whole space, wherein each avalanche starts from an event M site. The solid lines and black border indicate the known faults and the region with good resolution, respectively.

Figure 8

An activity field for the cross-section along (51.50°E, 35.87°N) to (53.30°E, 35.60°N) for the dynamics with avalanches in the whole space. All hypocenters of earthquakes within ±7 km distance from the profile are projected. The focal mechanism solution show an earthquake with M 4.9 occurred at (52.05°E, 35.78°N, 7km) on May 7th, 2020 at 20:18:21.00 which was reported by IrSC. The solid cyan, pink, and orange contours represent the P-wave velocity of 5.8, 6.1, 6.3 kms−1, respectively, calculated in Ref.[36]. The M.N.T.I abbreviation refers to the Mosha-North Tehran fault intersection (see Fig. 1).

I: completely random stimuli (the sites for external stimuli is completely random chosen),

II: fault stimuli (the sites for external stimuli are randomly chosen on the faults),

III: stimulation of low active regions, (the sites for the external stimuli are randomly chosen from sites of events with amplitudes of M ),

IV: stimulation of moderately active regions (the sites for the external stimuli are randomly chosen from sites of events with amplitudes of M ),

V: stimulation of highly active regions (the sites for the external stimuli are randomly chosen from sites of events with amplitudes of M ).

Although we applied the random event location (results depicted in Fig. 5), the potential to produce an earthquake mostly can be appeared around Kuh-Sorkh and Parchin faults. By exciting cells on the known faults (red lines in Figs. 1 and 2), this potential mostly appears around Atari-Firuzkuh, western of Mosha, Kuh-Sorkh, and Parchin faults. This potential can become apparent for cells excited by events’ location around Atari-Firuzkuh, Firuzkuh, and western of Mosha. A simple comparison for these activities (Fig. 5, 6, and 7) indicate the potential to produce an earthquake for superficial layers can be expected for micro-earthquakes (M 3.0) which is in agreement with seismicity in the study area. But, this potential can produce a larger earthquake at the greater depths (15 to 30 km) as experienced by historical earthquakes (see Fig. 2). Our results are consistent with previous studies (e.g.[36]) that P-wave tomography and gravity inversion models reveal a region at depth ranges of 12.5–17.5 km (follow the orange contour in Fig. 8 between distances of 60–90 km) for the potential of producing an earthquake with M 6.5.

The left-lateral Mosha fault, which is the most important internal fault to the central Alborz Mountains, can release a part of energy in the amount of the mm yr−1 range-parallel strike-slip motion. Vertical intersection distributions of the final activity models for Mosha fault are shown in Fig. 8. This profile along (°E, °N) to (°E, °N) is 165 km long. The random activity along this profile can clearly represent the upper crust (a thick layer up to km depth) in which most earthquakes can occur. This thickness is in agreement with the report by Abbassi et al.[18] obtained by the employment of a local seismic network. For the fault activity, the interface topography between the upper and middle crust improves realistically, which is in agreement with the bottom of event depths as reported in[50]. Inspection of this intersection reveals that high activities seem to be around the Firuzkuh, Atari faults, and also Mosha-North Tehran fault intersection (M.N.T.I) in an expectation area (see[36]) which have a potential of the future earthquake with M 6 to 6.5. As shown in activities of events profiles (8), two areas (distance 40–60 km , 100–120 km) clearly highlight most likely parts on the Mosha fault which can experience earthquakes larger than M 4. These anomalies can be consistent with the occurred events. An earthquake (52.05°E, 35.78°N, 7km) with M 4.9 has occurred on the Mosha fault around the resolved activity anomaly (distances 40–60 km) on May 7th, 2020 at 20:18:21.

The results that are shown in the above figures are graphical demonstration of the situation that the region have. It is now worthy to quantify the universal behaviors of the model for the five cases that we introduced. In Fig. 9, the distribution functions (P) of the avalanche duration D and size S are exhibited. In Fig. 9a and b, we show the log-log plot of P(S) and P(D) respectively where the linear decrease in the signature of power-law decay , . Interestingly the exponents depend of the taken situation, see TABLE 1. Noting that for the three-dimensional Bak-Tang-Weisenfeld (BTW) sandpile model, we have [44,47], we see that for all cases, the exponents and coincide well with 3D sandpile. The scaling dimension defined by is almost robust, i.e. , , , and for the cases I to V respectively.

Figure 9

The distribution function of (a) avalanche size S, (b) avalanche duration D. (c) The - plot of S-D scaling relation. (d) Activity dependent branching ratio b(x) for instantaneous avalanche size x. All plots are for the dynamics with avalanches in the whole space.

Table 1

Upper Table: The values of exponents , , , and for the dynamics with avalanches in the whole space, wherein the is defined as the width of , which is the length of the interval between two successive . Lower Table: The same for the dynamics with avalanches inside the high-resolution box.

	Case I (random)	Case II (faults)	Case III (M≥2.5)	Case IV (M≥3)	Case V (M≥4)
α¯	0.37	0.36	0.35	0.34	0.33
δα	0.46	0.45	0.42	0.42	0.41
τ_S	1.35±0.005	1.36±0.005	1.28±0.004	1.27±0.004	1.26±0.004
τ_D	1.55±0.016	1.55±0.016	1.45±0.01	1.44±0.009	1.44±0.009
α¯	0.34	0.32	0.31	0.30	0.30
δα	0.48	0.47	0.44	0.43	0.44
τ_S	1.49±0.008	1.44±0.008	1.39±0.007	1.38±0.007	1.37±0.007
τ_D	1.68±0.025	1.62±0.024	1.54±0.021	1.53±0.022	1.53±0.022

To test the criticality of the system, one can use the branching ratio function defined by the conditional expected value , where is expected value. For the critical systems or is in the vicinity of unity[4,51-54]. This function has been shown in Fig. 9d, for which , showing that although not exactly (super-critical), but we are in pretty in the vicinity of the critical region, which is confirmed by other power-law behaviors.

Multifractal detrended fluctuation analysis

Here we analyze the time series of our sandpile model and inspect the spectrum of the Hurst exponent. The method that we describe in this section is multifractal detrended fluctuation analysis (MDFA), a technique most commonly used in the analysis of self-affine fractal time series. For an uncorrelated time series the Hurst exponent is a single value [55,56]. For a mono-fractal time series, the Hurst exponent, H, is defined in terms of the asymptotic behaviour of the rescaled range. Consider a general time series , and the cumulative deviate (profile) series , where . The range is then defined by

The standard deviation is also defined as

Having defined these parameters, one can obtain the Hurst exponent H using the relation

This is applied for the case where there is an ensemble set. For the case where we have a single time series, we divide the time series into non-overlapping segments with the equal lengths s. Then the generalized variance is defined by the relation

This function shows often power-law behavior with s. For mono-fractal time series this function behaves like[57,58]

For multi-fractal time series, different Hurst exponents govern the system in different size scales. The spectrum (importantly the average and the width of the spectrum) of the Hurst exponent plays the dominant role in identifying the properties of the time series of earthquake. The multifractal analysis (MFA) is an important tool in such time series (which are not described by a single Hurst exponent and single scaling relations). Apparently, the inspecting such systems using single scaling relations, gives, at best, the average Hurst exponent, and sometimes a wrong result. Even if clean scaling relations are observed using this method, it does not say anything about the other exponents available in the system for different size scales. MFA instead gives us the spectrum of the Hurst exponent for these time series. Therefore, in this section we asses the multifractal analysis, which studies the earthquake time series (which is multifractal) with a spectrum of the Hurst exponent (for the mono-fractal time series, the spectrum is peaked with a zero width). This spectrum lets us know the mean as well as the fluctuations of the Hurst exponent, helping us to distinguish the type of correlations (and anti-correlations) existing in the system.

To describe MFA, we define a variance for each of the segments by[57]where represents the mean of Y over the segment v. The moment is then obtained using

The typical behavior of is as follows[58]where is the corresponding exponent that is related to the Hurst exponent. To extract the spectrum of the Hurst exponent we use the standard multifractal analysis. Therefore, we use the generalized variance Eq. (10), but now with a new exponentwhere represents the classical multifractal scaling exponent. This exponent is related to for stationary and normalized time series,The Legendre transform of the generalized scaling exponent gives the multifractal function as followswhere . Finally by a simple replacement one obtainswhich gives the spectrum of the Hurst exponent.

The function F(S) is shown in Fig. 10 for the cases I till V and for various amounts of , where a power-law behavior is evident in a large interval (nearly two decades). The exponents of these graphs are . If is the same for all values, then we have mono-fractal with Hurst exponent . For this case would be a peaked function around with zero width. Figure 11 shows that this is not the case, and we are facing with a strong multifractal time series for all cases. Figure 11a shows the result for which we let the avalanches go beyond the almost square region identified in Fig. 2b (high-resolution region), i.e. the region with high resolution that we are more confident about the weight field that we obtained. Figure 11b shows the results for the case where we restrict the avalanches to the square. For both cases, we see that the width of f is pretty high, and the peak position varies with the method of stimulation. Even for completely random stimulation, the peak is around (0.36) for the avalanches in the whole space (inside the high-resolution box), both being lower than 0.5, showing that the system is anticorrelated. The exponents are shown in the Table 1. Using the fact that the exponent of autocorrelation function is related to the Hurst exponent like [59] this result uncovers that the time series are anti-correlated, meaning that a large event is often followed by a small event and vice versa. This effect magnifies when the stimulation is more selective, i.e. for the fault stimuli, it is , and for highly active stimuli case () it is even smaller, for the case where avalanches are allowed to go through all the space. The main results of our analysis are gathered in Fig. 12a and b (in whole space and limited to the high-resolution box respectively), where the relation between various exponents are shown in terms of the strategy of stimuli. The decrease of the Hurst exponent (green diamond symbols) with the stimuli strategy is evident in these graphs, reflecting the fact that the time series becomes more and more anti-correlated. Along with this, the exponents and decrease, but are almost robust.

Figure 10

Calculation of the statistical function using Eq. (13) for the dynamics with avalanches in the whole space. The function of vs s display power laws , where h(q) depend on q. This feature demonstrates that the time series is a multifractal.

Figure 11

(a) Main panel: vs. for various dynamics for the dynamics with avalanches in the whole space. inset: h(q) vs. for various dynamics. (b) The same for the dynamics with avalanches inside the high-resolution box.

Figure 12

(a) Main panel: The value of different exponents vs. dynamics types for the dynamics with avalanches in the whole space. inset: The width of function in different values against various dynamics types. (b) The same for the dynamics with avalanches inside the high-resolution box.

Conclusion

In this study, we focused on the earthquakes in central Alborz, Iran. In the first part of the paper, we explored the properties of the earth in the region under study, as well as the rate of earthquakes. It helped us to construct a modified sandpile model which is much similar to the continuous dissipative sandpile model in which the energy dissipation is related to the quality factor and the velocity model of the earth. The weight function which was obtained using the signals-cross-correlation of the real seismic activities was used to estimate the weight field which was employed for distributing the energy to the neighboring sites in each toppling. Our model is based on external stimuli, the location of which can be (I) random, (II) on the faults, (III) on the low active points, (IV) on the moderately active points, and (V) on the highly active points in the region. The rate of earthquakes was shown to be related to the total activity field over the region of study. Some universal behaviors of the system are shown to be related to the scheme taken for the initial stimuli. The second part of the paper was devoted to the Multi-fractal analysis, which is exploited to extract the spectrum of the Hurst exponent of time series. The time series for each scheme was analyzed separately by multifractal analysis, for all of which the average Hurst exponent is shown to be lower than 0.5. This is an intrinsic property of anti-correlated time series, for which a large rare event is expected to be followed by a small event. The stimulation of highly active regions (in our study, the points with energies M ), a lowest average Hurst exponent is obtained, meaning that we have the strongest anti-correlated system in this case. An overall phase diagram for the model is sketched for all schemes that were considered in this paper.

A support for our model could be the location of the next earthquake in the study area on the Mosha Fault, which is in agreement with the location of the earthquake occurred on May 7th, 2020 (M 4.9). There are a lot of degrees of freedom that are ignored in our sandpile-based model, like the stochastic onset of movement (which was considered to be identified by a single threshold in each segment) and the conservation of energy/stress, which is violated in more realistic models, like the Olami-Feder-Christensen (OFC) model of earthquake[60]. Therefore, in case these simplifications fail, our model is not expected to work properly.