Unified Scaling Law for Earthquakes: Space-Time Dependent Assessment in Friuli-Venezia Giulia Region

The concept of the Unified Scaling Law for Earthquakes (USLE), which generalizes the Gutenberg-Richter relationship making use of the fractal distribution of earthquake sources in a seismic region, is applied to seismicity in the Friuli-Venezia Giulia region, FVG (Northeastern Italy) and its surroundings. In particular, the temporal variations of USLE coefficients are investigated, with the aim to get new insights in the evolving dynamics of seismicity within different tectonic domains of FVG. To this purpose, we consider all magnitude 2.0 or larger earthquakes that occurred in 1995–2019, as reported in the catalog compiled at the National Institute of Oceanography and Applied Geophysics (OGS catalog), within the territory of its homogeneous completeness. The observed variability of seismic dynamics for three sub-regions of the territory under investigation, delimited based on main geological and tectonic features, is characterized in terms of several moving averages, including: the inter-event time, τ; the cumulative Benioff strain release, Ʃ; the USLE coefficients estimated for moving six-years time intervals, and the USLE control parameter, η. We found that: 1) the USLE coefficients in FVG region are time-dependent and show up correlated; 2) the dynamical changes of τ, Ʃ, and η in the three sub-regions highlight a number of different seismic regimes; 3) seismic dynamics, prior and after the occurrence of the 1998 and 2004 Kobarid (Slovenia) strong main shocks, is characterized by different parameters in the related sub-region. The results obtained for the FVG region confirm similar analysis performed on a global scale, in advance and after the largest earthquakes worldwide. Moreover, our analysis highlights the spatially heterogeneous and non-stationary features of seismicity in the investigated territory, thus suggesting the opportunity of resorting to time-dependent estimates for improving local seismic hazard assessment. The applied methods and obtained parameters provide quantitative basis for developing suitable models and forecasting tools, toward a better characterization of future seismic hazard in the region.


INTRODUCTION
Capturing the main features of seismicity in a region is not an easy task. While average properties of earthquakes occurrence over large areas can be defined in pretty robust way, when zooming to a smaller territory, the specific features from limited data sets may prevail, and space and time regularities may be lost. Due to heuristic limitations, the smaller is the region and data time span, the less representative (and verifiable) is the resulting seismicity "portrait," with obvious consequences on earthquake hazard assessments and forecasting. Therefore in this study we aim to analyze, within the Friuli Venezia Giulia region (Northeastern Italy) and surrounding areas, a set of quantitative measures of seismic activity, which have been investigated also in other regions worldwide, so as to verify their possible similarity/generality and identify their peculiar features.
The study region is located at the Alps-Dinarides transition, along the compressional belt that marks the northern edge of Adria micro-plate; earthquakes in the region are shallow, mostly up to 10-15 km depth, and are prevalently of thrust type to the west and strike-slip to the east (Bressan et al., 2018 and references therein). Northeastern Italy is prone to large earthquakes occurrence (e.g., Gorshkov et al., 2009); the rather high seismic hazard in the region is testified by the several destructive events with occurred in the past, the most recent one being the M6.4 1976 Friuli earthquake (Slejko et al., 1999). However, the instrumental seismic activity recorded since then prevalently consists of low to moderate earthquakes, only sporadically exceeding magnitude 4.0, and with the largest recorded event having magnitude 5.6.
Several efforts have been already devoted to explore the different components of seismicity in the study region, namely the clustering properties Gentili, 2018, 2020;Varini et al., 2020) as well as the temporal variability of background seismicity component (Benali et al., 2020). Here we investigate changes in seismicity as a whole, without differentiating between clustered and background components. At the same time, the proposed analysis aims to account for the main geological and tectonic features of the territory under investigation, by assessing and cross-comparing the temporal variability of seismic dynamics within sub-regions sharing a common tectonic setting, based on seismic districts defined by Bressan et al. (2019).
In order to quantify the changes of observed seismicity in the time-space-energy domain, we make use of different parameters, including classical inter-event times and Benioff strain release (Benioff, 1951). To capture both variability and scaling properties of seismic energy release, we resort to the Unified Scaling Law for Earthquakes (USLE), which generalizes the Gutenberg-Richter (GR) relationship making use of the fractal distribution of earthquake sources in a seismic region (Kossobokov andMazhkenov, 1988, 1989;Bak et al., 2002;Kossobokov, 2020). In particular, besides the space-time variability of seismicity rate (A), earthquake magnitude exponent (B), and fractal dimension of epicenter loci (C), we consider the time series of the USLE control parameter, η (Kossobokov and Nekrasova, 2017;Kossobokov and Nekrasova, 2019).
An estimation of USLE coefficients for seismic hazard assessment at different space scales, namely from the local scale of FVG region and up to the national scale of Italy, was already performed by Nekrasova et al. (2011). This study, instead, focuses on the temporal variations of USLE parameters, rather than on their average time-independent estimates, with the aim to get new insights in the evolving dynamics of seismicity in different areas of FVG. Based on time-dependent estimates of USLE parameters, as well as on classical inter-event times and Benioff strain release, the time intervals of rather uniform seismic activity are identified for each of the investigated sub-regions; the features of seismicity, for different sub-regions and time windows, are then compared and their statistical significance is assessed by the Kolmogorov-Smirnov test. We anticipate that our analysis highlights significant differences and a number of changes in seismic regimes within the three outlined areas. Specifically, the long-term trend of seismic energy release is found to be different prior and after the occurrence of the largest main shocks within the related sub-region, evidencing non-stationary features of seismic activity that should be factored in seismic hazard assessment. Moreover the USLE coefficients in FVG subregions are found to be time-dependent and correlated, displaying interesting patterns in the dynamics of seismicity, which could be related with major earthquakes, as illustrated in some detail in Supplementary Material. Thus, the identified regions and time intervals, with homogeneous features of seismic activity, may provide valuable information for time-dependent seismic hazard assessment and may guide areas selection for earthquake forecasting.

METHODS
The seismic dynamics of the territory under investigation is characterized in terms of several moving averages, including 1) the inter-event time, τ, dual to seismic rate; 2) the Benioff strain release, Ʃ (Benioff, 1951) (i.e., the sum of square root of energy, thus, proportional to 10 0.75M ) 3) the Unified Scaling Law for Earthquake (USLE) control parameter, η and 4) the USLE parameters values estimated for the moving sexennial time intervals.
In agreement with earlier studies (e.g., Bak et al., 2002;Christensen et al., 2002) on the waiting times T between earthquakes with magnitude greater than M, occurring within a range L, the USLE generalizes the Gutenberg-Richter law (Gutenberg and Richter, 1944), and it provides the rate of occurrence N(M,L) defined as: where N(M,L) is the expected annual number of earthquakes of magnitude M within an earthquake prone area of diameter L; A, B, and C are constants. A and B characterize the annual rate of magnitude 3.5 events and the magnitude exponent respectively, which are analogous to the a-and b-values in the Gutenberg-Richter relationship (GR), while C estimates the fractal dimension df of the epicenters' loci at the site (Nekrasova and Kossobokov, 2005).
The accumulated data on hypocenters provides enough evidence for assuming that locally seismic locus generating earthquakes might have a self-similar structure of fractures of different size (e.g., Kossobokov and Mazhkenov, 1988;Bak et al., 2002). The multiscale analysis involved in evaluation of the USLE coefficients at a given site, where the count of events is performed in a set of cascading squares (i.e., a "telescope"), combines recurrences for earthquakes obtained from enveloping areas of different size, so as to get enough statistics on larger magnitude shocks from larger territories around, for a reliable confident estimation of scaling.
The results of global and regional studies (Kossobokov and Mazhkenov, 1989;Bak et al., 2002;Nekrasova and Kossobokov, 2002;Nekrosova and Kossobokov, 2003;Nekrosova and Kossobokov, 2005;Nekrasova et al., 2011;Parvez et al., 2014;Nekrasova et al., 2016;Nekrasova et al., 2015) confirm the validity of the USLE at different scales of analysis. The definition of the scales based on the preliminary analysis of the quality of data, the research goal, length dimensions of the regions under investigation and the objects of the investigation in the regions.
We can characterize the seismic rate in terms of USLE used the control parameter η: where τ is the time between two successive earthquakes, M is the magnitude of the second one, L is the distance between them, B and C the Global USLE values. Bak et al. (2002) provide a physical viewpoint on USLE: "To understand the Unified Law for Earthquakes, it is essential to see what the value of x represents. The quantity L df · S −b in the scaling function represents the average number of earthquakes per unit time, with seismic moment greater than S occurring in the area size L × L. Therefore, x is a measure of the number of earthquakes happening within a time interval T. The Unified Law states that the distribution of waiting times between earthquakes depends only on this value." Here x is equivalent to the control parameter η, and the USLE states that the distribution of inter-event times depends only on the value of the variable η. The Scaling coefficient estimation (SCE) algorithm (Nekrasova et al., 2015) is used to obtain the values of USLE coefficients. For a fine scale analysis, timely fit to the latest and best data registration period, from 1995 up to now, the OGS catalog permits calculation of the USLE parameters making use of the five bisecting steps of the spatial hierarchy, from the linear size L 0 1/2°down to L 4 1/32°. In addition, we also consider the global values of USLE (ABC GV ) parameters estimated by Nekrasova and Kossobokov (2019) to characterize the seismicity rate over the long time interval.

DATA
Earthquake data used in this study are published as monthly OGS bulletins (OGS, 2019), which are available via the following website http://www.crs.inogs.it/bollettino/RSFVG/. A general overview of data made available by the Friuli-Venezia Giulia (FVG) seismometric network, including the network characteristics, main changes related to network geometry and space-time completeness of the OGS catalog, is provided by Peruzza et al. (2015). A detailed analysis of data completeness in space and time is carried out in Peresan and Gentili (2018); based on their data analysis, we have selected the space-time and magnitude interval where data uniformly complete. Specifically, in our study, we used all earthquakes with duration magnitude M d ≥ 2.0, recorded during the time interval 1995/1/1-2019/3/21, and located within one of the three FVG sub-regions -Area 1 (A1), Area 2 (A2) and Area 3 (A3), outlined accounting for the tectonic setting of the study region, as shown in Figure 1.
The A1, A2 and A3 sub-regions are formed as combinations of 10 out of 11 FVG seismic districts as defined in Bressan et al. (2018;. The map of the 11 seismic districts of the Friuli Venezia Giulia Region is provided as Supplementary Material S1. Specifically: area A1 is the combination of ALP, CL, AMP and FOA districts; area A2 is the combination of MN, TOC and GE districts and area A3 is the combination of TAR, BOV and BA districts. The PL district is not considered in our analysis, due to its very low level of seismic activity (less than one event per year with magnitude above 2.0). The coordinates of the three subregions are provided in Supplementary Table S1.
OGS catalog reports 6640 earthquakes with duration magnitude M d ≥ 2.0 in the time span from 1995/1/1 to 2019/ 03/21. The epicenters of 3,011 seismic events fall within the selected sub-regions and 3,629 out of them ( Figure 1B). Table 1 lists the OGS catalog earthquake statistics information for the three selected FVG sub-regions and the b-value of the subregions earthquakes sets. The time series of earthquakes magnitude, latitude and longitude vs. time, within each of the three Friuli-Venezia Giulia sub-regions, are provided in Supplementary Figure S2 , which allow to verify the magnitude completeness of the data (e.g., according to Mignan and Woessner, 2012).
The spatial data distribution of events from OGS catalog and the borders of the three selected sub-regions are mapped in Figure 1B, while Figure 1C shows the empirical density distribution ρ of the OGS catalog data obtained from the regular geographic grid cells with size 1/16°× 1/16°. The temporal variability of seismicity rate parameters is compared with the occurrence of the largest earthquakes reported in the considered catalog. Although the FVG region and its surroundings experienced several destructive earthquakes in the past, including the 1976 Friuli earthquake; however, recent seismicity is characterized only by low to moderate magnitude events. The principal earthquakes occurred since 1977 are the 1998 April 12 (M d 5.6) earthquake and the 2004 July 12 (M d 5.1). Table 2 lists the information about these principal earthquakes. These events are marked as red stars in Figure 1B.
variation of seismicity within each sub-region, we consider the following parameters as a function of time: the inter-event times τ, the cumulative strain release Ʃ, the coefficients A, B, C of the USLE and USLE control parameter η. The empirical cumulative distribution functions of inter-event time τ and Benioff strain release Σ for the entire FVG Region and its three sub-regions are presented in Figure 2. Each distribution can be characterized by the c-values of the best fit trend. Specifically, the c-values of the natural exponential trend fit for τ p(>t) αexp (−γτ) are estimated as 0.48 (R 2 0.96), 0.064 (R 2 0.98), 0.048 (R 2 0.88), and 0.065 (R 2 0.83) in FVG, A1, A2, and A3, respectively. Actually, the best trend fit of τ in A3 sub-region appears to be logarithmic: p(>t) −0.091 ln(τ) + 0.41 (R 2 0.98). The c-values for the best power law trend fit for Σ, p(>Σ) αΣ −γ , are estimated as 1.23, 1.598, 1.407, 1.586 in FVG, A1, A2, and A3, respectively (R 2 0.98 for all the four cases). LogΣ is a linear function of magnitude M, therefore the distribution of Benioff strain release ( Figure 2B) differs from Gutenberg-Richter relation in the abscissa scale only. As is evident from Σ definition, which is proportional to 10 0.75M , the c-values for Σ are well in agreement with the Gutenberg-Richter the b-values (Table 1).  The sample size of the four regions considered may allow further comparison of the c-values with "q-exponential" or "qlogarithmic" functions of Non-Extensive Statistical Physics (Tsallis, 1988), as is done by (Chochlaki et al., 2018) for most of the 50 regions of the Globe. However, at this time the high R-squared values of the trend fits presented here make us to believe that in case of FVG Region and its sub-regions the Boltzmann-Gibbs theory of additivity might apply, an aspect to be further investigated in the future.  Table 3 lists the linear fit intervals of 〈τ〉 values shown in Figure 3, for each of the three FVG sub-regions; the inspection of linear fit intervals, in fact, may reveal clustered irregularity of the   Table 2) are marked with red triangles on the origin time scale. Benioff Strain Release Ʃ Time Changes Figure 4 shows the cumulative Benioff strain release vs. earthquake origin time, Ʃ(t), for the three FVG sub-regions. As in Figure 3 the origin time of the principal earthquakes is marked by red triangles. We may observe that, while the curve associated with A1 (blue) and A2 (green) sub-regions increase steadily, the red curve associated to the A3 sub-region displays evident jumps corresponding to both principal shocks, mostly reflecting the increased seismicity during aftershocks sequences. Table 4 lists the linear fit intervals of Ʃ(t) function for the three FVG sub-regions, shown in Figure 4. The time intervals listed in Table 4 are determined following the same procedure applied for 〈τ〉; the start/end times correspond to the first/last earthquakes that bound the intervals providing the best linear fit for Ʃ(t). Specifically, the sub-region A1 could be characterized by a uniform trend during the interval from 1997/03/05 to 2019/ 03/19 (i.e., up to the last event reported in sub-region), with Ʃ-value rate change (expressed in dimensionless arbitrary units, a. u.) of 2.027 × 10 3 a. u. per day (R 2 0.997). The A2 sub-region could be characterized by two intervals with uniform trend, from 1997/03/08 to 2001/08/18 and from 2003/05/01 to 2019/01/15 (last event in sub-region), which are characterized by Ʃ-value rate changes equal to 3.144 × 10 3 a. u. per day and R 2 0.984, and 0.902 × 10 3 a. u. per day and R 2 0.983 respectively. Similarly, for the A3 sub-region two intervals could be identified, from 1999/ 11/01 to 2004/04/04 and from 2006/04/01 to 2019/02/28 (i.e., up to the last event reported within the sub-region), with Ʃ-value rate changes equal to 2.204 × 10 3 a. u. per day and R 2 0.967 and 0.905 × 10 3 a. u. per day and R 2 0.992 respectively.

Inter-Event Time τ Temporal Variability
One can see that a stable rate of Ʃ-value over the observed time interval is obtained when no large event occurs (A1 subregion), while in case of a relatively large event, the rate before and after it displays essential differences. Specifically, the M4.9 Sernio earthquake (2002/02/14; Lat: 46.426°N Lon: 13.100°E) eventually marks the separation between the two linear seismic rate intervals for A2 sub-region; the ratio of the Ʃ(t) slope coefficients is 3.5 for A2 two intervals. In a similar way, the 2004 Kobarid event separates two linear seismic rate intervals for A3 sub-region;  the ratio of the Ʃ(t) slope coefficients is 2.5 for A3 two time intervals.
Given that Ʃ(t) function within sub-region A1 is characterized by a uniform trend during the whole investigated time span, we compare its slope coefficient with those of the other two regions, A2 and A3, for the different time intervals. Specifically, the ratios of A1 slope coefficient vs. slope A2 and A3 in advance of large events are 1.6 and 1.1, respectively, while after the large events they are both equal to 0.45 (i.e. the slope in A2 and A3 is 2.2 times lower than the slope in A1). Thus, after the large events, the increase of seismic rates in A2 and A3 sub-regions slows down, compared to the steady seismic rate in A1 sub-region.
In this analysis, the slope of Ʃ(t) linear fits essentially characterize the long-term steady trend of seismicity within each area, excluding pre-shock and aftershock related changes. Following Vallianatos and Chatzopoulos (2018), the generalized Benioff strain evolution (which here corresponds to the well-known cumulative Benioff strain, with exponent 1/2), during the initial part of a main shock preparation process is linear; then, as the time of the earthquake approaches, it deviates from linearity due to the beginning of an accelerating deformation stage. Our analysis shows that the slope of the linear part of cumulative Benioff strain release in areas A2 and A3 changes significantly after the occurrence of the largest earthquakes. According to a nonextensive statistical physics view (Vallianatos and Chatzopoulos, 2018), this observation suggests that the common exponent m associated with the steady (normal) time variation of the generalized Benioff strain within an area, may change after large events.

Unified Scaling Law for Earthquakes Coefficients Space-Time Variability
The available OGS catalog data allowed us to obtain 150 reliable estimates of the medium-term (24.2 years) USLE parameters. Specifically, we considered the hierarchy of square boxes, with linear size equal to 1/2°, 1/4°, 1/8°, 1/16°and 1/32°, respectively, centered at the nodes of a regular grid with 1/16°spacing, which include at least eight earthquakes from the OGS catalog in their 1/ 16°× 1/16°vicinity. The obtained A, B, and C coefficients are mapped in Figure 5, along with the squared sum of their standard errors σ A , σ B , and σ C . The error of determination of the USLE coefficients does not exceed 0.05, which confirms a rather high quality of the mapped values for the FVG territory. Specifically, there are 58 cells located within sub-region A1, 40 cells within sub-region A2, and 52 cells within A3, while 10 cells are out of three selected sub-regions (Table 1).
FIGURE 5 | The regional maps of A, B, and C coefficients (top panels) and spatial distribution of the sum of standard errors σ 2 σ 2 A + σ 2 B + σ 2 C . Color code symbols correspond to estimates associated with cells of size 1/16°× 1/16°. Frontiers in Earth Science | www.frontiersin.org January 2021 | Volume 8 | Article 590724 The seismic activity distribution (coefficient A), normalized to recurrence of a magnitude 3.5 earthquake, in a unit area of 1°× 1°a nd in a unit time of one year, varies mainly in the interval [−0.25; 0.45] for A1 sub-region, [0.22; 0.42] for A2 sub-region, and [−0.10; 0.45] for A3 sub-region, with median values 0.10, 0.36 and 0.33 for A1, A2 and A3, respectively. The B values, which characterize the slope of the frequency-magnitude graph, for A1 sub-region vary from 0.8 to 1.2, without any dominant value. For A2 sub-region B-values are well focused: they are equal to 1.1 in 45% cases, and to 1.2 in 55% cases. Within A3 sub-region B-values spread from 0.9 to 1.3.
The fractal dimension of spatial distribution of epicentres C in A1 sub-region has a sharp peak (55% of cells) around the value C 1.19; within A2 sub-region C varies from 0.90 to 1.47 (50% of cells), with the median value 1.40; in the A3 sub-region C displays a large variability, from below 0.9 to above 1.25. Figure 6A compares the empirical cumulative distributions functions F i (coef) of A, B and C coefficients, estimated over the fine scale grid for the three FVG sub-regions. Figure 6B shows the plots of the pairwise difference curves F m (coef) − F n (coef) for each USLE coefficient, and allows better understanding the differences in seismic activity in the three sub-regions. For instance, we may observe the remarkable differences between the distributions of C parameter obtained for A1, A2 and A3 sub-regions. Specifically the C-value distributions for A1 corresponds to a (fault) zone with a common dominant Alpine trend, while in A3 it corresponds to linear set of clusters, aligned with a Dinaric trend; finally it characterizes A2 as a highly complex, fractured (fault) zone, located at the Junction of the Alpine and Dinaric (fault) systems.
The maximum absolute difference between the empirical distributions is commonly used in the Kolmogorov-Smirnov two-sample criterion to distinguish whether or not the values from the two samples are drawn from the same statistical distribution of independent variables. Here the two sample Kolmogorov-Smirnov statistic (Smirnov, 1948) λ K−S is applied to pairwise sets of USLE coefficients distributions, as estimated for each of the different FVG sub-regions. λ K−S is defined as: λ K−S (D, n, m) [n × m/(n + m)] 1/2 × D, where D max |F m (coef) − F n (coef)| is the maximum value of the absolute difference between the empirical distributions F m (coef) and F n (coef), coef A, B, C, whose sample sizes are n and m respectively. Table 5 summarizes the results of comparison for each pairwise set of coefficients, in terms of D and λ K−S . It is possible to observe that, with 95%  probability, the three sub-regions are different in terms of USLE coefficients. The λ K−S statistic for F(A) allows us to conclude that the distributions of coefficient A for A1 and A2 sub-regions, as well as those for A1 and A3 sub-regions, are significantly different, with probability larger than 99.9%, The A coefficient distributions for A2 vs. A3 can also be marked as different, with probability larger than 95%. The pairwise comparison with λ K−S statistic for F(B) distributions provides similar results. Specifically the F(B) distributions difference for A1 vs. A2, A1 vs. A3 and A2 vs. A3 confirm that they are all significantly different, with more than 99.9% probability for A1 vs. A2, and 95% probability for the others. For the C-coefficient distributions, the two samples statistics for three pairwise differences allows us rejecting the assumption that coefficients follow the same statistical distribution with more than 99.999% probability.
Besides the medium-term average estimates described so far, the time-variable estimates of USLE coefficients have been performed for the fine grid, considering moving six-years time intervals with one year shift. As a results 19 sets ( Figures S3.1c,S3.3c), while the A-value slightly increases in sub-region A1. Accordingly, the number of low magnitude earthquakes decreases during the 2011-2019 all over the FVG territory, in A3 sub-region especially. C-values in the three sub-regions vary from one to less than 1.4 and, from 2000 to 2010 years, take relatively low values in A2 and A3 sub-regions due to the aftershocks effect; in addition, an interesting decreasing trend of C-values, down to 0.8, characterizes the A1 region from 2013 to 2019 years.
Finally, to investigate in more detail the space-time variations of USLE coefficients, the medians of A-values, B-values and C-values have been computed, for six-years time intervals with shift of one year, within each of nine out of 10 seismic districts of Friuli-Venezia Giulia Region (Supplementary Figure S1), as defined by Bressan et al. (2018;. The PL and FOA districts are not considered, because they do not include enough events for reliable estimation of USLE parameters. The values of the medians obtained for each of the nine seismic districts (Supplementary Figure S3.4) support the grouping into the three sub-regions used in this study: the districts that compose A3 sub-region (TAR, BA, BOV) display temporal trends of the A-B-and C-values much different that the other districts.

Control Parameter η Temporal Variability
The main parameter, which allows characterizing the seismic rate in terms of USLE, is the control parameter η defined in Eq. 2. In order to compute the η-values, necessary to investigate the temporal variability of this parameter within A1, A2 and A3 sub-regions, we used the global values of USLE coefficients listed in Table 1. Accordingly, the A, B and C coefficients used in Eq. 2 correspond to longterm robust estimates of USLE, which account for the moderate-large M > 4.0 earthquakes that are reported in USGS Global Hypocenter Data base system, 1964-2001. Note that for the purpose of this analysis we have extended the time window, including an additional time interval from 1988 to 1994, so as to be able and analyze seismicity changes before the principal 1998 Kobarid earthquake; according to Peruzza et al. (2015) this time interval is characterized by uniform OGS network conditions, and thus can be used for our analysis. Figure 7 shows the plots of the 50-events moving average of 〈η〉 vs. origin time, for the earthquakes located within the A1, A2 and A3 sub-regions, which occurred during the time interval from January 1988 to April 2019. It is possible to observe some nearly-flat portions of the graphs, where the moving averages vary within one decimal order; these periods of stability are characterized by a low rate of seismicity. To quantify this feature, let us define the periods of stability by the condition that {t: 〈η〉(t) is larger than 〈η〉 max /10}, where 〈η〉 max is the maximal value of 〈η〉 over the entire considered time interval. Table 6 lists the periods of 〈η〉 stability identified for the three FVG sub-regions. For the A3 sub-region we identified three periods of stability of 〈η〉 interrupted by bursts of activity, associated with the origin times of principal earthquakes and their aftershocks. Namely, before the Kobarid April 12, 1998 event, a 5.8 years time interval with very low seismic rate 〈η〉 10.02 ± 0.39 is detected from 1992/07/13 to 1998/04/12; between the two principal events, precisely from 2001/10/26 to 2004/07/ 13, a 2.8 years interval with much larger seismic rate 〈η〉 3.59 ± 0.11 has been identified; during the last 13.4 years interval, from 2005/10/02 to the end of the OGS catalog, an intermediate level seismic rate, with 〈η〉 6.92 ± 0.20, is quantified. For A1 subregion we determined four periods of 〈η〉 stability: a low seismic rate interval with 〈η〉 7. 87 ± 0.49 in 1990/09/03-1996/06/29; two comparatively higher rate intervals in 1998/02/07-2000/09/ 19 and 2001/03/18-2003/12/27, with 〈η〉 equal 2.71 ± 0.05 and 3.31 ± 0.11 respectively; a 15 years interval of moderate seismic rate, with 〈η〉 value 4.71 ± 0.08, from 2004/03/24 to the end of catalog. Finally, according to η control parameter, a single 〈η〉 stable interval can be identified for A2 sub-region, from 1988/09/ 18 to 2019/01/15 (30.3 years), with 〈η〉 value 2.51 ± 0.05, which apparently characterizes the relatively high and constant seismic rate of this area.  The discrete and cumulative empirical distribution functions of η for each of three sub-regions are compared (left top panel, Figure 8 and central top panel, Figure 8), confirming a broad spread over about six or more decimal orders of the bulk density distribution of η. According to the nonparametric two-sample Kolmogorov-Smirnov statistic λ K-S , applied to the pairwise differences of η distribution functions η(A3), η(A2) and η(A3) for the three sub-regions, all of them have different probability distribution, with probability above 99% (Tabure 7). Finally, the comparison is restricted to the most recent time interval, from 2005/10/02 to the end of OGS catalog, when according to Table 6, stable intervals of 〈η〉 are obtained in all FVG sub-regions ( Figure 8, bottom panels). The pairwise differences of η distribution functions for 2005-2019 time interval (Figure 8, right bottom panel) and the λ K-S values for pairwise differences of η distributions (Tab. 7), confirm that even for the time period without principal earthquakes η(A3), η(A2) and η(A3) do not come from the same distribution, with probability above 99% in almost all cases (except for η(A2) vs. η(A3), for which probability is still above 95%).

DISCUSSION AND CONCLUSIONS
In order to get new insights in the evolving dynamics of seismicity, which controls its temporal variations within different areas of FVG territory, a multi-parametric analysis was performed. Different parameters, including inter-event times τ, Benioff strain release Ʃ, and the control parameter η of the USLE, were used to assess heterogeneity in spatial and temporal patterns of earthquakes occurrence, as reported in OGS catalog for Northeastern Italy and its surroundings.
As a preliminary step, the space variability of the USLE coefficients, namely the seismicity rate (A), the earthquake magnitude exponent (B), and the fractal dimension of epicenter loci (C), was examined in some detail, comparing the values obtained within three sub-regions defined grouping seismic districts with similar tectonic features (Bressan et al., 2019). Different tectonic domains, in fact, can be identified within the study region, which is located at the junction between the E-W oriented Alpine and the NW-SE oriented Dinaric fault systems ( Figure 1A). According to Bressan et al. (2018), the superposition of different tectonic phases caused high fragmentation and heterogeneity of the upper crust in this area, which is characterized by sharp lateral heterogeneities of the elastic moduli in the upper crustal structure (down to 10 km depth). The analysis performed in this study confirms that the features of seismic energy release, including their temporal variations, are statistically different within the three outlined subregions A1, A2 and A3. The spatial distribution of USLE coefficients is shown in Figure 5, while their empirical probability functions are given in Figure 6. Specifically, the coefficient of magnitude balance B, which is analogous to the b-value of the Gutenberg-Richter law and quantifies the relative proportion between small and large earthquakes, is characterized by well focused B values within [1.1-1.2], in the central zone (area A2); this area hosted the M6.4 1976 Friuli earthquake, but since then instrumental seismicity has been characterized only by moderate size events. At the same time, the western zone (area A1), corresponding to the Carnie Pre-Alps, is characterized by higher variability, which indicates higher spatial heterogeneity in the seismic energy release, and by low-intermediate B values, within the range [0.8-1.2], and thus by a comparatively low ratio of small-moderate magnitude events. Instead, the Dinaric eastern zone (area A3) is associated with relatively higher B values, in the range [0.9-1.3], evidencing the occurrence of a comparatively large number small to moderate seismic events, mostly related with aftershocks occurrence. Significant differences have been observed between the distributions of the fractal dimension of earthquake epicenters C ( Figure 5), obtained for A1, A2 and A3 sub-regions ( Figure 1A). Specifically, within area A3, C appears characterized by quite low values, mostly up to 1.0, which can be related with a linear set of clusters, aligned with a Dinaric trend. This is possibly due to the occurrence, within area A3 of the two largest magnitude events reported in the considered data set (the 1998 and 2004 Kobarid earthquakes) that, along with their highly clustered aftershocks, dominate seismicity in the area. The western area (area A1) displays prevailing intermediate values, corresponding to a (fault) zone with a common dominant Alpine trend and associated with rather complex swarm-like earthquake sequences (Peresan and Gentili, 2018). Finally, the central zone (area A2) is characterized by higher fractal values, corresponding to a highly complex, fractured (fault) zone, located at the junction of the Alpine and Dinaric (fault) systems ( Figure 1A).
In the long-term of the considered dataset (i.e., from 1988 to 2019) we found different intervals of rather steady seismic activity, which are characterized by a near constant value of η, with switches at times of transition associated with the relatively large M d > 5.0 events. As long as the temporal features of seismicity are concerned, a time interval of rather stable seismic activity could be determined according to the different parameters; during such interval, starting on 2005 and up 2019 (i.e., up to the end of available data) no major earthquakes (i.e., M d > 5.0) are reported in the catalog, thus providing presumably independent (background) parameters values. Although the temporal pattern of activity rate changes identified in this study reflects trends at the sub-regions scale, this result appears compatible with the observation of a period of rather low background seismicity rate (i.e. for the declustered catalog) ongoing since more than a decade, detected by Benali et al. (2020) for the whole FVG region.
The results obtained for Northeastern Italy and surrounding areas confirm similar analysis performed on a global scale, in advance and after the largest earthquakes worldwide. Specifically, we found that: 1) the dynamical changes of τ, Ʃ, and η in the three sub-regions highlight a number of different seismic regimes; 2) the seismic activity prior and after the occurrence of strong main shocks (e.g., the 1998 and 2004 Kobarid earthquakes) is characterized by significantly different parameters within the related sub-region; 3) the USLE coefficients in FVG region are time-dependent (as observed in Nekrasova and Kossobokov, 2005;Nekrasova, 2007;Nekrasova et al., 2011) and show up correlated, displaying interesting features in dynamics of seismicity that can be related with major earthquakes (see Supplementary Figures S3.1-S3.3).
The temporal changes of the USLE coefficients estimated for three FVG sub-regions exposed correlated, though complex behaviors in dynamics of the Earth crust hierarchical system of blocks-and-faults. In addition, the analysis of time variations of the cumulative Benioff strain release, Ʃ(t), evidenced that the slope of its linear long-term trend may change significantly after the occurrence of a major earthquake. Although the number of the moderate earthquakes in the FVG region is too small for contributing to "the hypothesis that many large earthquakes are preceded by accelerating-decelerating seismic release rates which are described by a power law time to failure relation" (Vallianatos and Chatzopoulos, 2018), according to a non-extensive statistical physics view, our observations suggest that the steady (normal) state of the system (as described by the common exponent m associated with the time variation of the Benioff strain within an area), may change after the occurrence of a large earthquake.
The obtained results highlight non-stationarity of seismic activity, at a time-scale of several years and up to decades, in agreement with earlier findings by Benali et al. (2020), an element that should be taken into account for improving local seismic hazard assessment. The regions and time intervals identified in this study, which display homogeneous features of seismic activity, may supply valuable information toward timedependent seismic hazard assessment (e.g., Kossobokov et al., 2015 and references therein), while providing new constraints for earthquakes forecasting in Northeastern Italy.

DATA AVAILABILITY STATEMENT
Publicly available datasets were analyzed in this study. The earthquake datasets used in this study are published as yearly and monthly bulletins by the National Institute of Oceanography and Applied Geophysics (OGS) and are publicly available via the OGS website (http://www.crs. inogs.it/bollettino/RSFVG/).

AUTHOR CONTRIBUTIONS
AN performed time series computations, and ABC coefficients analysis, while AP provided the data, performed preliminary data analysis and defined the setting and regions for investigation. Both AN and AP collaborated at data interpretation and manuscript write up.

FUNDING
The IEPT RAS researcher part of this study is supported by grant of the Russian Science Foundation .