Statistical Investigation of the 2020–2023 Micro-Seismicity in Enguri Area (Georgia)

Abstract

In this study, we analyzed the microearthquake seismicity in the Enguri area (Georgia) recorded between 2020 and 2023 using a newly installed seismic network developed within the DAMAST project. The high sensitivity of the network allowed the detection of even very small seismic events, enabling a detailed investigation of the temporal dynamics of local seismicity. Statistical analyses suggest that the seismic activity around the Enguri Dam is influenced by a combination of natural tectonic processes and subtle reservoir-induced stress changes. While the dam does not appear to exert strong seismic forcing, the observed ≈7-month delay between water level variations and seismicity may indicate a triggering effect. Localized stress variations and temporal clustering further support the hypothesis that water level fluctuations modulate seismic activity. Additionally, the mild persistence in interoccurrence times is consistent with a stress accumulation and delayed triggering mechanism associated with reservoir loading.

1. Introduction

The Enguri Dam is situated along the Enguri River, which flows through the northwest part of Georgian territories, including the region of Abkhazia, and plays a crucial role in the country’s power supply. It is situated in the northwestern part of Georgia, within the southern section of the Greater Caucasus Mountains, where the ongoing convergence between the Eurasian and African–Arabian plates results in folding, thrusting, and crustal shortening [1,2,3,4,5]. Previous studies have extensively examined the seismotectonic characteristics of the Enguri area, including GPS-derived convergence rates [1,6,7,8,9,10], regional fault systems [11,12,13], and seismicity patterns [14,15].

Standing at 271 m, the Enguri Dam is among the tallest arch dams in the world. Its reservoir, with an mean volume of 1 billion cubic meters and yearly water level fluctuations of approximately 100 m, is situated in close proximity to active faults. This unique setting makes it a compelling case study for reservoir-triggered seismicity and the complex interactions between dam reservoirs and fault systems [16]. According to Basili et al. [17], the active fault structures around the Enguri Dam pose the most significant seismic hazard to the Enguri Hydropower Plant (HPP).

In the framework of the DAMAST (Dams and Induced Seismicity Transfer for Sustainability) project (https://www.damast-caucasus.de/68.php), a local modern seismic network was installed to monitor the seismic activity around the Enguri Dam. The newly established local seismic network consists of four surface stations, three stations located in shallow boreholes and one in a deep borehole. Such significant advancements in the local seismic network, developed within the DAMAST project, enhanced the detection capabilities, allowing even minor seismic events to be recorded. This improvement enables a more detailed examination of the catalog’s properties. Although the recorded seismicity shows activity along existing faults, ongoing studies are focused on the detailed identification of active fault branches, their geometry, and slip mechanisms, as well as exploring the potential connection between seismic activity and fluctuations in the dam-reservoir water levels.

In this work, we aim to provide a comprehensive statistical analysis of the temporal and magnitude distribution characteristics of the seismicity of the Enguri area recorded from 2020 to 2023 within the DAMAST project. In addition to well-established seismological and statistical metrics—such as the b-value of the Gutenberg–Richter law, the magnitude of completeness, and the coefficient of variation (both global and local)—this study examines the clustering in the time distribution of the Enguri seismicity by using the Allan factor analysis, as well as the correlation properties of the interevent times and magnitude by means of the detrended fluctuation analysis, and, moreover, introduces a novel set of parameters to characterize the seismicity of the region under study. These parameters are derived from the visibility graph (VG) method, originally developed by Lacasa et al. [18], and then increasingly utilized in recent studies to analyze seismic sequences [19,20].

Our study focuses also on the correlation between the seismicity and water level of the Enguri Dam. Information on water level fluctuations at the Enguri reservoir was obtained in collaboration with LTD Engurhesi, the operator of the Enguri Hydropower Plant (HPP). The dataset consists of daily reservoir elevation records, systematically collected since the dam’s commissioning in the 1970s and continuously maintained to the present. For this study, we focused on the 2020–2023 period to facilitate comparison with the local seismicity catalog.

2. Seismo-Tectonic Settings

The Enguri Dam is located directly on the Mukhuri segment (MS) of the reverse Racha-Lechkhumi fault system, which extends with a strike of 288°. The average (mean) maximum magnitude for this fault is estimated at 6.61, with a mean slip rate of 0.375 mm/year. The largest historical earthquake, the Kvira event (Mw = 6.9; 1250), occurred along this fault. The magnitude of this earthquake corresponds closely to the 98th percentile estimate of the maximum magnitude (MWMAXP98 = 6.92) for the Mukhuri segment, meaning there is only a 2% probability that the actual maximum magnitude exceeds this value. To the south, the Enguri Dam area is surrounded by three segments of the Odishi active fault system. The Achigvara (OA) segment (strike 350°) and the Martvili (OM) segment (strike 233°) represent right-lateral and left-lateral faults, respectively, while the Khobi (OKH)segment (strike 298°) is characterized as a reverse fault. The Odishi fault has an average maximum magnitude of 6.43 and a slip rate of 0.1 mm/year. At the junction of the Achigvara and Khobi faults, a historical earthquake with Mw = 6.1 occurred in the year 1614. Another active fault that could significantly impact the Enguri HPP is the Larakvakva reverse fault, which consists of two branches—North and South—crossing the Jvari Reservoir. The Larakvakva South branch (LS) has a strike of 288° and an average maximum magnitude (Mwmax) of 6.61. A historical earthquake, the 1750 Akiba event (Mw = 7.0), occurred along this fault. The magnitude of this earthquake closely corresponds to the 98th percentile estimate of the maximum magnitude (MWMAXP98 = 7.02). The average slip rate of the fault is 0.17 mm/year, with the North branch (LN) being more active, exhibiting an average Mw of 6.74 and a slip rate of 0.34 mm/year.

Figure 1 presents the surface projection of active faults around the Enguri Dam [17], along with the investigated seismicity (yellow dots) recorded by the new local seismic network developed within the DAMAST project within 33 km radius around the Enguri Dam and down to 30 km depth. The dataset covers the period from 1 October 2020 to 31 December 2023, and the earthquake magnitudes range from −1.3 to 3.3. Seismic activity is notably high along the Mukhuri segment, affecting both branches of the Larakvakva reverse fault and the Achigvara segment. Additionally, significant seismicity is observed along the pressure tunnel that connects the Jvari reservoir to the powerhouse at the Gali Reservoir (Abkhazia), located approximately 15 km southwest of the Enguri High Dam. Old literature published by Murusidze in 1980 [21] references the so-called ‘Enguri faults,’ marked by red dashed lines in Figure 1, which are clearly visible in the topography. However, these faults are not included in the European active fault database [17], primarily due to insufficient investigation. Based on the seismic activity pattern, we can infer that the Ingrishi faults are quite active and may extend farther southwest toward the Gali Reservoir (GR).

3. Seismic Network and Data

3.1. Description of the Seismic Network

A ten-station local seismic network was established to monitor low-magnitude seismic activity near the Enguri Dam. Seven of the stations are located within a 10 km radius of the dam, while the remaining three are within a 10 km radius of the Gali Reservoir, ensuring comprehensive coverage with enhanced sensitivity near active fault segments, as illustrated in Figure 1.

In this study, we used a dedicated seismic catalog independently constructed from continuous waveform data recorded by the local seismic network deployed around the Enguri Dam within the DAMAST project. This catalog was not derived from any official or national database. The network, specifically designed to monitor microseismicity near the dam, enables the detection of events with magnitudes well below the threshold of the national seismic monitoring system [15].

The finalized configuration of the deployed network around the Enguri Dam, comprising both surface and borehole stations, includes Centaur data loggers, routers, and SIM cards for data transmission at each active station. The network was deployed in five major phases: October 2020, December 2020, September 2021, January 2022, and October 2023. The initial deployment consisted of five stations, followed by the addition of two more. Subsequently, one station was installed in the Abkhazia region, with two additional stations added later. The coordinates and elevations of these stations are provided in Table 1.

Among the stations, KIT1 and HPP are powered by AC electricity, while the others rely on solar panels combined with battery storage. Most stations (KIT1, NIKA, BUFF, BRID, KETI, and DOG) have maintained continuous data recording and transmission since their installation. Site selection was significantly influenced by challenging terrain on the western side of the dam reservoir and restrictions near the Georgia–Abkhazia border. Despite these obstacles, the GULB station was successfully installed in a high-altitude mountainous area. However, due to unreliable real-time data transmission and extremely difficult—sometimes impossible—access, the GULB station was decommissioned in June 2023.

Borehole stations are equipped with Isometrics Trillium Compact Posthole 20s (TC20PH2) broadband sensors installed at depths of 17–20 m (BUFF: 19 m, NIKA: 17 m, DOG: 20 m) and 250 m (KIT1), effectively minimizing surface noise. These sensors provide a flat velocity response over a broad frequency range, from approximately 0.008 Hz to beyond 100 Hz, in accordance with manufacturer specifications and confirmed by instrument response spectra analysis. Their self-noise levels approach the New Low Noise Model (NLNM) within the most sensitive portion of the seismic band, making them ideal for low-noise borehole environments. Additionally, they offer a dynamic range exceeding 156 dB at 1 Hz. Data acquisition is handled by Centaur CTR4-3 digitizers, which deliver 24-bit resolution, GPS-based timing, and onboard data storage.

Surface broadband stations (KETI, BRID, and GULB) are equipped with Kinemetrics MBB-2 sensors connected to Centaur digitizers. These sensors exhibit a flat velocity response between approximately 0.03 Hz and 80 Hz, as confirmed by pole-zero analysis. Although they are nominally specified down to 0.0083 Hz (120 s), the effective usable bandwidth is more limited in practice. Their dynamic range exceeds 155 dB at 1 Hz, and careful mechanical coupling during field installation helps minimize noise.

Short-period installations (e.g., OKM, HPP, and GAL), especially in remote or access-restricted areas, employ three-component 4.5 Hz geophones (type PE-6/B) connected to compact DATA-CUBE³ Type 2 recorders. These sensors provide a flat velocity response above 4.5 Hz, with a usable bandwidth extending up to 80 Hz. The geophone response spectra were modeled using simplified empirical two-pole transfer functions, validated through SEED RESP files. This configuration offers a dynamic range of approximately 100–110 dB, ensuring high sensitivity for detecting shallow microseismic events.

The instrument responses at all stations were verified using SEED RESP files, ensuring consistent response spectra across each component (HHE, HHN, and HHZ) of the same sensor model. These responses include the complete pole-zero structure, gain, and normalization factors, and were uniformly applied during data correction and processing. A summary is provided in Table 2.

3.2. Seismic Ambient Noise Analysis and Station Performance

To evaluate the ambient noise environment at each seismic station, we analyzed the root-mean-square (Vrms) velocity amplitudes of both vertical and horizontal components by integrating the power spectral density (PSD). Two frequency bands were selected to characterize the noise: 1–10 Hz, which captures the dominant energy of local and regional events and is critical for local magnitude ($M_L$) estimation, and 5–15 Hz, which highlights high-frequency noise that may impact the detectability of small-magnitude ($M<2$) events.

The Vrms values (see Table 3) exhibit considerable variability across the network. Borehole stations—including KIT1 and the shallower boreholes BUFF, NIKA, and DOG—consistently show the lowest Vrms levels, confirming their superior performance for microseismic monitoring and reliable magnitude calibration.

The surface broadband station GULB, located in a high-altitude alpine environment, also demonstrated excellent noise performance prior to its removal in mid-2023, providing valuable high-quality data during its operational period. In contrast, surface stations such as KETI and GAL exhibit elevated Vrms levels, with GAL particularly noisy in the 5–15 Hz band—likely due to local anthropogenic activity or wind-induced noise. The OKM and HPP stations show intermediate noise levels, likely benefiting from more favorable site conditions.

To further evaluate station performance, we calculated the median vertical-component power spectral density (PPSD) [22] for each station using the ObsPy PPSD module. The resulting spectra were plotted against the New High Noise Model (NHNM) and New Low Noise Model (NLNM) [23] defined by Peterson (1993) [24], which serve as the global reference standards for seismic noise levels (Figure 2). PPSD values are reported in decibels (dB) relative to $\left[{\mathrm{m}}^2/{\mathrm{s}}^4/\mathrm{Hz}\right]$, representing acceleration power spectral density.

The PPSD analysis confirms and complements the V_rms-based observations: stations KIT1, BUFF, NIKA, DOG, and GULB exhibit spectral noise levels that closely approach the New Low Noise Model (NLNM) across the entire seismic frequency band. In contrast, the surface stations GAL, KETI, OKM, and HPP show elevated noise levels over a broad portion of the spectrum (until 10 s), reflecting more exposed or less shielded installation environments. Nevertheless, all stations remain below the New High Noise Model (NHNM), confirming that even the noisiest sites in the Enguri network operate within globally accepted noise limits.

Together, the Vrms and PPSD analyses provide a coherent and detailed characterization of ambient noise variability across the Enguri seismic network, underpinning data quality assessment, magnitude completeness evaluation, and microseismic detectability throughout the monitored area.

Although ambient noise studies commonly rely on power spectral density (PSD) values in acceleration units, root-mean-square velocity (Vrms) estimates in $\mu$m/s provide an intuitive and practical complement for evaluating station performance, particularly within local networks. Similar methods have been employed to assess broadband station quality and sensor performance [25,26,27]. While formal Vrms-based noise thresholds are not yet standardized in the literature, the values reported in this study fall within the expected range for well-installed broadband and borehole sensors. This analytical framework is consistent with widely adopted methodologies for ambient noise benchmarking and instrumental self-noise characterization [28].

3.3. Local Magnitude Estimation and Uncertainty

The local magnitude (ML) for each event was calculated using the SEISAN software, which applies a standard Richter-type formulation [29] and was later adapted for regional seismic networks [30]. The implementation used in this study follows the methods described in [31], as part of the SEISAN processing framework. Although this method is not specifically calibrated for the Enguri region, it was applied consistently across all events. This uniform application ensures that any systematic bias in absolute magnitudes remains constant, thereby preserving the internal relative scaling essential for statistical analyses such as Gutenberg–Richter b-value estimation, Allan factor calculations, and time clustering metrics.

For each event, ML was independently computed at each station that recorded it. The final event magnitude was defined as the average of the station-specific ML values. The associated magnitude uncertainty was quantified as the standard deviation of these station estimates, reflecting intra-event variability arising from local noise conditions and site amplification effects. Across the analyzed events, this standard deviation typically ranged from 0.1 to 0.3, with a network-wide average uncertainty of approximately 0.23 magnitude units.

To ensure robustness in magnitude determination, only events recorded by at least four stations were included. The number of contributing stations varied between four and eight, depending on event magnitude, location, and signal-to-noise ratio.

Ground conditions—such as surface versus borehole installations and soil versus rock coupling—can significantly affect recorded amplitudes. No site-specific corrections were applied prior to ML estimation. Instead, to mitigate biases from individual stations with potentially amplified ground motions, event magnitudes were computed as averages across all contributing stations. Borehole sites (e.g., KIT1, BUFF, DOG, and NIKA) typically exhibited lower noise levels and more stable amplitude readings, thereby contributing to more reliable magnitude estimates. Future improvements will incorporate station-specific corrections based on empirical site response functions or $V_{S30}$ values, pending further site characterization.

Two illustrative examples (Figure 3 and Figure 4) highlight how site conditions affect amplitude and magnitude variability. We present recordings of two local microearthquakes, observed at six stations, which occurred on 13 February 2023 (ML = 0.8, Lat: 42.762, Lon: 41.987) and 13 March 2023 (ML = 1.0, Lat: 42.764, Lon: 42.031). The recordings come from the following stations: KIT1, a deep borehole station (250 m, Trillium Compact 20s); GULB, a surface broadband station on a concrete pedestal (MBB-2); and KETI, a surface broadband station on soft soil (MBB-2). Waveforms were bandpass-filtered between 5 and 10 Hz.

For the first event, KIT1 exhibited the lowest ambient noise and the most stable amplitudes, making it ideal for reliable ML estimation. In contrast, KETI recorded significantly higher amplitudes, likely due to site amplification and resonance effects associated with its soft-soil installation, despite being at a slightly greater epicentral distance. GULB showed intermediate waveform characteristics. The resulting magnitude uncertainty for this event, calculated as the standard deviation across stations, was 0.28.

For the second event, signal amplitudes at KETI were not markedly amplified, resulting in reduced inter-station variability. The associated magnitude uncertainty in this case was 0.23.

These examples illustrate the variability of amplitude-based magnitude estimates across stations with differing site conditions. By applying consistent preprocessing, multi-station averaging, and maintaining a well-distributed network geometry, this study ensures that such variability does not significantly affect the statistical integrity of the catalog. Given that the primary objectives focus on temporal and structural seismicity analyses rather than absolute source scaling, the adopted approach provides sufficient accuracy and robustness. Regional calibration of the ML scale remains a target for future work.

4. Methods

4.1. The Frequency–Magnitude Distribution

The Gutenberg–Richter (GR) law [32] describes the relationship between the number of earthquakes, N, with a magnitude greater than or equal to a threshold magnitude $M_{th}$ and the threshold. It is given by the equation

(1)${log}_{10}N=a-b{M}_{th}.$

The coefficient b plays a crucial role as it reflects the proportion between smaller and larger events, and it also informs about the stress state within the Earth’s crust. Estimating precisely the b-value is essential for an accurate assessment of seismic hazard [33,34].

The b-value can be estimated by using the least square linear regression (LSR) method on the logarithm of the cumulative number of earthquakes with $M_{\mathrm{th}}\ge M_c$ [35], where $M_c$ is the completeness magnitude of the earthquake catalog corresponding to the smallest magnitude that can be detected by the seismic network [36], as follows:

(2)$b={\displaystyle \frac{N{\sum }_{i=1}^N{M}_{\mathrm{th},i}{log}_{10}\left(N_{M\ge M_{\mathrm{th},i}}\right)-{\sum }_{i=1}^N{M}_{\mathrm{th},i}{\sum }_{i=1}^N{log}_{10}\left(N_{M\ge M_{\mathrm{th},i}}\right)}{N{\sum }_{i=1}^N{M}_{\mathrm{th},i}^2-{\left({\sum }_{i=1}^N{M}_{\mathrm{th},i}\right)}^2}}$

Its error is given by the following formula:

(3)${\sigma }_b={\displaystyle \frac{1}{\sqrt{N-2}}}\sqrt{{\displaystyle \frac{{\sum }_{i=1}^N{\left(lo{g}_{10}{N}_{M\ge M_{\mathrm{th},i}}-lo{g}_{10}{\widehat{N}}_{M\ge M_{\mathrm{th},i}}\right)}^2}{{\sum }_{i=1}^N{\left(M_{\mathrm{th},i}-<M>\right)}^2}}}$

where $lo{g}_{10}{\widehat{N}}_{M\ge M_{\mathrm{th},i}}$ represents the linear regression estimate of $lo{g}_{10}{N}_{M\ge M_{\mathrm{th},i}}$.

$M_c$ plays an important role in the estimation of b. $M_c$ can be calculated by using several methods. In particular, by the method of maximum curvature (MAXC) [37], $M_c$ is the magnitude corresponding to the largest frequency of seismic events in the non-cumulative binned frequency–magnitude distribution. In the goodness-of-fit test (GFT) method [37], $M_c$ is estimated by comparing the frequency–magnitude distribution of the observed seismicity with those of synthetic earthquake sequences. The synthetic distributions are generated based on the a and b values estimated from the observed dataset for magnitudes $M\ge M_{co}$, evaluated as the cutoff magnitude $M_{co}$ increases. Let R represent the absolute difference in earthquake counts per magnitude bin between the observed and synthetic distributions. The magnitude of completeness, $M_c$, is defined as the smallest $M_{co}$ for which the observed data for $M\ge M_{co}$ align with a straight line on a log-linear plot, within a specified confidence threshold given by $R=R_c$.

4.2. The Allan Factor

A set of earthquakes can be modeled by a discrete-time point process marked by the magnitude

(4)$x\left(t\right)=\sum _i^N{M}_i\delta (t-t_i)$

where $t_i$ is the time of occurrence of the event i with magnitude $M_i$.

By partitioning the time axis into non-overlapping intervals of duration $\tau$, referred to as the timescale, we construct the sequence $N_k\left(\tau \right)$, which is the series the number of earthquakes occurring within the k-th time window [38].

The AF is defined by the following formula:

(5)$AF\left(\tau \right)={\displaystyle \frac{<{(N_{k+1}\left(\tau \right)-N_k\left(\tau \right))}^2>}{2<N_k\left(\tau \right)>}}.$

The AF was employed in the investigation the time dynamics of several phenomena [39,40,41,42]. In a Poissonian seismic sequence, where events are independent and uncorrelated, the AF remains relatively flat across all considered timescales, oscillating around a value of 1, except at very large timescales where finite-size effects become significant [43,44]. In contrast, a clustered earthquake sequence exhibits an AF that varies with the timescale $\tau$, typically increasing following a power-law trend, as expected in a fractal (self-similar) process.

(6)$AF\left(\tau \right)\approx {\left({\displaystyle \frac{\tau }{{\tau }_0}}\right)}^{\alpha }.$

In this case, the exponent of the power-law $\alpha$, known as fractal exponent, measures the strength of time clusterization; ${\tau }_0$ is called fractal onset time and marks the smaller limit of significant scaling behavior in the AF [38]. Thus, if $\alpha \approx 0$, the earthquake sequence can be considered Poissonian, while if $\alpha >0$, it is clusterized.

4.3. The Global and Local Coefficient of Variation

If an earthquake series is represented by the sequence of the interoccurrence intervals, the following quantities can be defined: the global coefficient of variation $C_v$ [45]

(7)$C_v={\displaystyle \frac{{\sigma }_T}{{\mu }_T}}$

and the local coefficient of variation $L_v$ [46]

(8)$L_v={\displaystyle \frac{1}{N-1}}\sum _{i=1}^{N-1}3{\displaystyle \frac{{(T_i-T_{i+1})}^2}{{(T_i+T_{i+1})}^2}}$

where $\{T_i,i=1,...,N\}$ is the series of the interoccurrence times and ${\sigma }_T$ and ${\mu }_T$ are their standard deviation and mean.

These two metrics are both employed to assess the degree of time clustering within an earthquake sequence. In the case of a Poissonian process, where events are randomly distributed, both $C_v$ and $L_v$ take on values close to 1. For regular or periodic sequences, these values drop below 1, indicating reduced variability, whereas for clustered sequences, they exceed 1. The coefficient of variation, $C_v$, reflects the overall clustering behavior of the sequence and can be influenced by fluctuations in the event rate. On the other hand, the local variation, $L_v$, captures the variability between successive interevent times and is less sensitive to gradual changes in the average rate. For example, if two periodic sequences are concatenated, the combined sequence may appear strongly clustered on a global scale—resulting in $C_v>1$—but remains locally periodic, leading to $L_v=0$.

4.4. The Detrended Fluctuation Analysis

Detrended fluctuation analysis (DFA) was introduced as a method to identify long-range correlations in non-stationary signals [47], and it has since been extensively employed [48,49,50,51,52]. Telesca et al. [53] used DFA to analyze earthquake magnitudes and discovered that the scaling exponent of the fluctuation function augmented during the reactivation phases of volcanic activity at El Hierro (Spain) from 2011 to 2014. Varotsos et al. [54] considered DFA results for earthquake magnitudes as a tool for forecasting earthquakes. Additionally, Lennartz et al. [55] studied the long-range correlations in the earthquake magnitudes in Northern and Southern California using DFA. Varotsos et al. [56] demonstrated that the DFA scaling exponent in California’s seismicity could indicate disruptions in long-range correlations preceding major seismic events.

The scaling behavior of seismic interevent time series has been examined in various studies using DFA [57]. Investigating long-range correlations in both earthquake magnitudes and interevent times is therefore a key approach for understanding the underlying dynamics of seismic activity.

The DFA method is described below:

The series $x_i$, where $i=1,\dots ,N$, and N is the size of the series, is integrated and the integrated series $y_k$ is divided into boxes of equal size n.

For each n-size box, a line $y_{n,k}$ fits $y_k$ by the least square method and is subtracted from $y_k$.

The fluctuation, $F_n$, is calculated as

(9)$F_n=\sqrt{{\displaystyle \frac{1}{N}}\sum _{k=1}^N{\left(y_k-y_{n,k}\right)}^2}$

Steps (i)–(iii) are repeated for all the available box sizes n. A power-law relationship between $F_n$ and n indicates the presence of long-range correlation in the series, as follows:

(10)$F_n\propto n^{\alpha }$

The numerical value of the scaling exponent $\alpha$ provides insight into the temporal correlation structure of the seismic series, as follows:

If $\alpha =0.5$, the series is uncorrelated.

If $\alpha >0.5$, the series exhibits persistent correlations, meaning that an increase (or decrease) in one period is likely to be followed by a similar trend in the next.

If $\alpha <0.5$, the series displays antipersistent correlations, where an increase is typically followed by a decrease, and vice versa.

4.5. The Natural Visibility Graph

The Natural Visibility Graph (NVG) transforms a time series into a graph or network. Two arbitrary data points, $(t_a,y_a)$ and $(t_b,y_b)$, are considered connected nodes in the corresponding graph if any other data point $(t_c,y_c)$ located between them satisfies the condition given by [18]

(11)$y\left(t_c\right)<y\left(t_b\right)+(y\left(t_a\right)-y\left(t_b\right)){\displaystyle \frac{t_b-t_c}{t_b-t_a}}.$

The NVG method has been utilized across various research areas, including economics [58,59,60], weather forecasting [61,62], medicine [63,64,65], earthquakes [66], biology [67], among others. One of the key parameters in the NVG method is the degree, which indicates the number of connections originating from a node. Typically, higher values in the series function as hubs within the graph, as they tend to “attract” more connections compared with nodes associated with lower values.

An earthquake sequence can be examined using the NVG method by treating the magnitudes as nodes of a graph and establishing links between magnitudes that satisfy the “visibility” criterion (Equation (11)). In this way, a seismic sequence can be transformed into a graph using the NVG method.

The adjacency matrix $A=\left\{a_{ij}\right\}$ describes the connections between nodes as follows:

(12)$a_{ij}=\left\{\begin{array}{c}1\mathrm{if} \mathrm{nodes} \mathrm{i} \mathrm{and} \mathrm{j} \mathrm{are} \mathrm{connected}\hfill \\ 0\mathrm{otherwise}\hfill \end{array}\right..$

Considering two time series with adjacency matrices $A^1=\left\{a_{ij}^1\right\}$ and $A^2=\left\{a_{ij}^2\right\}$, the $Average$$Edge$$Overlap$ (AEO) [68,69] is given by

(13)${\omega }_{1,2}={\displaystyle \frac{{\sum }_i{\sum }_{j>i}{\sum }_{\nu =1}^2{a}_{ij}^{\nu }}{M{\sum }_i{\sum }_{j>i}\left(1-{\delta }_{0,{\sum }_{\nu =1}^2{a}_{ij}^{\nu }}\right)}}.$

where ${\delta }_{0,{\sum }_{\nu =1}^2{a}_{ij}^{\nu }}$ is the Kronecker delta. ${\omega }_{1,2}$ can vary between 1 (when both the adiacency matrices are identical) and ${\displaystyle \frac{1}{2}}$ (when each edge exists just for one VG); we can interpret these two conditions like the maximum of correlation and the absence of correlation among the variables, respectively.

The AEO defined in Equation (13) can be used to estimate the time-lag between two series. Given two time series $\left\{x_t\right\}$ and $\left\{y_t\right\}$, keeping the first fixed, we construct the time-shifted copies of the second one, $\left\{y_{t-{\tau }_1}\right\},\left\{y_{t-{\tau }_2}\right\},\left\{y_{t-{\tau }_3}\right\},...$. For $\left\{x_t\right\}$ and each of the time-shifted copies of $\left\{y_t\right\}$, we calculate the adjacency matrix, obtained by applying Equation (11). Indicating as $A^x$ the adiacency matrix of $x_t$ and as $A_{{\tau }_1}^y$, $A_{{\tau }_2}^y$, $A_{{\tau }_3}^y$, etc., the adjacency matrices of the delayed copies of $y_t$, from the maximum of sequence $\left\{{\omega }_{x,y}\left({\tau }_1\right)\right),{\omega }_{x,y}\left({\tau }_2\right)),{\omega }_{x,y}\left({\tau }_3\right))...\}$, the time-lag between the series $\left\{x_t\right\}$ and $\left\{y_t\right\}$ can be derived. As an example, consider two sinusoids with delay $d=5$, $y\left(t\right)=100sin\left({\displaystyle \frac{2\pi t}{20}}\right)+w(0,5)$ and $y_d\left(t\right)=100sin\left({\displaystyle \frac{2\pi (t-5)}{20}}\right)+w(0,40)$, where $w(0,5)$ and $w(0,40)$ represent white noise with ${\mu }_{y,y_d}=0$ and variances ${\sigma }_y^2=5$ and ${\sigma }_{y_d}^2=40$, respectively (Figure 5a). Plotting ${\omega }_{y,y_d}$ between the two sinusoids, the maximum is at $\tau =5$, which is exactly the delay of the second sinusoid with respect to the first one (Figure 5b).

4.6. The Correlogram-Based Periodogram

For short time series, an effective algorithm for computing the power spectrum is the correlogram-based periodogram.

Consider the simple model of periodic series, as follows:

(14)$y_n=\beta cos(\omega t+\varphi )+{ϵ}_n$

with $\beta$ > 0, $0<\omega <\pi$, $\varphi$ distributed as a uniform variable in $(-\pi ,\pi ]$, and ${ϵ}_n$ being a series of uncorrelated random variables with mean 0 and standard deviation $\sigma$, independent of $\varphi$.

The classical periodogram is given by

(15)$I\left(\omega \right)={\displaystyle \frac{1}{N}}|\sum _{n=1}^N{y}_n{e}^{-i\omega n}{|}^2,0\le \omega \le \pi$

where N is the length of the series, and it is calculated for

(16)${\omega }_l={\displaystyle \frac{2\pi l}{N}},l=0,1,...,a$

where $a=\left[{\displaystyle \frac{N-1}{2}}\right]$ and $\left[x\right]$ indicates the integer part of x. If the series is modulated by a strong cycle with frequency ${\omega }_0$, the periodogram will display a peak at ${\omega }_0$. In contrast, for an uncorrelated series ($\beta =0$), the periodogram remains nearly flat, with no distinct peak [70]. Ahdesmäki et al. [71] introduced the correlogram-based periodogram, which is more effective than the classical periodogram in identifying the main cycle of a series.

The classical periodogram $I\left(\omega \right)$ is equivalent to the correlogram

(17)$S\left(\omega \right)=\sum _{k=-N+1}^{N-1}\widehat{r}\left(k\right)e^{-i\omega k}$

where

(18)$\widehat{r}\left(m\right)={\displaystyle \frac{1}{N}}\sum _{k=1}^{N-m}{y}_k{y}_{k+m}$

is the estimator of the autocorrelation function. Using the sample autocorrelation function

(19)$\rho \left(m\right)={\displaystyle \frac{\frac{1}{N}{\sum }_{i=1}^N(x_i-\overline{x})(y_i-\overline{y})}{{\sigma }_x{\sigma }_y}}$

where $\overline{x}$ indicates the mean, Ahdesmäki et al. [71] proposed the following correlogram-based periodogram to estimate the power spectrum of a series:

(20)$\tilde{S}\left(\omega \right)=2Re\left(\sum _{k=0}^L\tilde{\rho }\left(k\right)e^{-i\omega k}\right)-\tilde{\rho }\left(0\right).$

where $\mathit{Re}\left(x\right)$ denotes the real part of x, N is the length of the series, and $\tilde{\rho }\left(m\right)$ is the biased Spearman’s rank correlation coefficient. In contrast with the standard periodogram, the correlogram-based periodogram can take on negative values, so the absolute value is typically considered.

To test a frequency ${\omega }_l$, the g-statistics can be used, as follows:

(21)$g={\displaystyle \frac{\left|\tilde{S}\left({\omega }_l\right)\right|}{{\sum }_{i=1}^a\left|\tilde{S}\left({\omega }_i\right)\right|}}.$

If g is large, the tested frequency is considered highly significant. To assess the significance of the tested frequency, a simulation-based method can be applied [71]: a series of random sequences, modeled as in Equation (14), are generated. For each simulated series, the g-value is computed using Equation (21). The distribution of the g-statistics is then estimated from all the calculated g-values, and the p-value is determined. A low p-value indicates the presence of a strong periodicity.

5. Results

A total of 474 seismic events were recorded in the vicinity of the Enguri Dam, with magnitudes ranging from −1.3 to 2.8. Figure 6 shows the time distribution of the analyzed seismicity. The frequency–magnitude distribution of the investigated seismicity is shown in Figure 7. The completeness magnitude $M_c$ was found to be 0.4 for both the MAXC and GFT (with $R_c$ = 90%) methods. The b-value was calculated as $1.02\pm 0.019$ using Equation (2). Although the value of b suggests that the Enguri Dam area may be in a relatively stable state, there could exist possible localized stress accumulation possibly near active faults or stress concentration zones influenced by the dam.

We selected only those earthquakes with magnitude $M\ge 0.4$; thus, the number of analyzed earthquakes reduced from 474 to 279 (represented by the red triangles in Figure 6) and occurred at depths of up to 30 km and within a maximum distance of 33 km from the center of the dam. Various statistical analyses were applied to this dataset with the objective of achieving the most comprehensive characterization possible of the temporal dynamics of the seismic activity around the Enguri Dam.

Figure 8 presents the sequence of interoccurrence times (in days). The global and local coefficients of variation are $C_v=1.283$ and $L_v=1.404$, respectively. The 95% confidence intervals for these estimates, [0.888, 1.138] and [0.875, 1.116], were determined based on 100 simulated Poissonian seismic sequences, which match the original dataset in size and mean interoccurrence time. These intervals were derived as the 2.5% and 97.5% percentiles of the $C_v$ and $L_v$ distributions calculated from the Poissonian sequences. Since the $C_v$ and $L_v$ values for the Enguri seismicity exceed their respective 95% confidence intervals, the analyzed earthquakes exhibit a time-clustered process on both global and local scales, meaning that earthquakes tend to occur in clusters rather than randomly over time. The local coefficient of variation is even larger than the global one, which could suggest stress accumulation and release cycles in certain fault segments, as well as hydro-seismic effects due to water level changes.

Figure 9 displays the Allan factor (AF) for the analyzed seismic activity alongside the 95% confidence bands corresponding to two types of surrogate sequences: 100 Poissonian sequences sharing the same rate as the original dataset (in green) and 100 randomized sequences produced by shuffling the original interoccurrence times, preserving the same probability density function of interoccurrence times, $p\left(\tau \right)$ (in red). The AF slightly increases with the timescale $\tau$, indicating that the seismicity exhibits a certain degree of time clusterization. However, from the numerical value of the scaling exponent $\alpha$, which quantifies the strength of the clusterization, estimated over the range of the large timescales, such clustering appears rather weak. Except the range of very large timescales, where the AF is characterized by a certain fluctuating behavior, it significantly exceeds the Poissonian 95% confidence band, indicating non-Poissonian behavior across the considered timescales. However, the AF remains within the 95% confidence interval derived from the 100 random shuffles, suggesting that the clustering behavior could be very likely due to the behavior of the $p\left(\tau \right)$. The AF results might suggest that microearthquakes near the high dam could exhibit a noticeable time clustering that is weak, meaning that there is no strong external trigger. However, the results leave open the possibility that subtle processes like pore pressure changes or small stress redistributions due to reservoir loading/unloading may play a role in shaping the seismicity pattern.

Figure 10 illustrates the fluctuation function for both the interoccurrence times (Figure 10a) and the magnitudes (Figure 10b). The interoccurrence times exhibit a mildly persistent behavior, suggesting that an increase or decrease in the interoccurrence time during one period is likely to be followed by a similar trend in the subsequent period. This suggests that past seismic activity slightly influences future events, potentially due to periodic stress accumulation and release near the dam. Fluctuations in reservoir water levels can alter the stress state of surrounding rock, gradually increasing pressure until microearthquakes are triggered. This process appears not entirely random but exhibits some memory of prior seismic events—implying that once stress is accumulated, it may require time to be released again, resulting in clustered seismic activity. While the timing of earthquakes shows such memory effects, the magnitudes do not exhibit any clear temporal pattern. This supports the notion that earthquake sizes follow statistical distributions rather than being governed by a predictable process.

We analyzed the seismic dataset by means of the NVG method. Figure 11a shows the links among the magnitudes drawn on the base of the “visibility” criterion of Equation (11).

Following Telesca et al. [72], who identified an intriguing relationship between the degree k in a seismic sequence and earthquake magnitude M, the slope of the linear fit describing this relationship is referred to as the k–M slope. Several studies using observational [73], experimental [20], and physics-based synthetic [74] seismic catalogs have reported a positive correlation between the k–M slope and the Gutenberg–Richter b-value. More recently, Li et al. [75] demonstrated that this correlation is both universal and stable.

Figure 11b illustrates the relationship between the degree k and magnitude M in our earthquake dataset, along with the corresponding linear regression line. The slope of this fit, known as the k–M slope, is positive ($8.78$), indicating a direct correlation between degree and magnitude. In other words, larger-magnitude earthquakes tend to have higher degrees—they function as “hubs” in the seismic network, drawing more connections than smaller events, and are thus “visible” to more earthquakes.

Moreover, Figure 12 shows that the Enguri seismic dataset aligns well with the general trend observed in other seismic datasets, further supporting the robustness of the k–M slope–b-value relationship.

Figure 13 illustrates the monthly changes in the number of earthquakes (with $M\ge M_c$) alongside the monthly average water level of the Enguri reservoir. The monthly variation in earthquake frequency remains below 20 events per month. Meanwhile, the average water level is characterized by a distinct annual cycle. This cyclical behavior is further supported by correlogram-based periodogram analysis. Figure 14 displays the correlogram-based periodogram for the monthly mean water level (in yellow) and the monthly earthquake counts (in blue); the 99% significant periods are indicated by circles. The water level exhibits a clear and significant quasi-annual periodicity, which is also mirrored in the monthly earthquake counts, though with less intensity. However, given the availability of the water level time series starting from 1978, we applied the correlogram-based method to the entire dataset in order to accurately determine the periodicity of the monthly water level. The resulting spectrum, shown in Figure 15, indicates that the dominant periodicity is 12 months. A secondary peak at 6 months is present and represents a harmonic of the annual cycle.

We applied the NVG to both monthly time series and computed the AEO while varying the time lag $\tau$ from −12 to 12 months.

To evaluate the significance of the obtained $\omega$ values, we constructed a null distribution representing the case of no correlation between the two graphs. This was implemented using a permutation-based approach, structured as follows:

The adjacency matrices of both graphs were randomly permuted, while preserving the null elements on the diagonal and maintaining symmetry. These constraints ensure that node connections are shuffled randomly, but the overall graph structure remains consistent with the visibility graph rules.

For each permutation and time lag $\tau$, the $\omega$ value was computed.

This process was repeated 100 times, resulting in a null distribution of $\omega$ values for each time lag $\tau$.

The 95th percentile of each distribution was calculated to determine the threshold for statistical significance.

These percentiles were then used to construct a 95% confidence envelope for $\omega$. If the observed $\omega$ value exceeds this envelope at a given $\tau$, it is considered statistically significant.

The results shown in Figure 16 indicate a significant correlation between the two time series, with $\omega$ consistently remaining above the 95% confidence curve for all values of $\tau$. Notably, the highest value of $\omega$ is observed at $\tau =7$ months, indicating that the monthly mean water level leads the monthly earthquake counts by approximately 7 months (±0.5 months, considering the monthly resolution of the time series). This result further supports the possibility that seismic activity around the Enguri reservoir may be influenced by the loading and unloading cycles associated with reservoir operations, which, although not dominant, could play a contributing role in triggering observed seismic events. Furthermore, the found delay would strengthen the hypothesis that the triggering mechanism by which the dam influences fault activation—already critically stressed by the tectonic stress field—is not the immediate undrained elastic response to loading, but rather a delayed seismic response driven by the diffusion of pore fluid pressures.

The effect of the annual loading/unloading of the dam can be approximated by a periodic pore-pressure variation at the base of the dam. This can be expressed as

(22)$\omega ={\displaystyle \frac{2\pi }{T}}$

where T is the period of the loading/unloading cycle. Following Saar and Manga [79], the hydraulic diffusivity D can be expressed as

(23)$D={\displaystyle \frac{h^2}{2\omega {\tau }^2}}$

where h is the mean depth of the seismic events, $\omega$ is given by Equation (22), and $\tau$ is the time delay between water level and seismicity (which in our case are the monthly mean water level and earthquake counts, respectively).

In our case, $T=11$ months (as shown in Figure 15), $\tau =7\mathrm{months}$ (as shown in Figure 16), and $h=8.59$. Substituting these values into Equation (23) gives

$D=0.56{\mathrm{m}}^2/\mathrm{s}$

This estimated value of D falls within the typical range of seismic hydraulic diffusivity (0.1–100 m²/s) reported by Talwani and Acree [80] in their study of seismicity near the Monticello Reservoir (Italy).

6. Conclusions

In this study, we conducted several statistical analyses on microearthquake seismicity in the Enguri area (Georgia) that occurred between 2020 and 2023. These events were recorded by a newly installed seismic network developed as part of the DAMAST project. The enhanced sensitivity of this network enabled the detection of even very small seismic events, allowing a more detailed picture of the temporal dynamics of seismicity in the area. Our main findings can be summarized as follows: (1) the seismicity around the Enguri Dam is likely influenced by a combination of natural tectonic processes and subtle reservoir-induced stress changes; (2) although the dam does not appear to be causing strong seismic forcing, it could have a role in triggering the seismicity that appears to be delayed by about 7 months from the water level; (3) localized stress variations and time clustering suggest that water level changes may modulate seismic activity over time; and (4) the mild persistence of interoccurrence times aligns with the idea of stress buildup and delayed earthquake triggering due to reservoir loading effect.

Author Contributions

Conceptualization, L.T.; methodology, L.T.; software, L.T.; formal analysis, L.T.; investigation, L.T., T.C. and N.T. (Nino Tsereteli); data curation, N.T. (Nino Tsereteli) and N.T. (Nazi Tugushi); writing—original draft preparation, L.T. and N.T. (Nino Tsereteli); writing—review and editing, L.T., T.C. and N.T. (Nino Tsereteli); visualization, L.T., N.T. (Nino Tsereteli) and N.T. (Nazi Tugushi). All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

The data presented in this study are available on request.

Acknowledgments

We acknowledge LTD Engurhesi for providing the data on the water level of the Enguri Dam and the DAMAST project.

Conflicts of Interest

The authors declare no conflicts of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1 Seismicity and active fault structures in the Enguri Dam region. Yellow circles indicate microearthquakes recorded between 2020 and 2023; large blue circles mark historical earthquakes. DAMAST seismic stations are shown as black triangles. Active faults from Basili et al. [17] (hatched areas) and locally identified Ingrishi fault segments (red dashed lines) are overlaid. Major rivers, the Jvari Reservoir (JR), Gali Reservoir (GR), and the location of the Enguri Dam (white triangle) are labeled. Inset: Regional map showing the location of the study area (red box) within Georgia and surrounding countries.

Figure 2 Median power spectral density (PPSD) for vertical component across all stations compared with Peterson noise models.

Figure 3 Event recordings on three seismic stations: (a) KITI, (b) GULB, (c) KETI. Records of the 13 February 2023 event (ML = 0.8; Lat: 42.762, Lon: 41.987, Depth = 8.7 km).

Figure 4 Event recordings on three seismic stations: (a) KITI, (b) KETI, (c) GULB. Records of the 13 March 2023 event (ML = 1.0; Lat: 42.764, Lon: 42.031, Depth = 3.1 km).

Figure 5 (a) Comparison between two sinusoids with delay $d=5$, $y\left(t\right)=100sin\left({\displaystyle \frac{2\pi t}{20}}\right)+w(0,5)$ and $y_d\left(t\right)=100sin\left({\displaystyle \frac{2\pi (t-5)}{20}}\right)+w(0,40)$, where $w(0,5)$ and $w(0,40)$ represent white noise with ${\mu }_{y,y_d}=0$ and variances ${\sigma }_y^2=5$ and ${\sigma }_{y_d}^2=40$, respectively. (b) Variation of ${\omega }_{y,y_d}$ with the time lag $\tau$; the maximum of ${\omega }_{y,y_d}$ indicates the time lag of $d=5$ between the two sinusoids.

Figure 6 Time distribution of the earthquakes in the investigated area. The red triangles represent the earthquakes above the completeness magnitude.

Figure 7 Cumulative (red) and non-cumulative (blue) frequency–magnitude distribution of earthquakes occurring at depths of up to 30 km and within a 33 km radius from the center of the dam. The slope of the regression line (black) is the estimate of the b-value.

Figure 8 Sequence of the interoccurrence times.

Figure 9 The Allan factor of the earthquake sequence is shown in blue. The red dotted lines indicate the 95% confidence interval derived from random shuffles of the original interevent times, while the green dotted lines represent the 95% confidence interval based on Poissonian surrogate sequences. The black solid line shows the linear regression of the Allan factor over the larger timescales, with its slope corresponding to the fractal exponent.

Figure 10 Detrended fluctuation analysis of the (a) interoccurrence times and (b) magnitudes. the red dotted lines show the linear regression of the fluctuation function over the scale n.

Figure 11 (a) Sequence of analyzed earthquakes. The links among the magnitudes are defined by the Natural Visibility Graph. (b) Relationship between the degree k and the magnitude M.

Figure 12 Relationship between the k-M slope and the b-value for the Enguri dataset compared with that of the seismic datasets analyzed in previous studies (Iran [19], Italy and Taiwan [20], Mexico [72], Pannonia [76], R1 and R2 [77], and ST2 [78]).

Figure 13 Time variation of the monthly earthquake counts (blue) and mean water level (red).

Figure 14 Correlogram-based periodogram of the monthly earthquake counts (blue) and the monthly mean water level (yellow). The circles indicates the 99% significant periods in both periodograms.

Figure 15 Correlogram-based periodogram of the monthly mean water level from 1978 to 2023. The dominant periodicity occurs at 12 months, corresponding to the main cycle, while the 6-month peak represents its first harmonic.

Figure 16 Average edge overlap $\omega$ between the monthly time series of earthquake counts and mean water level. The red dotted line represents the 95% confidence curve.

Table 1

Coordinates and elevations of seismic stations.

Station	Lat	Lon	Elevation (m)	Installation Depth (m)	Starting Time	Ending Time
BRID	42.7424	41.9827	496	Surface	21/10/2020	currently
BUFF	42.7052	41.9330	369	19	21/10/2020	currently
NIKA	42.6583	41.9701	277	17	21/10/2020	currently
KIT1	42.7563	42.0336	427	250	23/10/2020	currently
DOG	42.7901	42.0518	740	20	15/12/2020	currently
GULB	42.7734	41.9329	1877	Surface	22/09/2021	05/06/23
KETI	42.7185	42.0841	418	Surface	06/10/21	currently
HPP	42.6672	41.8504	300	Surface	01/01/22	currently
GAL	42.6316	41.7356	54	Surface	23/10/2023	currently
OKM	42.7262	41.7703	227	Surface	31/10/2023	currently

Table 2

Instrument response characteristics by sensor type and station class.

Sensor Type	Digitizer	Stations	Response Type	Bandwidth (Flat Velocity Response)	Dynamic Range
Trillium Compact Posthole 20s	Centaur-4 [CTR4-3s] (3-CH)	KIT1, BUFF, NIKA, DOG	Broadband (velocity)	0.008–100+ Hz	>156 dB @ 1 Hz
Kinemetrics MBB-2	Centaur-4 [CTR4-3s] (3-CH)	KETI, BRID, GULB	Broadband (velocity)	0.03–80 Hz	>155 dB @ 1 Hz
4.5 Hz Geophone PE-6/B	DATA-CUBE3 Type 2	HPP, OKM, GAL	Short-period (velocity)	4.5–80 Hz	100–110 dB

Table 3

Root-mean-square (Vrms) velocity amplitudes across vertical and horizontal components in two frequency bands (1–10 Hz and 5–15 Hz) for all stations in the Enguri seismic network.

Station	Vrms Vertical ($\mathsf{\mu }$m/s)			Vrms HHX ($\mathsf{\mu }$m/s)			Vrms HHY ($\mathsf{\mu }$m/s)			Vrms Horizontal Avg ($\mathsf{\mu }$m/s)
	1–10 Hz	5–15 Hz		1–10 Hz	5–15 Hz		1–10 Hz	5–15 Hz		1–10 Hz	5–15 Hz
BRID	0.0230	0.0107		0.0324	0.0172		0.0383	0.0195		0.0355	0.0184
BUFF	0.0136	0.0082		0.0201	0.0135		0.0150	0.0081		0.0177	0.0112
DOG	0.0168	0.0031		0.0189	0.0040		0.0178	0.0040		0.0183	0.0040
KETI	0.0637	0.0306		0.1315	0.0392		0.1333	0.0357		0.1324	0.0375
KIT1	0.0092	0.0018		0.0073	0.0016		0.0082	0.0015		0.0078	0.0015
GULB	0.0092	0.0022		0.0080	0.0020		0.0112	0.0040		0.0096	0.0030
NIKA	0.0247	0.0050		0.0240	0.0031		0.0290	0.0036		0.0266	0.0033
GAL	0.0355	0.0400		0.1509	0.1735		0.0967	0.1309		0.1267	0.1537
HPP	0.0591	0.0585		0.0723	0.0561		0.0693	0.0552		0.0708	0.0556
OKM	0.0250	0.0163		0.0248	0.0144		0.0253	0.0138		0.0251	0.0141

References

1. Reilinger, R.; McClusky, S.; Vernant, P.; Lawrence, S.; Ergintav, S.; Cakmak, R.; Ozener, H.; Kadirov, F.; Guliev, I.; Stepanyan, R. . GPS constraints on continental deformation in the Africa-Arabia-Eurasia continental collision zone and implications for the dynamics of plate interactions. J. Geophys. Res. Solid Earth; 2006; 111, pp. 1-26. [DOI: https://dx.doi.org/10.1029/2005JB004051]

2. Allen, M.; Jackson, J.; Walker, R. Late Cenozoic reorganization of the Arabia-Eurasia collision and the comparison of short-term and long-term deformation rates. Tectonics; 2004; 23, pp. 1-16. [DOI: https://dx.doi.org/10.1029/2003TC001530]

3. Mosar, J.; Kangarli, T.; Bochud, M.; Glasmacher, U.A.; Rast, A.; Brunet, M.F.; Sosson, M. Cenozoic-recent tectonics and uplift in the greater Caucasus: A perspective from Azerbaijan. Geol. Soc. Spec. Publ.; 2010; 340, pp. 261-280. [DOI: https://dx.doi.org/10.1144/SP340.12]

4. Zor, E.; Sandvol, E.; Gürbüz, C.; Türkelli, N.; Seber, D.; Barazangi, M. The crustal structure of the East Anatolian Plateau (Turkey) from receiver functions. Geophys. Res. Lett.; 2006; 33, L08301. [DOI: https://dx.doi.org/10.1029/2003GL018192]

5. Copley, A.; Jackson, J. Active tectonics of the Turkish–Iranian Plateau. Tectonics; 2006; 25, TC6006. [DOI: https://dx.doi.org/10.1029/2005TC001906]

6. Mcclusky, S.; Balassanian, S.; Barka, A.; Demir, C.; Ergintav, S.; Georgiev, I.; Gurkan, O.; Hamburger, M.; Hurst, K.; Kahle, H. . Global positioning system constraints on plate kinematics and dynamics in the eastern Mediterranean and Caucasus. J. Geophys. Res. Solid Earth; 2000; 105, pp. 5695-5719. [DOI: https://dx.doi.org/10.1029/1999JB900351]

7. Aktuğ, B.; Meherremov, E.; Kurt, M.; Özdemir, S.; Esedov, N.; Lenk, O. GPS constraints on the deformation of Azerbaijan and surrounding regions. J. Geodyn.; 2013; 67, pp. 40-45. [DOI: https://dx.doi.org/10.1016/j.jog.2012.05.007]

8. Kadirov, F.; Gadirov, A.; Abdullayev, N. Gravity modelling of the regional profile across South Caspian basin and tectonic implications. The Modern Problems of Geology and Geophysics of Eastern Caucasus and the South Caspian Depression; Nafta-Press: Baku, Azerbaijan, 2012; Volume 231, 251.

9. Müller, M.; Asch, G.; Ratschbacher, L.; Stöckhert, B.; Zor, E.; Akal, C.; Sipahioglu, S.; Lueth, S.; Kind, R. Present-day deformation across the Lesser and Greater Caucasus from geodetic and seismological data. Earth Planet. Sci. Lett.; 2022; 579, 117369.

10. Sokhadze, G.; Floyd, M.; Godoladze, T.; King, R.; Cowgill, E.; Javakhishvili, Z.; Hahubia, G.; Reilinger, R. Active convergence between the lesser and greater Caucasus in Georgia: Constraints on the tectonic evolution of the lesser-greater Caucasus continental collision. Earth Planet. Sci. Lett.; 2018; 481, pp. 154-161. [DOI: https://dx.doi.org/10.1016/j.epsl.2017.10.007]

11. Caputo, M.; Gamkrelidze, I.; Malvezzi, V.; Sgrigna, V.; Shengelaia, G.; Zilpimiani, D. Geostructural basis and geophysical investigations for the seismic hazard assessment and prediction in the Caucasus. Il Nuovo Cimento C; 2000; 23, pp. 191-216.

12. Philip, H.; Cisternas, A.; Gvishiani, A.; Gorshkov, A. The Caucasus: An actual example of the initial stages of continental collision. Tectonophysics; 1989; 161, pp. 1-21. [DOI: https://dx.doi.org/10.1016/0040-1951(89)90297-7]

13. Tan, O.; Taymaz, T. Active tectonics of the Caucasus: Earthquake source mechanisms and rupture histories obtained from inversion of teleseismic body waveforms. Spec. Pap. Geol. Soc. Am.; 2006; 409, 531.

14. Telesca, L.; Tsereteli, N.; Chelidze, T.; Lapenna, V. Spectral investigation of the relationship between seismicity and water level in the Enguri high dam area (Georgia). Geosciences; 2024; 14, 22. [DOI: https://dx.doi.org/10.3390/geosciences14010022]

15. Karamzadeh, N.; Tsereteli, N.; Gaucher, E.; Tugushi, N.; Shubladze, T.; Varazanashvili, O.; Rietbrock, A. Seismological study around the Enguri dam reservoir (Georgia) based on old catalogs and ongoing monitoring. J. Seismol.; 2023; 27, pp. 953-977. [DOI: https://dx.doi.org/10.1007/s10950-023-10173-5]

16. Braun, T.; Cesca, S.; Kuhn, D.; Martirosian-Janssen, A.; Dahm, T. Anthropogenic seismicity in Italy and its relation to tectonics: State of the art and perspectives. Anthropocene; 2018; 21, pp. 80-94. [DOI: https://dx.doi.org/10.1016/j.ancene.2018.02.001]

17. Basili, R.; Danciu, L.; Beauval, C.; Sesetyan, K.; Vilanova, S.P.; Adamia, S.; Arroucau, P.; Atanackov, J.; Baize, S.; Canora, C. . The European Fault-Source Model 2020 (EFSM20): Geologic input data for the European Seismic Hazard Model 2020. Nat. Hazards Earth Syst. Sci.; 2024; 24, pp. 3945-3976. [DOI: https://dx.doi.org/10.5194/nhess-24-3945-2024]

18. Lacasa, L.; Luque, B.; Ballesteros, F.; Luque, J.; Nuno, J.C. From time series to complex networks: The visibility graph. Proc. Natl. Acad. Sci. USA; 2008; 105, pp. 4972-4975. [DOI: https://dx.doi.org/10.1073/pnas.0709247105]

19. Khoshnevis, N.; Taborda, R.; Azizzadeh-Roodpish, S.; Telesca, L. Analysis of the 2005–2016 earthquake sequence in Northern Iran using the visibility graph method. Pure Appl. Geophys.; 2017; 174, pp. 4003-4019. [DOI: https://dx.doi.org/10.1007/s00024-017-1617-8]

20. Telesca, L.; Chen, C.C.; Lovallo, M. Investigating the relationship between seismological and topological properties of seismicity in Italy and Taiwan. Pure Appl. Geophys.; 2020; 177, pp. 4119-4126. [DOI: https://dx.doi.org/10.1007/s00024-020-02470-8]

21. Murusidze, G. Seismicity and the possible position of strong earthquakes sources in the area of the Inguri HPP. Seismic Impacts on Hydrotechnical and Power Structures; Savarensky, E. Pub, House Nauka: Moscow, Russia, 1980; pp. 73-80.

22. McNamara, D.E.; Buland, R.P. Ambient noise levels in the continental United States. Bull. Seismol. Soc. Am.; 2004; 94, pp. 1517-1527. [DOI: https://dx.doi.org/10.1785/012003001]

23. Krischer, L.; Megies, T.; Barsch, R.; Beyreuther, M.; Lecocq, T.; Caudron, C.; Wassermann, J. ObsPy: A bridge for seismology into the scientific Python ecosystem. Comput. Sci. Discov.; 2015; 8, 014003. [DOI: https://dx.doi.org/10.1088/1749-4699/8/1/014003]

24. Peterson, J. Observations and Modeling of Seismic Background Noise; U.S. Geological Survey Open-File Report 93-322 U.S. Geological Survey: Reston, VA, USA, 1993; [DOI: https://dx.doi.org/10.3133/ofr93322]

25. Ringler, A.T.; Hutt, C.R. Self-Noise Models of Seismic Instruments. Seismol. Res. Lett.; 2010; 81, pp. 972-983. [DOI: https://dx.doi.org/10.1785/gssrl.81.6.972]

26. Anthony, R.E.; Ringler, A.T.; Wilson, D.C.; McNamara, D.E. Ambient Noise Levels in the Continental United States. Seismol. Res. Lett.; 2020; 91, pp. 2199-2210.

27. Casey, A.T.; Ringler, A.T.; Anthony, R.E. Ambient Seismic Noise as a Measure of Station Performance. Seismol. Res. Lett.; 2018; 89, pp. 1453-1460. [DOI: https://dx.doi.org/10.1785/0220180045]

28. Sleeman, R.; van Wettum, A.; Trampert, J. Three-channel correlation analysis: A new technique to measure instrumental noise of digitizers and seismic sensors. Bull. Seismol. Soc. Am.; 2006; 96, pp. 258-271. [DOI: https://dx.doi.org/10.1785/0120050032]

29. Richter, C.F. An Instrumental Earthquake Magnitude Scale. Bull. Seismol. Soc. Am.; 1935; 25, pp. 1-32. [DOI: https://dx.doi.org/10.1785/BSSA0250010001]

30. Hutton, L.K.; Boore, D.M. The ML Scale in Southern California. Bull. Seismol. Soc. Am.; 1987; 77, pp. 2074-2094. [DOI: https://dx.doi.org/10.1785/BSSA0770062074]

31. Havskov, J.; Ottemöller, L. SEISAN: The Earthquake Analysis Software for Windows, Solaris, Linux and MacOSX. Seismol. Res. Lett.; 2010; 81, pp. 532-540.

32. Gutenberg, R.; Richter, C.F. Frequency of earthquakes in California. Bull. Seismol. Soc. Am.; 1944; 34, pp. 185-188. [DOI: https://dx.doi.org/10.1785/BSSA0340040185]

33. Scholz, C.H. The frequency-magnitude relation of microfracturing in rock and its relation to earthquakes. Bull. Seismol. Soc. Am.; 1968; 58, pp. 399-415. [DOI: https://dx.doi.org/10.1785/BSSA0580010399]

34. Wyss, M. Towards a physical understanding of the earthquake frequency distribution. Geophys. J. R. Astron. Soc.; 1973; 31, pp. 341-359. [DOI: https://dx.doi.org/10.1111/j.1365-246X.1973.tb06506.x]

35. Kárník, V.; Klíma, K. Magnitude-frequency distribution in the European-Mediterranean earthquake regions. Tectonophysics; 1993; 220, pp. 309-323. [DOI: https://dx.doi.org/10.1016/0040-1951(93)90238-F]

36. Utsu, T. Representation and analysis of the earthquake size distribution: A historical review and some new approaches. Pageoph; 1999; 155, pp. 509-535. [DOI: https://dx.doi.org/10.1007/s000240050276]

37. Wiemer, S.; Wyss, M. Minimum magnitude of completeness in earthquake catalogs: Examples from Alaska, the Western United States, and Japan. Bull. Seismol. Soc. Am.; 2000; 90, pp. 859-869. [DOI: https://dx.doi.org/10.1785/0119990114]

38. Thurner, S.; Lowen, S.B.; Feurstein, M.C.; Heneghan, C.; Feichtinger, H.G.; Teich, M.C. Analysis, synthesis, and estimation of fractal-rate stochastic point processes. Fractals; 1997; 5, pp. 565-596. [DOI: https://dx.doi.org/10.1142/S0218348X97000462]

39. Fadel, P.; Barman, S.; Phillips, S.; Gebber, G. Fractal fluctuations in human respiration. J. Appl. Physiol.; 2004; 97, pp. 2056-2064. [DOI: https://dx.doi.org/10.1152/japplphysiol.00657.2004]

40. Fadel, P.; Orer, H.; Barman, S.; Vongpatanasin, W.; Victor, R.; Gebber, G. Fractal properties of human muscle sympathetic nerve activity. Am. J. Physiol. Circ. Physiol.; 2004; 286, pp. H1076-H1087. [DOI: https://dx.doi.org/10.1152/ajpheart.00577.2003]

41. Yuan, Y.; Zhuang, X.T.; Liu, Z.Y.; Huang, W.G. Time–clustering behavior of sharp fluctuation sequences in Chinese stock markets. Chaos Solitons Fractals; 2012; 45, pp. 838-845. [DOI: https://dx.doi.org/10.1016/j.chaos.2012.02.020]

42. Anteneodo, C.; Malmgren, R.; Chialvo, D. Poissonian bursts in e–mail correspondence. Eur. Phys. J. B—Condensed Matter Complex Syst.; 2010; 75, pp. 389-394. [DOI: https://dx.doi.org/10.1140/epjb/e2010-00139-9]

43. Telesca, L.; ElShafey Fat ElBary, R.; Amin Mohamed, A.E.E.; ElGabry, M. Analysis of the cross-correlation between seismicity and water level in the Aswan area (Egypt) from 1982 to 2010. Nat. Hazards Earth Syst. Sci.; 2012; 12, pp. 2203-2207. [DOI: https://dx.doi.org/10.5194/nhess-12-2203-2012]

44. Serinaldi, F.; Kilsby, C.G. On the sampling distribution of Allan factor estimator for a homogeneous Poisson process and its use to test inhomogeneities at multiple scales. Phys. A Stat. Mech. Its Appl.; 2013; 392, pp. 1080-1089. [DOI: https://dx.doi.org/10.1016/j.physa.2012.11.015]

45. Kagan, Y.Y.; Jackson, D.D. Long-Term Earthquake Clustering. Geophys. J. Int.; 1991; 104, pp. 117-134. [DOI: https://dx.doi.org/10.1111/j.1365-246X.1991.tb02498.x]

46. Shinomoto, S.; Miura, K.; Koyama, S. A measure of local variation of inter-spike intervals. Biosystems; 2005; 79, pp. 67-72. [DOI: https://dx.doi.org/10.1016/j.biosystems.2004.09.023]

47. Peng, C.K.; Havlin, S.; Stanley, H.E.; Goldberger, A.L. Quantification of scaling exponents and crossover phenomena in nonstationary heartbeat time series. Chaos; 1995; 5, pp. 82-87. [DOI: https://dx.doi.org/10.1063/1.166141]

48. Heneghan, C.; McDarby, G. Establishing the relation between detrended fluctuation analysis and power spectral density analysis for stochastic processes. Phys. Rev. E; 2000; 62, pp. 6103-6110. [DOI: https://dx.doi.org/10.1103/PhysRevE.62.6103]

49. Talkner, P.; Weber, R.O. Power spectrum and detrended fluctuation analysis: Application to daily temperatures. Phys. Rev. E; 2000; 62, pp. 150-160. [DOI: https://dx.doi.org/10.1103/PhysRevE.62.150]

50. Kantelhardt, J.W.; Koscielny-Bunde, E.; Rego, H.H.A.; Havlin, S.; Bunde, A. Detecting Long-range Correlations with Detrended Fluctuation Analysis. Phys. A Stat. Mech. Its Appl.; 2001; 295, pp. 441-454. [DOI: https://dx.doi.org/10.1016/S0378-4371(01)00144-3]

51. Colás, A.; Vigil, L.; Vargas, B.; Cuesta-Frau, D.; Varela, M. Detrended Fluctuation Analysis in the prediction of type 2 diabetes mellitus in patients at risk: Model optimization and comparison with other metrics. PLoS ONE; 2019; 14, e0225817. [DOI: https://dx.doi.org/10.1371/journal.pone.0225817]

52. Zhao, S.; Jiang, Y.; He, W.; Mei, Y.; Xie, X.; Wan, S. Detrended fluctuation analysis based on best-fit polynomial. Front. Environ. Sci.; 2022; 10, 1054689. [DOI: https://dx.doi.org/10.3389/fenvs.2022.1054689]

53. Telesca, L.; Lovallo, M.; Lopez, M.C.; Molist, J.M. Multiparametric statistical investigation of seismicity occurred at El Hierro (Canary Islands) from 2011 to 2014. Tectonophysics; 2016; 672–673, pp. 121-128. [DOI: https://dx.doi.org/10.1016/j.tecto.2016.01.045]

54. Varotsos, P.A.; Sarlis, N.V.; Skordas, E.S. Study of the temporal correlations in the magnitude time series before major earthquakes in Japan. J. Geophys. Res. Space Phys.; 2014; 119, pp. 9192-9206. [DOI: https://dx.doi.org/10.1002/2014JA020580]

55. Lennartz, S.; Livina, V.N.; Bunde, A.; Havlin, S. Long-term memory in earthquakes and the distribution of interoccurrence times. Europhys. Lett.; 2008; 81, 69001. [DOI: https://dx.doi.org/10.1209/0295-5075/81/69001]

56. Varotsos, P.; Sarlis, N.; Skordas, E. Scale-specific order parameter fluctuations of seismicity before mainshocks: Natural time and detrended fluctuation analysis. Europhys. Lett.; 2012; 99, 001. [DOI: https://dx.doi.org/10.1209/0295-5075/99/59001]

57. Telesca, L.; Chen, C.C. Fractal and spectral investigation of the shallow seismicity in Taiwan. J. Asian Earth Sci.; 2019; [DOI: https://dx.doi.org/10.1016/j.jseaes.2018.10.009]

58. Zhang, R.; Ashuri, B.; Shyr, Y.; Deng, Y. Forecasting construction cost index based on visibility graph: A network approach. Physica A; 2018; 493, pp. 239-252. [DOI: https://dx.doi.org/10.1016/j.physa.2017.10.052]

59. Long, Y. Visibility graph network analysis of gold price time series. Physica A; 2013; 392, pp. 3374-3384. [DOI: https://dx.doi.org/10.1016/j.physa.2013.03.063]

60. Wang, N.; Li, D.; Wang, Q. Visibility graph analysis on quarterly macroeconomic series of China based on complex network theory. Physica A; 2012; 391, pp. 6543-6555. [DOI: https://dx.doi.org/10.1016/j.physa.2012.07.054]

61. Jiang, W.; Wei, B.; Zhan, J.; Xie, C.; Zhou, D. A visibility graph power averaging aggregation operator: A methodology based on network analysis. Comput. Ind. Eng.; 2016; 101, pp. 260-268. [DOI: https://dx.doi.org/10.1016/j.cie.2016.09.009]

62. Chen, S.; Hu, Y.; Mahadevan, S.; Deng, Y. A visibility graph averaging aggregation operator. Physica A; 2014; 403, pp. 1-12. [DOI: https://dx.doi.org/10.1016/j.physa.2014.02.015]

63. Yu, M.; Hillebrand, A.; Gouw, A.A.; Stam, C.J. Horizontal visibility graph transfer entropy (HVG-TE): A novel metric to characterize directed connectivity in large-scale brain networks. NeuroImage; 2017; 156, pp. 249-264. [DOI: https://dx.doi.org/10.1016/j.neuroimage.2017.05.047]

64. Bhaduri, S.; Ghosh, D. Electroencephalographic data analysis with visibility graph technique for quantitative assessment of brain dysfunction. Clin. EEG Neurosci.; 2015; 46, pp. 218-223. [DOI: https://dx.doi.org/10.1177/1550059414526186] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/24781371]

65. Ahmadlou, M.; Adeli, H.; Adeli, A. New diagnostic EEG markers of the Alzheimer’s disease using visibility graph. J. Neural Transm.; 2010; 117, pp. 1099-1109. [DOI: https://dx.doi.org/10.1007/s00702-010-0450-3] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/20714909]

66. Aguilar-San Juan, B.; Guzmán-Vargas, L. Earthquake magnitude time series: Scaling behavior of visibility networks. Eur. Phys. J. B; 2013; 86, 454. [DOI: https://dx.doi.org/10.1140/epjb/e2013-40762-2]

67. Zheng, M.; Domanskyi, S.; Piermarocchi, C.; Mias, G.I. Visibility graph based temporal community detection with applications in biological time series. Sci. Rep.; 2021; 11, 5623. [DOI: https://dx.doi.org/10.1038/s41598-021-84838-x]

68. Lacasa, L.; Nicosia, V.; Latora, V. Network structure of multivariate time series. Sci. Rep.; 2015; 5, 15508. [DOI: https://dx.doi.org/10.1038/srep15508]

69. Carmona-Cabezas, R.; Gómez-Gómez, J.; Ariza-Villaverde, A.B.; Gutiérrez de Ravé, E.; Jiménez-Hornero, F.J. Multiplex Visibility Graphs as a complementary tool for describing the relation between ground level O₃ and NO₂. Atmos. Pollut. Res.; 2020; 11, pp. 205-212. [DOI: https://dx.doi.org/10.1016/j.apr.2019.10.011]

70. Priestley, M.B. Spectral Analysis and Time Series; Academic Press: London, UK, New York, NY, USA, 1981; 890.

71. Ahdesmäki, M.; Lähdesmäki, H.; Pearson, R.; Huttunen, H.; Yli-Harja, O. Robust detection of periodic time series measured from biological systems. BMC Bioinform.; 2005; 6, 117. [DOI: https://dx.doi.org/10.1186/1471-2105-6-117]

72. Telesca, L.; Lovallo, M.; Ramirez-Rojas, A.; Flores-Marquez, L. Investigating the time dynamics of seismicity by using the visibility graph approach: Application to seismicity of Mexican subduction zone. Physica A; 2013; 391, pp. 1553-1562. [DOI: https://dx.doi.org/10.1016/j.physa.2011.08.042]

73. Azizzadeh-Roodpish, S.; Cramer, C.H. Visibility graph analysis of Alaska crustal and Aleutian subduction zone seismicity: An investigation of the correlation between b value and k–M slope. Pure Appl. Geophys.; 2018; 175, pp. 4241-4252. [DOI: https://dx.doi.org/10.1007/s00024-018-1947-1]

74. Perez-Oregon, J.; Lovallo, M.; Telesca, L. Visibility graph analysis of synthetic earthquakes generated by the Olami–Feder–Christensen spring-block model. Chaos; 2020; 30, 093111. [DOI: https://dx.doi.org/10.1063/5.0007480]

75. Li, L.; Luo, G.; Liu, M. The K-M Slope: A Potential Supplement for b-Value. Seismol. Res. Lett.; 2023; 94, pp. 1892-1899. [DOI: https://dx.doi.org/10.1785/0220220268]

76. Telesca, L.; Lovallo, M.; Toth, L. Visibility graph analysis of 2002–2011 Pannonian seismicity. Physica A; 2014; 416, pp. 219-224. [DOI: https://dx.doi.org/10.1016/j.physa.2014.08.048]

77. Telesca, L.; Lovallo, M.; Ramirez-Rojas, A.; Flores-Marquez, L. Relationship between the frequency magnitude distribution and the visibility graph in the synthetic seismicity generated by a simple stick-slip system with asperities. PLoS ONE; 2014; 9, e106233. [DOI: https://dx.doi.org/10.1371/journal.pone.0106233] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/25162728]

78. Telesca, L.; Thai, A.T.; Lovallo, M.; Cao, D.T.; Nguyen, L.M. Shannon Entropy Analysis of Reservoir-Triggered Seis-micity at Song Tranh 2 Hydropower Plant, Vietnam. Appl. Sci.; 2022; 12, 8873. [DOI: https://dx.doi.org/10.3390/app12178873]

79. Saar, M.; Manga, M. Seismicity induced by seasonal groundwater recharge at Mt. Hood, Oregon. Earth Planet. Sci. Lett.; 2003; 214, pp. 605-618. [DOI: https://dx.doi.org/10.1016/S0012-821X(03)00418-7]

80. Talwani, P.; Acree, S. Pore pressure diffusion and the mechanism of reservoir-induced seismicity. Pure Appl. Geophys.; 1984; 122, pp. 947-965. [DOI: https://dx.doi.org/10.1007/BF00876395]