Nonlinear Dynamics of Seismicity and Fault Zone Dynamics Around Large Dams: The Case of Enguri Dam, Caucasus

The 271 m high Enguri arch dam, still one of the highest arch dams in operation in the world, was built in the canyon of the Enguri river (West Georgia) in the 1970s. It is located in a zone of high seismicity (MSK intensity IX) and close to the Ingirishi active fault. The high seismic and geodynamical activities together with the large number of people living downstream of the dam make the Enguri dam a potential source of a major catastrophe in Georgia. Thus, the Enguri Dam with its 1 billion cubic meter water reservoir should be under permanent monitoring. At the same time, this area is an amazing natural laboratory, where one can investigate both tectonic and geotechnical strains and processes as well as their response to the lake load–unload impact (i.e., the reaction to a controllable loading of the Earth crust). This is an important scientific issue, connected with a fundamental problem of reservoir-induced earthquakes as well as with environmental geotechnical problems related to the safety of large dams. Application of nonlinear dynamics methods allows for the division of events, ordered by reservoir water regular strain impact from the background seismicity.


Introduction
Monitoring of strains and seismic activity in the area of a large dam is a unique tool for understanding the intimate connections between earthquake generation and man-made regular quasi-periodic strains in the Earth created by seasonal water load-unload in the reservoir. We can consider the area of large dams as a natural laboratory, providing the possibility of studying the seismic process in almost controlled (repeated) conditions.
The 271 m high Enguri arc dam (still one of the largest in the world) was built in the canyon of Enguri river in West Georgia. It is located close to the Ingirishi active fault system, in a zone of high seismicity, intensity IX by MSK scale. The volume of the lake at Enguri dam is 10 9 cubic meters and the high water level in the lake varies seasonally by 100 m, which means that Enguri reservoir can activate reservoir-triggered seismicity (RTS). The dominant tectonic feature of the region is the active east-west-oriented Ingirishi fault, located to the north of the dam; its branch fault crosses the foundation of the Enguri dam [1].
Taking into account the high potential danger of the object, a geophysical monitoring system was organized even before construction works for providing secure exploitation of the large Enguri dam. Due to a high seismic activity of the region (Figure 1), the seismic station's network was installed in the area of Enguri dam also well before its construction, with the aim of studying possible reservoir-triggered activity [2]. The monitoring system of Enguri Dam and its foundation includes network of tiltmeters, piezometers, and reverse plumblines in the dam body [1], a weather station, water level gauge for monitoring water level in the lake, as well as a complex of strainmeters and tiltmeters, installed in the dam body and its foundation [3]. The problem of human-induced earthquakes, including RTS, became a reality in the last decades [4][5][6]. The RTS pattern in the Enguri area should depend on the water level (WL) variation regime in the lake [7,8]. The main goal of the paper is to apply new methods of complexity analysis in order to assess, in a quantitative way, the correlation between WL variations and local seismicity and define the scale of man-made activity on the local (natural) seismicity pattern.

Data
The branch fault of the main Ingirishi fault crosses the foundation of Enguri dam, thus poses a hazard to its safety. In order to permanently monitor the fault behavior, two years before the first filling of the reservoir, in December 1974, a quartz strainmeter crossing the fault zone (FZ) was installed in the adit, located 100 m downstream from the foundation of the dam. The strainmeter's fixed and free parts are located on the intact rocks on the opposite sides of the FZ and are separated from this 10 m wide zone by a 5 m distance (the full length of the quartz tube is 22.5 m). This means that the device records displacement of the intact blocks, divided by the fault zone in the normal to the fault plane direction, so it shows the fault zone's extension-contraction. The free end of the tube is equipped with a photo-optical recording system [3]. The displacements' sensitivity of this system is of the order of 0.18 μm/mm, which allows also a tidal component of the fault zone strain to be recorded. At present, the laser system (Laser model R-39568, Green HeNe Laser, 633 nm and Laser Position Sensor OBP-A-9L) doubles the photo-optical registration. The laser is attached to the free end of the same quartz tube. Sensitivity of the strainmeter with the laser sensor is one μm/mm. The earthquake time series (ETS) for the Enguri area from 3 January 1974 to 31 December 2016 was compiled using catalogs from the Institute of Geophysics and International Seismological Centre. Our study area includes events located at the distances of 50 or 100 km from the lake. The completeness magnitude (CM) for the whole used catalog was around M 1.7 (Figure 1), but in some cases we confine ourselves by magnitude 2.2 for confidence, as in some periods the CM value increased to M2.2 due to non-stable functioning of the national seismic network.
In Figure 2 we present a spatial distribution of seismicity in the Enguri Dam region within 100 km distance from the dam with two magnitude gradations, from 0 to 5 and from 5 to 7, as well as the scheme of active tectonic faults, according to [9].   Figure 3, the dam area experiences stresses of different origin, acting on the different time scales, from decades to months and days. Actually, the object under study is a natural large-scale laboratory for investigation of geotectonic, man-made, and environmental impacts on the fault zone deformation. The summary contributions of these processes are reflected in the time series of fault zone strain. It is evident that the fault dynamics reflect the joint influence of two main factors: One leads to piecewise linear (in time) displacement (trend component) and the other one to quasiperiodic oscillations, decorating the main trend. where the coefficient a, the slope of the linear component of the FZE or the strain rate, differs from one period to another ( Table 1). As the trend component with the same strain rate was recorded even before dam construction and lake filling, we attribute it to the long-term regional tectonic stress action.
At the same time, the fault zone extension rate (FZER) changes significantly with time, reflecting action of some non-stationary factors. In Table 1 we show the periodization of FZE behavior following the pattern of data evolution according to Figure 3, taking into consideration both components of strain-tectonic and anthropogenic. The most probable source of quasiperiodic changes in the FZE dynamics in Enguri area is the variation of the water load in the lake. We elucidate six periods with appreciable differences in the WL regime: (i) May 1974-June 1978, period before water fill, which we consider as a reference,    To separate out more clearly the man-made impact, we present in Figure 4b the detailed seismic rate during the first 200 months after January 1974. We present the corresponding water level regime in Figure 3: the recharge began in April 1978 (arrow 1) and WL abruptly rose to 100 m in November 1978 (arrow 2). Almost simultaneously the abrupt compaction of the fault zone crossing the foundation of dam by approximately 90 μm was registered by the strainmeter installed on the fault (Figure 3). Almost a year later, in December 1979, series of EQs of magnitudes from 3.7 to 4.3 occurred close to the reservoir. These effects follow (with a year lag) the time of WL rising to the critical high of 100 m, when, according to existing data [5,7,8] the water load can generate reservoir-triggered seismicity (RTS).
In order to better resolve the seismic events related to filling and exploitation of the dam reservoir, we used the cellular approach [10], namely, we plotted separately the time sequences of all registered EQs (all EQs of magnitudes M > 1) in the near to dam zone, in the radius R = 50 km from the dam, and in the larger area, in the radius R = 100 km (Figure 5b). In Figure 5 we show the same data for the EQs of magnitudes M ≥ 2.2. Considering Figure 5, we can conclude that the EQ statistics in the near zone R = 50 km (Figure 5a) is dominated by events connected with reservoir impoundment and the swarm of EQs with ∆ of the order of several km (compare with Figure 4b), when the major part of seismic activity in the larger zone R = 100 km (Figure 5b) is due to the relatively strong remote events of magnitude M5.2-5.9 with ∆ of the order of 80-100 km, which are too far from Enguri dam and belong to the class of regional tectonic events (see Figures 2 and 3b). It follows that, in order to distinguish RIS events, it is better to analyze the seismic catalog in the near zone (R = 50 km), where seismic events dominate, located close to Enguri dam. This restriction does not work for analysis of complexity, especially when we analyze waiting times of EQ, because RIS is characterized by a quasiperiodic recurrence property due to regularity of reservoir load-unload. In turn it means that the role of random seismic events when studying regularity in waiting times, even at the distance 100 km, is relatively small.

Methodology
The earthquake time series (ETS) presented in Figure 4a,b presents a complex mix of background seismicity, characteristic for the seismotectonics of the test area with a seismic response to the lake impoundment and, further, to WL quasiperiodic regulation. To single out the dynamical patterns of the seismic data sets connected with WL variation, we used new effective methods of complexity analysis [11] applied to magnitudes and waiting times of earthquake time series (ETS). The complexity analysis allow recognition of periods with different levels of ordering/determinism, which we connect with transition in WL regime from disordered (1978)(1979)(1980)(1981)(1982)(1983)(1984) to more ordered (1985)(1986) and finally, to quasi-periodic loading (1986 until now). We give a short description of these methods below.
Earlier several studies were devoted to complexity analysis of seismic regime in the Enguri dam area, namely to variation of the phase diffusion coefficient of phase differences between daily released seismic energy and water level daily variations [12,13], Visibility Graph Analysis [14], and Singular Spectrum Analysis [15]. In the present paper we analyze recurrent patterns of local seismicity, using such methods as recurrence plots, detrended fluctuation analysis, and Lempel and Ziv complexity measure.
Recurrence Plots (RP). Recurrence Plots allow visualization of the recurrent behavior of a dynamical system by plotting the arbitrary close (after some time lag) states in the two-dimensional projection of the high-dimensional phase space trajectory [16,17]. The recurrence of the same state after some time lag is plotted on the square matrix, where both axes represent time, by zeros and ones or by differently colored dots. The time lag between recurrent points i and j of the trajectory is defined as the threshold time interval (threshold distance) εi . The RP revealed some structural patterns, which are different for different degrees of determinism in the phase space of the system.
Detrended Fluctuation Analysis (DFA). Long-range time correlations in the investigated data sets were assessed by the method of detrended fluctuation analysis (DFA) [18,19]. Method of DFA permits the detection of long-range correlations embedded in a nonstationary time series through calculation of a quantitative parameter-DFA scaling exponent. This analysis technique is widely accepted and often used for different types of time series including geophysical data sets [14,15].
The basics of DFA are well known and described in series of often cited articles, so we will just briefly discuss its main steps. At first, a given time series of N samples is integrated. After, the integrated time series is divided into boxes of length n, and in each box the polynomial local trend is calculated and removed. Then N/n mean squared residuals-detrended fluctuation functions (F(n))-should be calculated for each box of size n.
[ ] Since F(n) increases with the box size n, in case of fractal or self-similar properties of analyzed data, a power-law behavior α n n F) ( can be revealed. If a power law scaling exists, the F(n) vs. n relationship, in a double logarithmic fluctuation plot, will be linear or close to linear and the scaling exponent can be estimated. If the scaling exponent 5 . 0 = α , we deal with the uncorrelated dynamics of random walk type (Peng et al. 1994;Liu et al. 1999). In this case, the time series is identical to white noise. If is different from 0.5, then the time series is regarded as long-range correlated or anticorrelated, with [18,19]. The scaling exponent is considered as an indicator of the nature of the fluctuations giving the information about the long-range power law correlation properties in the analyzed data sets. DFA can be accomplished for different orders of the polynomial fitting in order to eliminate trends of certain origin.
Recurrent quantification analysis. In order to further quantify changes in dynamical structure of analyzed data sets, we used the Recurrence Quantification Analysis (RQA) approach [17,20,21]. In general, RQA is a quantitative extension of the recurrent plot (RP) construction method, which is based on the fact that returns (recurrence) to the certain system condition or state space location is a fundamental property of any dynamical system with quantifiable extent of determinism in underlying laws [16]. In order for RQA calculations to be successfully fulfilled, at first the phase space trajectory should be reconstructed from the given scalar data sets, the proximity of points of the phase trajectory should be tested and marked by the condition that the distance between them is less than a specified threshold [16]. In this way, a two-dimensional representation of the recurrence features of dynamics embedded in high-dimensional phase space can be obtained. Then small-scale structure of recurrence plots can be quantified [20][21][22][23][24]. RQA technique quantifies visual features in a N N × distance matrix recurrence plot and defines several measures of complexity. Exactly, RQA provides several measures of complexity based on the quantification of diagonally-and vertically-oriented lines in the recurrence plot. In this research, we present one of such measures (%determinism), which often is used to reveal changes in the extent of regularity in analyzed data sets.
Lempel and Ziv complexity measure. The Lempel and Ziv algorithmic complexity (LZC) calculation [25,26] is another often used method for quantification of the extent of order in analyzed data sets of different origin. LZC is based on the transformation of a given data sequence into a new symbolic sequence. For this original data are converted into a (0, 1) sequence by comparing them to a certain threshold value (usually median of the original data set). Once the symbolic sequence is obtained, it is parsed to obtain distinct words, and the words are encoded. Denoting the length of the encoded sequence for those words, the LZ complexity can be defined as where L(n) is the length of the encoded sequence and n is the total length of sequence [26]. Parsing methods can be different [26,27]. In this work, we used the scheme described in [27]. Data sequences with a certain regularity are less complex, and the LZ complexity increases as the sequence grows in length and irregularity. In our case, the sequence length is constant and LZC depends only on the level of regularity.

Results
Results of Complexity analysis. In this study, we analyzed two types of ETS related to the seismicity of Enguri area: The waiting times and the magnitude sequences to study the changes possibly linked with dynamical changes in the RTS characteristics of the investigated seismic area related with the loading regime of the dam.
As we mentioned earlier, though the number of events in the far zone can be spoiled by the EQs, not related to the reservoir-induced strain, the waiting times (WTs) distribution is less sensitive to (background) random events, even those relatively far from the dam (e.g., for R = 100 km).
Results for the near zone R = 50 km. In this case (for R = 50 km), we included into complexity analysis, namely RP, namely the Visual Recurrence Analysis methods [28] and DFA methods, EQs below representative magnitude in order to fulfill the condition of used methods-to have at least 500 events. The distinct transitions from a less regular to more regular pattern in seismic regime occurred in 1985 to 1986 (a), when the WL change in reservoir became quasiperiodic (i.e., at the day mark 2000 in Figure 6b, which corresponds to the year 1985 (b); see also Figure 3). Note the light yellow diagonal lines in (a) after 1986, which is a mark of quasi-periodicity in RIS. Similarly, the DFA analysis ( Figure 6c) points to the beginning of ordering in seismic events after 1980 and the transition to strong recurrence regime after 1986.
Results for M2.2 the far zone R = 100 km. We used DFA, RP, RQA, and LZC methods for analysis of interevent times and magnitude data sequences from the Enguri seismic catalogue (1974-2017) involving 913 events above M2.2 that occurred within 100 km distance from the dam.
In Figure 7, we present results of DFA exponent calculation of waiting time sequences. Calculation was done for 500 data length windows shifted by one data. DFA exponents for interevent times (Figure 7) indicate a gradual DFA exponent increase toward the period of reservoir water level periodic variation. A stronger increase took place in the period starting from 1984 and lasts till 2017. According to these results under the influence of water level periodic variation in the reservoir, long-range correlation clearly increased in the time distribution of earthquakes, while the magnitude distribution of earthquakes is characterized by slight or negligible changes in long-range features just at the beginning of observation period and after the 200th window. The DFA for the magnitude distribution of earthquakes is characterized by slight or negligible changes in long-range features.  We also applied RQA and LZC methods to waiting times and magnitude data sequences from the Enguri seismic catalogue . Calculation was done for 500 data length windows shifted by one data (Figures 8-10).

Discussion: Nonlinear Dynamics Patterns in Seismicity and Water Level Variations in Enguri Lake
Let us consider results, obtained by different complexity analysis approaches: The methods abbreviation with a subscript m for magnitude time series and by a subscript wt for the waiting times, for example, (RQA)m and (RQA)wt.

•
We can elucidate the transition in seismic regime by the intensive increase of the exponent for the waiting times (DFA) wt after 1984.

•
Using the RQA approach, we reveal drastic changes in ETS dynamics for (%DET)m after 1986 and for (%DET)wt around 1986 and 2004, when the waiting times became maximal and stable. • Using LZC, we see that there are no changes in (LZC)m, but (LZC)wt undergoes drastic changes in waiting times dynamics around 1986 and 2004.
Resuming, we can mark the dates of significant changes in the seismic time series dynamics around 1986 and in 2004 by both RQA and LZC methods.
Results of calculations presented in Figures 6-10 convinces us that the changes that occurred in waiting times data sets are much stronger than in magnitude sequences. At the same time %DET of the waiting times sequence essentially increases and LZC noticeably decreases: Both these effects point to the growth of order (recurrence) in ETS in the period 1985-1986. According to these results under influence of water level periodic variation in reservoir, long-range correlation clearly increased in earthquakes time distribution, while earthquakes magnitude distribution is characterized by slight or negligible change in long-range features.
Comparing dates of WL regime change and fault zone deformation patterns with transitions in the ETS dynamics, we can conclude that the transition in RQA and LZC in the period 1984-1986 is connected with the beginning of the quasiperiodic load-unload process of the reservoir (see Figure 3). Note, that from 1985 to 1986 the strain in the fault zone under the Enguri Dam also reveals quasi-periodic decoration of the summary strain line (Figure 3).
Thus, in the last period, beginning from the period 1985-1986, the dynamics of local seismicity, especially, waiting times is much more ordered due probably to synchronization of seismic activity with the regular pattern of WL variation. This conclusion is confirmed by our earlier work where we carried out analysis of Enguri area seismic activity using the singular spectrum analysis (SSA) technique in order to investigate the relationship of local seismicity with the reservoir water variations [12][13][14][15]. We revealed the dominant one-year period in seismicity, which corresponds to seasonal load-unload of the Enguri dam lake; this period was absent in ETS of the area in the reference period before lake impoundment.

Conclusions
On the basis of the recurrent plots, recurrent quantification analysis, and Lempel-Ziv complexity analysis carried out on interevent and magnitude sequences of the Enguri area seismic catalogue, we conclude that influence of water level periodic variation makes time distribution of local earthquakes more regular (synchronized with water level variation), compared to the period without such weak periodic influences. This means that nonlinear dynamics methods are effective in detection and quantitative analysis of reservoir-induced seismicity near large dams, as they make it possible to divide events, ordered by the impact of reservoir water regular strain from the background seismicity.

Conflicts of Interest:
The authors declare no conflicts of interest.