Slow slip-driven seismic signature triggered by water-reservoir impoundment

One of the most important and widely used renewable energy sources is hydroelectric energy produced via Water Reservoir Impoundment (WRI). WRI can trigger strong earthquakes under favourable geological conditions. Thus, the socio-economic impact of reservoir triggered seismicity is very significant. Although many studies have investigated the relationship between the pore pressure changes due to WRI and the observed seismicity, hydromechanical models that explain the observed processes are rare. Here, we investigate the role of hydromechanical interactions during fault deformation to understand earthquake swarm bursts under pore pressure changes due to WRI. As a natural laboratory, we selected the Song Tranh 2 Reservoir in Vietnam. Because the analysed triggered seismicity has swarm characteristics, our work contributes to the further investigation of the physical mechanisms responsible for earthquake swarms and their relationship to slow slip. We conclude that the small high-frequency seismic swarms accompanying WRI are driven by slow slip along a fault; they occur due to the temperature-controlled frictional fault heterogeneity, and their rate and magnitude depend on the sizes of these heterogeneities. Swarm earthquakes are the effect of slip acceleration on the seismic radiation level. The nucleation fronts expand the nucleation regime and may transition into stronger earthquakes. These results provide insights into the physical mechanisms of seismic processes triggered by WRI, which may have implications for assessing the seismic hazards associated with hydroelectric energy production. Teaser Seismic swarms accompanying water reservoir impoundment are driven by slow slip controlled by frictional fault heterogeneity. Introduction The basis of the European Union's climate and energy policy is to reduce greenhouse gas emissions, increase the proportion of renewable energy sources, and increase energy efficiency. Thus, geo-energy projects are an essential element in society's sustainable development strategy. One renewable energy source is hydroelectric energy production via the usage of Artificial Water Reservoirs (AWRs). This technique is also crucial for natural water retention to control floods and mitigate droughts (1, the EU Water Framework Directive and the EU Floods Directive: https://ec.europa.eu/environment/water/index_en.htm). Hydropower represents about 17% of the total electricity production worldwide (source: International Energy Agency).


Introduction
The basis of the European Union's climate and energy policy is to reduce greenhouse gas emissions, increase the proportion of renewable energy sources, and increase energy efficiency. Thus, geo-energy projects are an essential element in society's sustainable development strategy. One renewable energy source is hydroelectric energy production via the usage of Artificial Water Reservoirs (AWRs). This technique is also crucial for natural water retention to control floods and mitigate droughts (1, the EU Water Framework Directive and the EU Floods Directive: https://ec.europa.eu/environment/water/index_en.htm). Hydropower represents about 17% of the total electricity production worldwide (source: International Energy Agency).
Water reservoir impoundment (WRI) is known to trigger earthquakes under favourable geological conditions. The socio-economic impact of Reservoir Triggered Seismicity (RTS) is very significant (e.g., 2, 3). Globally, hundreds of sites have been reported to have experienced RTS (4). The largest WRI triggered earthquakes, with magnitudes of M7.1 and M6.5, occurred on 17 August 1959 and 15 May 1995 in the USA and Greece, respectively (5,6). Another large WRI triggered earthquake, with a magnitude of M6.3, occurred on 10 December 1967 in the vicinity of Koyna Dam near India's west coast. In this region, a 3 km deep pilot borehole has been implemented as a prelude to setting up a borehole laboratory for near field studies of RTS (3,7). There are also earthquakes for which the genesis-associated with either tectonic loading or WRI-is still under debate. For example, the M7.9 earthquake in Sichuan, China, which occurred on 12 May 2008 after the nearby Zipingpu Reservoir was filled, claimed over 80,000 lives. A detailed review of RTS can be found in Gupta (3).
The different temporal seismic responses to AWR impoundment have been observed and classified as follows: rapid response of seismicity, which starts with the filling of the reservoir; delayed response of seismicity, which occurs after a few cycles of the water level increasing and decreasing; and continued or protracted seismicity, which lasts throughout the reservoir exploitation (8)(9)(10)(11). The related seismicity is often swarm type seismicity, in which multiple similar-sized earthquakes occur over periods of minutes to days. Earthquake swarms occur in various tectonic settings as a result of gas or fluid flow [12][13][14] and as a result of aseismic fault slip (e.g., 15,16). Recent studies using highquality geophysical data have provided insights into the physical processes that trigger earthquake swarms (e.g., 17, 18), but the details of this mechanism remain unclear.
Here, we investigate the role of hydromechanical interactions during fault deformation to understand swarm bursts in the time evolution of the seismic moment release under pore pressure changes due to AWR impoundment. As a natural laboratory, we selected the Song Tranh 2 (ST2) AWR in Vietnam (Methods; Fig. 1). Hydropower projects in Southeast Asia have accelerated in the last 10 years, because of the growing need for electricity in developing countries like Vietnam, Laos and Cambodia; many projects are currently in development or planned. The impoundment of reservoirs triggered the strongest anthropogenic seismic event in Vietnam: a M4.9 in Hoa Binh Province, Northern Vietnam, in 1989 (19). The selected ST2 reservoir exhibits typical seasonal water level changes and is located in a tectonically complex area. Such hydro-tectonic conditions are common for AWR projects; thus, the results of our analysis are of broad applicability in terms of seismic process evolution.
Although many studies have reported a relationship between the pore pressure diffusion due to AWR impoundment and observed seismicity, hydromechanical models that explain the observed processes are limited (3). Many studies have presented results on rock samples and in situ experiments with the fluid injection (e.g., 20-22); however, no analogue for the seismicity triggered by water reservoir impoundment has been reported.
Here, we compile the results of recent studies involving fluid-rock interaction experiments and numerical modelling of frictional sliding, utilize these conclusions to explain the hydromechanical mechanism of earthquake triggering, and provide a model for seismicity triggered by reservoir impoundment. Because the analysed RTS has swarm characteristics, our work also contributes to the further investigation of the physical mechanisms that cause earthquake swarms and the exploration of the possibility that swarms are associated with slow slip.

Study area
To conduct this analysis, we launched a seismic monitoring campaign, during which we recorded more than 7000 earthquakes with magnitudes ranging from −0.6 to 3.6 in the ST2 region from August 2013 to June 2017 (see Methods for details). This experiment is one of three projects targeting AWR impoundment seismicity in the ST2, Lai Chau, and Hue regions of Vietnam (23). The ST2 Reservoir's volume is 740 million cubic meters, and the water level varies from 140 to 175 m. The seismic activity in the reservoir region began in 2011, soon after the reservoir was filled. The natural seismic activity in the vicinity of the ST2 AWR is low. This area represents the transition zone between the Indochina Block and the East (South China) Sea, with right-lateral strike-slip movement and minor seismicity (24). Only 13 events were reported in this area between 1775 and 1992. One M4.7 earthquake (in 1715) occurred near the planned hydropower reservoir site (25). The seismic activity in the region increased significantly after the reservoir was filled in November 2010, and it has continued to the present day. The seismicity in the ST2 area reached more than 250 events per month (M0. 9

Stress changes due to water reservoir impoundment
The critical information needed to mitigate potential seismic hazards associated with AWR is knowledge of the potential factors governing the seismic moment release in response to the operational parameters of the AWR. It is known that the impoundment of a reservoir can trigger earthquakes due to (i) the elastic stress increase following the filling of the reservoir; (ii) the increase in the pore fluid pressure, , in saturated rocks, which is usually caused by a decrease in the pore volume due to compaction in response to an increase in the elastic stress; and (iii) the pore pressure changes, ∆ , related to fluid migration (3 and references therein). Talwani & Acree (27) pointed out that the pore water's effect on the earthquake process may also act as a chemical process due to stressaided corrosion. There is evidence to suggest that the diffuses along pre-existing fractures (28), or it diffusion can be associated with new crack propagation through stress corrosion (29). Moreover, Talwani & Acree (27) suggested that the mechanical effects of the control the spatial and temporal pattern of the RTS. In contrast, the actual onset of seismicity may be influenced by the water's chemical effect in reducing friction in clays.
The fault-stability study of Gahalaut et al. (30) of the ST2 AWR revealed that the seismicity, which started in 2011 after the ST2 Reservoir's impoundment, was a rapid response due to the load of the reservoir and then continued in time due to pore pressure diffusion. They concluded that the load of the ST2 reservoir's water destabilized nearby faults soon after its impoundment and triggered earthquakes. The range of the pore pressure changes due to the compression caused by the load did not exceed 5 • 10 −3 for the nearby faults. The pore pressure diffusion further expanded the area of destabilization and reached 20 • 10 −3 and 0.8 • 10 −3 in the northern and southern seismicity clusters, respectively, triggering two major earthquakes with local magnitudes of 4.6 and 4.7 on 22 October and 15 November 2012, respectively (31).
Since Gahalaut et al. (30) demonstrated that the magnitude of the stress changes arising from the poroelastic coupling between the load of the reservoir's water and the rock matrix is relatively small and diminishes with distance, we consider the ∆ due to compression and diffusion after two years from that period herein. For the analysed dataset, the hydraulic diffusivity, c, was estimated to be around 15 2 / . Based on this information, the estimated Green's function solutions of Kalpna & Chander (32) were used to estimate the induced pore pressure throughout a half-space with a given permeability (Methods). We calculated the ∆ at the time and location of the earthquakes. The range of the diffused ∆ in the ST2 region did not exceed 0.4 and 0.2 in the northern and southern seismicity clusters, respectively. The mean depth of the seismicity in the northern cluster is 3 km, while it is 6 km in the southern cluster. Fig. 2 shows that the pore pressure's response to the water level fluctuations was delayed by about three and six months in the northern and southern seismicity clusters, respectively. A similar delay in the ∆ caused by the water reservoir impoundment was observed by Gahalaut et al. (30) at the beginning of the ST2 operation. The results of the high-resolution seismic tomography of the ST2 region (details are provided in the Methods) also indicate the presence of fluids. Fig. 3 presents the ratio of the velocities of the shear and compressional waves, , and the compressional wave velocity, , along cross section A-A', which is delineated by the location of the strongest earthquakes in the analysed period, at the depths selected to satisfy the high data resolution criterion. The values are moderately high (1.73), indicating that there is increased pore pressure in this region. There are low ratios beneath the westernnorthern part of the water reservoir at a depth of 1.5 km, and there is a simultaneous increase in . This decrease in and increase in continues in the deeper layers, with a shift to the right following the water reservoir's extent. This may imply that the fluids leaking from the part affected by the overpressure due to the compression of the load of the reservoir's water migrate downward along the underlying shear zones, which function as fluid pathways. These migrating fluids may reduce the strength of the crustal rocks and enhance the aseismic slip (e.g., 33). The ratio of the seismic velocities decreases with increasing P-wave velocity, indicating that the rocks are water-saturated (e.g., 34,35). At a depth of 3.5 km, the minimum and maximum correlate with the location of the three strongest earthquakes (turquoise circles in Fig. 3), which occurred below 3.5 km. The low and high body above the zone of higher (solid blue curve in Fig. 3) may act as an impermeable barrier, trapping fluids in this zone and generating fluid overpressure (e.g., 33, 36).
We calculated the moment tensors of 148 events with local magnitudes of 0.9 to 3.6 through inversion of the P-and S-wave amplitudes in the time domain (see Methods for details). To perform a temporal analysis of the stress field changes, we applied a stress inversion to the focal mechanisms of the moving window of 15 earthquakes with a moving step size of 1 event. The setup of the stress inversion was selected to detect the small variations in the relative stress magnitude and to ensure the requirement of a certain variety of focal mechanisms was met. The stress inversion analysis revealed an extensional environment characterized by a subvertical, effective maximum principal stress with a strike/plunge of 1 = 289°/79°. The 1 axis is oriented vertically, and its magnitude can be assessed based on the overburden weight (e.g., 40), i.e., 1 ≅ , where is the density of the crustal material assumed (2800 / 3 ), is the acceleration due to gravity, and is the focal depth. Assuming that the plane is critically stressed, the least principal stress, 3 , can be estimated as ( 1 − )/( 3 − ) = (( 2 + 1) 1 2 − ) 2 (e.g., 41,42), where is the coefficient of friction estimated to produce the highest overall instability of the fault (Methods). Finally, using the -value, the absolute magnitude of the principal stresses and the normal, , and shear, , stresses on the individual fault plane can be estimated (Fig. 2). Based on these calculations, the average differential stress was 70 MPa, which varied from 30 MPa to 140 MPa in the two short summer periods in 2014.  (43,44). The timing of the stiffness drops is shown in light orange. The dashed yellow and blue rectangles mark the two different parameter signatures of the process indicating potential different mechanisms responsible for the stiffness drop. The first parameter leads to the porosity loss resulting in the pore pressure increase (compression, shown by yellow rectangle), and the second leads to the volume expansion resulting in the pore pressure decrease (dilatation, shown by blue rectangle). All the recorded seismicity is marked by the grey circles.

Fault strength evolution with fluid pressure variations
The variations in ∆ within the fault zone influence the elastic stiffness of the loading system, . If the value of , which describes the fault strength loss, ∆ , with slip, , is smaller than the critical fault rheological stiffness, , and thus, slip instability may occur (e.g., 45,46). is defined by the effective normal stress, ′ = − , and the frictional constitutive properties of the fault: the friction rate parameter of the rate-and state-dependent constitutive law (R&S) (47)(48)(49), ( − ), positive for slip instability, and the critical slip distance, , over which the strength breaks down during earthquake nucleation. If ~, the fault is conditionally stable although dynamic instabilities may arise when the elastic unloading of the surrounding medium matches the frictional weakening rate of the fault, allowing slip acceleration at slow speeds due to the small force imbalance (e.g., 50,51). For > , only aseismic creep is possible. An increase in the fluid pressure reduces , favouring stable sliding rather than earthquake slip (52). Assuming that is 10 −3 and corresponds to natural earthquakes (e.g., 43,44), then on average, ( − ) has to be smaller than 0.008 (green line in Fig. 2) to satisfy Eq. (1). The evolution of the fault's strength indicates cyclic strength loss. The periods of the total stiffness and shear traction evolution are framed in light orange in Fig. 2 to better visualize the periodicity of the process. Each curve in the cycle was mapped using the parameters of the seismic frictional instabilities (seismic events). We observed two signatures for this process. The most often observed is that the fault stiffness and shear stress decrease with time, ∆ decreases with time, and the evolution of ( − ) shifts slightly towards neutral velocity behaviour (blue dashed rectangle in Fig. 2). The correlation coefficient between ( − ) and ∆ is 0.4, and it is statistically significant. However, ( − ) usually remains at the neutral velocity level. The second behaviour also indicates strength loss, but ∆ increases with time and the degree of velocity weakening slightly increases with time (yellow dashed rectangle in Fig. 2).
The temporal evolutions of the estimated quantities, such as the stiffness, shear stress, friction, slip velocity (Eq. (2)), and slip for the seismic events representing one cycle of the process, indicate the underlying slower evolution of the weakening behaviour (Figs. 4A, B), which is similar to the temporal evolution of the dynamic traction within the cohesive zone during the propagation of an earthquake rupture (e.g., 53 and references therein). Moreover, this slip weakening behaviour agrees with the fault slip behaviour observed during slow-slip events (S-SE) generated in fault gouges (e.g., 37,51).
To better understand the mechanical processes that are at the origin of such slips, we stacked the curves framing the quantities for all of the cycles (Fig. 4C) and analysed the relationships between the evolution parameters and the pore pressure changes. There are six regular distances over which the weakening behaviour was completed (Fig. 4C), which is considered to be a proxy of the local critical weakening distance of the slow-slip events, − . The values of − depend on the relationship = ∆ ′ ( Fig. 4D; colours of dots are related to the stacked group of cycles). This logarithmic dependence is strong, and for 17 cases (two cases with low were removed as outliers), the nonparametric Spearman's correlation coefficient is -0.7 and the correlation is statistically significant. The −dependency is inversely proportional to the maximum observed magnitude in the cycle (Fig 4D; colours of stars are related to the stacked group of cycles). The Spearman's correlation coefficient for and the maximum observed magnitude is -0.87 (calculated for 17 cases). According to Dublanchet's (54) study of the impact of the fault's frictional heterogeneity on fault behaviour, our results indicate that the loading of the ST2 fault matches the conditions for the acceleration of the background aseismic slip along the fault, and the associated seismic events are coupled with an increase in the rate of the background slip (regime IV in 54). The required conditions correspond to a neutral velocity regime with velocity-weakening patches larger than the local nucleation length. The fault patch with the accelerated slip becomes larger as the slip rate increases, and this can be seen as subcritical crack growth (55), e.g., growth with a sub-critical energy release rate. One important feature revealed by the numerical simulations of Dublanchet (54) is that the evolution of the background slip rate under such a regime seems to be unaffected by the occurrence of seismic events. Thus, the observed seismic events in the S-SE cycles can be considered by-products of the acceleration of the background slip, which can be regarded as aseismic slip acceleration (54). If the time it takes to reach failure is long enough (e.g., weeks to months), these seismic event swarms may occur as repeating events in the same velocity weakening patches of the fault. In some cases, the nucleation develops as very narrow slip pulses travelling through the fault plane (54).

Slip rate of S-SEs and its relationship with the rate-weakening properties and fault length
Based on the relationship between the rate dependence of the strength of the R&S friction fault and the rate dependence of a single crack's growth rate (56, 57), we approximated the S-SEs slip rate from the slip rate of the networks formed by the fault velocity weakening patches (FN) ruptured within the S-SE cycles. For this purpose, we applied the subcritical fracture growth (SFG) law (58), which describes the crack tip velocity in a subcritical regime. Charles' SFG law can be used in the following form, which describes the observed average increase in the FN length per unit time (59): where is the total length of the FN at time ; is the growth exponent; is a parameter that depends on the stress state; and = 1, … is the number of observed increases in the FN in time period ∆ due to constant loading, which is assumed in ∆ . The FN propagates with a mean velocity of _ in ∆ . The parameters and are determined based on the slope and intercept of the linear regression of the plot of all of the log ( _ ,∆ ) values versus log ( ) values. The value of parameter is controlled by the rate of the applied stress (57,59,60). The rupture length is proportional to the cubed root of the seismic moment, 0 . Assuming that the average stress drop is roughly constant for earthquakes, the slip on the fault scales with the fault's length.
We found that the estimated values for most of the S-SE cycles were positive, ranging from 0.4 to 0.9, indicating that the S-SE rate increases with increasing stress. For 6 of the cycles, the values were negative, which may imply predominantly aseismic behaviour (e.g., 61). The estimated S-SE velocities of the growth of the velocity weakening patches in the cycles with positive values ranged from 10 −3 m s to 10 −2 m s (from Eq. (2), for the rupture length approximated using the cubed root of 0 ). These velocities correspond to the rupture front velocity of the slow slip events, e.g., in the Cascadia region (e.g., 62) and in Mexico (63), as well as to the numerically modelled rupture velocities of spontaneously occurring slow slip events (e.g., 64). This may indicate that the seismic swarms in these cycles were driven by aseismic slow slip.
To determine the fault geometry conditions conducive to generating subcritical growth of seismic events, we used the synthetic catalogues produced by Romanet et al. (65), which were derived from the numerical rate-weakening fault models of the seismic cycle that account for the geometry and stress transfers of faults. In these models, two overlapping faults of the same length are considered to study the effect of the complex stress interactions between neighbouring fault segments. The friction on both faults is controlled by the R&S constitutive friction with aging state evolution (66,67). All of the details of the models, simulations, and results have been reported by Romanet et al. (65). Here, we analysed the evolution of the seismic moment release for all the models according to Eq. (2).
We found that SFG behaviour was only observed under one specific set of conditions. Subcritical events emerge when the faults are long, i.e., the fault length is much larger than the nucleation length (L ≫ L nuc ) with a small relative distance between the rupturing fault segments, and the faults have velocity neutral properties ( a b ≥ 0.9, where a and b are the R&S constitutive friction parameters). These faults host a mixture of rapid slip earthquakes and slow slip events, with a slip velocity in the range of 1 μm/s to 1 mm/s and a rupture velocity of less than 0.001 , where is the shear wave speed (65). For these models, we found that the values range from 0.6 to 0.9. In other cases, where the faults are mild to strongly rate-weakening and have small lengths, there is no relationship between the average rupture velocity and the length of the rupture according to Eq. (2).
These results confirm our interpretations, which were derived based on the study by Dublanchet (54). The numerical simulation-based results imply that the fault segments in the ST2 region have neutral velocity conditions conducive to generating slow slip. Moreover, their geometric properties and the heterogeneity of the frictional conditions result in the growth rate of the velocity weakening patches being dependent on the rupture length.

Indications of slip pulses nucleating within the S-SEs
The seismic events within the S-SEs in the ST2 region are characterized by the radiation of high-frequency seismic waves, which can be linked to slip pulses. Slip pulses are a primary source of the high-frequency seismic waves that radiate during a rupture (64). The frequency of the near-field ground motions depends on the rise time, that is, the duration of the slip at a particular location along the fault. Due to the conditions of the seismic network, we cannot accurately estimate the rise time (e.g., 68). Instead, following the method of Boore (69) and Beresnev (70), we rely on the inverse of the corner frequency to draw conclusions about the approximated source duration ( _ = −1 , where is the corner frequency). We estimated the _ of the seismic events in each S-SE cycle (see Methods for details) and compared them to the _ estimated for other earthquakes induced by geo-resource exploitation. Fig. 5 shows the _ for a particular series of seismic events during the S-SE cycles in the ST2 region and the results obtained for the seismicity induced by other geo-energy production activities. We tested the internal variability of the results for the ST2 region with respect to the variability between earthquakes in the ST2 region and the datasets representing wet (geothermal activity) and dry conditions (mining environment) using analysis of variance (ANOVA). ANOVA identifies any differences in the tested groups' means. We also tested the influence of the local magnitude, , on the approximated rise time. The results of the ANOVA (Methods) and the post hoc tests show that the _ values of the seismic events in the ST2 region are statistically significantly lower than the _ values of the other induced seismicity. does not influence _ .
Heaton (71) reported that short rise times may be linked to earthquakes that are composed of a sequence of small-scale events that are separated by locked regions or to ruptures propagating as a narrow self-healing slip pulse along a fault surface with velocityweakening friction. The theoretical and laboratory studies of dynamic ruptures and the numerical simulations of the evolution of R&S faults with velocity-weakening regions have revealed that the rupture mode depends on the fault's pre-stress and the velocityweakening properties of the faults (72). Pulse-like ruptures may occur along homogenous interfaces when the frictional resistance is considerably reduced at a rupture's tip (e.g., [72][73][74]. Pulse-like ruptures can also result from a complex fault geometry or local heterogeneities (e.g., 64,75,76) and due to normal stress variations due to differences in the material properties across the interface (e.g., [77][78][79]. Moreover, their velocities may range from slow to fast seismic speeds. Beyond these theoretical studies (72), it has been experimentally demonstrated that pulse-like ruptures may transition into crack-like and supershear ruptures (74). Passelègue et al. (80) pointed out that sudden accelerations and decelerations of rupture speeds while transitioning from a sub-Rayleigh rupture to a supershear rupture could play an important role in the generation of high-frequency radiation.
The seismic events accompanying the S-SEs in the ST2 region are primarily characterized by slow rupture velocities, , with mean shear wave velocities, , of 0.2 to 0.8 m/s. Moreover, we found that increases as the ∆ and the water level increase (Figs. 2 and 5). The statistically significant Spearman correlation coefficient for the correlation between the ⁄ ratio and ∆ is equal to 0.6. The values were estimated through inversion of the peak ground acceleration (PGA) and the peak ground velocity (PGV) data by applying the methodology proposed by Boatwright (81) (Methods). This methodology also allowed us to account for the rupture directions, which provide insight into the fault mechanics model presented below. The observed rupture directions are up-and down dip.
The up-dip ruptures may be conditioned by their nucleation at the interface between a slipping aseismic zone beneath a slip weakening seismic zone (82,83), and they usually exhibit slower rupture velocities (84). We noticed that there was a decrease in the number of rupture directions up-dip proceeding the strong earthquakes, indicating significant frictional heterogeneity along the fault interface (a higher resistance barrier); and shortly before the strongest earthquake ( 3.6), there was a significant increase in the number of down-dip earthquakes, which may be linked to the breaking of this barrier.

Interpretation of the fault slip mechanism
The analysed seismic events in the ST2 region provide a window into the slip processes triggered by the water reservoir impoundment. Our observations indicate that the seismic swarms were driven by the progressive acceleration of the slow slip front towards the radiative threshold. The nucleation and evolution of the slow slip is governed by the pore pressure changes caused by the fluctuations in the water level of the reservoir. The rate of the generation of the swarms during the slow-slip cycles corresponds to the rupture velocity of the slow slip events observed in the various seismogenic areas. The accumulation of small ruptures forms a larger nucleation path, and dynamic rupture occurs when the stress concentration ahead of the nucleation crack becomes large enough to drive the crack to failure. This rupturing process is governed by the subcritical rupture growth law.
Our findings also illustrate that swarm evolution is controlled by the frictional heterogeneity of the fault. The presence of dynamic events during the S-SE cycles requires a transition from the conditional stability of the fault path to the velocity-weakening state. We conclude that the S-SEs occur when an increase in fluid pressure promotes a slow pulse-like slip, which passes the velocity-neutral fault segments, becomes faster, and finally reaches the velocity-weakening conditions and generates the observed seismic slip rate. The approximated weakening slip distance of the slow-slip patches is inversely proportional to the ratio of the pore pressure changes to the normal stress. This means that the increase in the fluid pressure along a fault reduces the applied normal stress, which promotes the shorter weakening distance of the S-SEs. This agrees with other laboratory experiments (e.g., 20, 21). These conditions, together with the evolution from velocitystrengthening to velocity neutral, lead to a transition from stable sliding to frictional instability. Subsequently, the slip may spontaneously accelerate and proceed as slow slip events or even accelerate to seismic rupture rates under slip-favourable conditions (65,76,85).
In the considered case of the S-SEs, the pulse-mode rupture requires the slip area to be bounded, and the rupture only propagates along strike. The event is then arrested in the area of unfavourable pre-stress. This geometric limitation and slip arrest may be linked to the frictional heterogeneity and the transition from a velocity-weakening domain to a velocity-strengthening domain, which are controlled by the geometric complexities, curvature, and roughness of the fault ( To further investigate the processes behind the fault weakening, we computed the average fault fracture energy, , which is defined as the frictional work associated with fault weakening (e.g., 90-92) based on the slip evolution in the identified S-SE cycles (the shaded orange area below the slip weakening curve in Fig. 2B): where ( ) is the frictional resisting stress; and = ( ) is the dynamic stress.
These results were compared to the scaling relationship between the fracture energy and the slips derived for a wide range of earthquake sizes and to an elastic-dynamic numerical model that accounts for possible weakening mechanisms, which was developed by Viesca & Garagash (93) (Fig. 6). These authors found two different relationships between and for small and large slips. Furthermore, they concluded that thermal pressurization of the pore fluid by the shear heating of the fault gouge can account for the observed scaling of the fracture energy during small and large earthquakes. For small slips, the scaling of the fracture energy with the slip is similar to 2 , and it is consistent with linear slip weakening under undrained-adiabatic pressurization (dashed line in Fig. 6). For large slips, the scaling breaks down and is replaced by a sublinear scaling of ∝ 2 3 (dotted line in Fig.  6). The heating is effectively reduced with the slip along a plane and the pressurization, which is limited to the slip plane, is conditioned by the hydrothermal diffusivity of the gouge and the width of the actively shearing gouge (93). Such a process leads to a gradual strength drop, and it is less efficient than the exponential slip weakening of the undrainedadiabatic shear in the small-slip case. The estimates of the fracture energy of the S-SEs follow the scaling of the fracture energy with slip attributed to thermal pressurization, and they are larger than those of the swarm events. This may be related to the scale of the fault heterogeneities, i.e., indicating larger sizes (e.g., 54,95). The relationship shown in Fig. 5 implies that the thermal pressurization process took place during the nucleation and growth of the S-SEs in the ST2 region.
Based on our observations, we propose a plausible hydro-mechanical scenario for the slip along the fault in the ST2 region. (i) The pore pressure changes related to the fluid migration induced by the water reservoir impoundment reach the fault and decrease the normal stress to a level suitable for slip nucleation. (ii) During the slip, the gouge is compacted and the and friction increase, which both increases and decreases the frictional strength. Due to the thermal pressurization, the pore pressure increases further, leading to the abrupt acceleration of the slip and causing a significant reduction in the fault's frictional strength at the asperity contacts (frictional heterogeneity along the fault), which is manifested as high-frequency seismic events that are observed as swarms in the case of small weakening patches. If the slip acceleration rate is high, the nucleation length increases, and a stronger seismic event may occur. (iii) Multiple seismic events were observed to break the strengthening barriers (asperities) along the fault and promote large slip. Due to the fault's frictional heterogeneity at the asperity contacts, we observed seismic ruptures in two directions, i.e., opposite to the motion (up-dip) and in the same direction as the slip (down-dip), whose growth is governed by the subcritical fracture growth law. In both cases, the rate and magnitude of the reduction of the normal stress at the rupture tip is different and may drive the slip-like pulses. (iv) During the seismic events embedded in the domain of the slow slip, decreases due to the dilatancy. This process increases the effective normal stress and reduces the seismic slip (52,87,96). (v) The reduction of the slip rate to slow velocities leads to fault compaction and thermal pressurization, and recovers to close to the initial value. The increase in enhances the intrinsic weakening of the mechanism, but when the fault slips due to the S-SE cycle, this mechanism is suppressed by the drop in . (vi) When the fault is slipping, the implied fault dilatancy decreases and the effective friction coefficient, which corresponds to the increase in the shear strength before the next slip. The next slow slip cycle occurs under further continuation of the load of the pore pressure diffusion. A schematic model of this process is presented in Fig. 7. Based on the results of this study, we conclude that the small high-frequency seismic swarms accompanying artificial water reservoir impoundment are driven by slow slip along a fault, which is controlled by the effective normal stress and the rate and state of the friction parameters. The swarms occur due to the fault temperature-controlled frictional heterogeneity. Their rate and magnitude depend on the size of these heterogeneities. The nucleating swarm fronts expand the nucleation regime and may transition into stronger earthquakes. In the study area, i.e., the Song Tranh 2 region, this process triggered five earthquakes with magnitudes of ≥ 3.0 during the analysed period. The strongest seismic events in the ST2 region were associated with high pore pressure levels. This can be explained by the increase in the critical nucleation length, which is necessary for dynamic slip, due to the higher pore pressure or velocity-neutral frictional properties of the fault. However, these strong earthquakes did not occur during the highest pore pressure changes, which can be interpreted as the increase in the critical nucleation length being larger than the length of the entire fault segment. These results provide insights into the physical mechanisms of the seismic processes triggered by water reservoir impoundment, which may have implications for assessing the seismic hazards associated with this energy production technique.

Methods
Seismic monitoring and earthquakes parameters. The monitoring of the seismicity in the ST2 area was undertaken in three phases (97). At the beginning of the reservoir filling, the area was monitored by two seismic stations located away from the reservoir. Phase II of the monitoring began in October 2012 when a five station seismometer network operated by the Institute of Geophysics of the Vietnam Academy of Sciences and Technology (IG VAST) was deployed in the ST2 region. Since then, it has been possible to both record events with magnitudes of less than 2.0 and to determine their epicentre locations more accurately. Phase III began in August 2013 when the IGP VAST and the Institute of Geophysics, Polish Academy of Sciences (IG PAS), modernized the existing local stations and extended the number of stations to 10. At this time, the network was renamed VERIS (ViEtnam Reservoir Induced Seismicity). Five of the stations were equipped with long-period Guralp CMG-6TD instruments (30 s), while the other five were equipped with short-period (1 s) Lennartz LE-3DLite seismometers (26). All of the details about the equipment, catalogue, 1D velocity model, and water levels of the Song Tranh 2 reservoir are available on the IS-EPOS e-research platform (https://tcs.ahepos.eu/#episode:SONG_TRANH). The moment tensors (MT) of 148 events with local magnitudes of 0.9 to 3.6 were obtained through inversion of the P-and S-wave amplitudes in the time domain (98). The initial catalogue of events used to our investigation includes the hypocentre location, focal mechanism parameters (MT components and nodal plane orientations), and local magnitude. Additional seismic source parameters and the static stress drops were calculated for the entire input catalogue via spectral investigation. The analyses were conducted following the methodology described by Boore and Boatwright 99), using the formulas of Brune (100), and assuming a circular fault model (101). During the analyses, we used 3-component seismograms obtained from the stations located around the sources. The original velocity signals were corrected for the instrument response to obtain the real displacement recordings. To obtain the spectral level and corner frequency, each spectrum was examined separately.
Hydraulic diffusivity and pore pressure changes. Following Talwani & Acree (27), we assumed that the diffusion of the pore pressure is the transmitter of the ∆ to the hypocentre depths. The observed seismicity is associated with the pore pressure front. Because the rate of the diffusion can be estimated directly from the distance, L, between the source of the pressure front (i.e., the reservoir) and the location of the seismicity, and from the time delay, t, between generating this front (i.e., filling or draining the reservoir) and the onset of seismicity, the following relationship was used (27): The hydraulic diffusivity c estimated using Eq. (4) is called the seismic hydraulic diffusivity, and it has been found to be a good estimate of the material hydraulic diffusivity (e.g., 103,104). For the analysed dataset for the ST2 region, the c value is around 15 2 / , which is within an order of magnitude of the material hydraulic diffusivity and is in agreement with other values observed for fractured crystalline rocks Our studies show that in the analysed period (05/2014-06/2015), the hydraulic diffusivity c was higher than that calculated by Gahalaut et al. (30). This could be the result of fracturing. Witherspoon & Gale (107) have shown that water flow in fractured rocks principally occurs through joints, faults, and other planar features. This is also in agreement with the observations of other authors, i.e., that pore pressure flow occurs through discrete fractures (e.g., 102, 103, 108). As a result of the reservoir exploitation over time, the hydraulic diffusivity may increase due to the fracturing of the crystalline rocks. The induced change in the pore pressure due to the reservoir throughout the halfspace of a given permeability was calculated by solving the diffusion equation using the Green's function approach with the determined boundary and initial conditions (32). The algorithm can be implemented for the actual time varying water level changes of the load of the reservoir's water. The homogeneous diffusion equation for a porous elastic medium (106,109) is where ∆ is the change in the mean stress; = ∆ + ∆ + is the sum of the stress tensor components; is Skempton's coefficient; and the pore pressure change, ∆ , is the sum of the change in the pore pressure due to the instant compression, ∆ , caused by the load of the reservoir's water and the change in the pore pressure due to the diffusion of the load of the reservoir's water, ∆ ( where , , and ′, ′, ′ refer to the earthquake's location and the point on the reservoir dam closest to the earthquake's hypocentre, respectively; and the , , and axes refer to the north, east, and downward (vertical) directions, respectively. is the Green's function defined as follows: )]. (7) The method used and examples of its application are presented in Kalpna & Chander (32); and their extension to the case with heterogeneity, i.e., different permeability/hydraulic diffusivities, is presented in Gahalaut & Gupta (110).

Magnitude of the relative stress R.
We calculated the magnitude of the relative stress as R = (σ1 − σ2)/(σ1 − σ3), where s1, s2, and s3 are the principal stresses, retrieved using the stress inversion methodology used in the STRESSINVERSE software package (39). STRESSINVERSE calculates the stress orientation based on the focal mechanisms (strike/dip direction/dip angle). We used 1000 noise samplings of the input focal mechanisms to estimate the 95% confidence intervals of the magnitude of the relative stress. We performed a temporal analysis of the changes in the stress field. We applied the stress inversion to the entire dataset. Then, we performed a stress inversion of the focal mechanisms from the moving windows of 15 earthquakes using a step size of 1 event to detect the small variations in the magnitude of the relative stress. The number of events in the moving window was selected to achieve balance via a trade-off between the discrimination of the different values of ∆ and the requirement of a certain variety of focal mechanisms.
ANOVA results. Analysis of variance (ANOVA) provides a global assessment of the statistical difference between the considered independent means (111). In ANOVA, the observed variance in a particular variable is partitioned into components, which are attributable to different sources of variation. In its simplest form, ANOVA provides a statistical assessment of whether the means of two or more populations are equal, and the alternative hypothesis is that at least one of the populations' means is different. The F-test statistic is used to draw conclusions about the null hypothesis in ANOVA. In the present study, we tested three hypotheses: 01 states that the mean values of the _ calculated for the earthquakes in the S-SE cycles in the ST2 region are equal for all of the S-SE cycles; 02 states that the mean values of the _ calculated for the S-SE cycle earthquakes and for other earthquakes induced by mining and geothermal energy production are equal; and 03 states that the magnitudes of the events do not influence the mean values of the _ . The tested null hypothesis is rejected, and the alternative hypothesis is accepted if the p value, which is the highest level of significance at which the null hypothesis can still be rejected, is less than the assumed significance level, usually 0.05. ANOVA also allows us to estimate and test interaction effects. The interaction effects represent the combined effects of factors on the dependent measure. When an interaction effect is present, the impact of one factor depends on the level of the other factor.
To follow the assumption of ANOVA about the variance's homogeneity, we removed one S-SE cycle from the first test. For the remaining 18 S-SE cycles, the p value of the F-test statistic was 0.24, indicating that there was no significant differences in the _ values of the seismic events in the ST2 region. For 02 , the p value of the F-test statistic was 0.0, indicating that there was a significant difference in the _ values of the seismic events recorded in the ST2 region and those of the other mining and geothermal energy production induced seismic events. To assess the effect of the magnitude, we selected the data with ≥ 2.0. The p value of the F-test statistic for 02 as 0.0, which implies a significant difference in the _ value means. The p value of the F-test statistic for 03 was 0.14, which confirms that does not influence the _ . In the present case, the interaction effect is insignificant. The p value of the F-test statistic for the effect of the interaction of the inducing factor (AWR, mining and geothermal energy production) on was 0.65.
To obtain more information about where the differences between the groups lie, we ran post hoc tests. The Tukey honest significant difference (HSD) test is used to determine if the interaction between two sets of data is statistically significant. The value of the Tukey test is obtained by dividing the absolute value of the difference between the pairs of means by the standard error of the mean. The post hoc Tukey's test confirmed that the _ values of the S-SE cycle earthquakes in the ST2 region are significantly lower than the _ values calculated for the seismic events in LGCD and USCB, with p value of less than 0.0002.
Inversion of the peak ground acceleration and the peak ground velocity from the rupture direction and rupture velocity. A simple inversion of the peak ground acceleration (PGA) and peak ground velocity (PGV) from the rupture direction and rupture velocity was calculated based on the methodology developed by Boatwright (81). This method is based on the directivity function derived from unilateral ruptures, which was developed by Ben-Menahem (112). The residuals from each event are fitted to the logarithm of the Ben-Menahem directivity function, which depends on the rapture direction and rapture velocity. To calculate the inversion, in the first step, the appropriate ground motion parameters (PGA and PGV) were determined based on seismograms of selected events, and then, the Ground Motion Prediction Equations (GMPE) were constructed.
Data. The ground motion parameters were obtained based on the event-related waveforms produced by 147 events. Due to the lack of information about the azimuth and S-wave propagation direction, one event was excluded from the catalogue. The PGV was calculated, and the PGA was estimated as the total ground motion peak based on the measured seismograms from every station. All the measured seismograms and the accelerograms obtained via numerical differentiation were examined through visual inspection. Therefore, to ensure the quality of the information, the signals of two events were excluded. Furthermore, the parameters from station TNVB were also excluded from the analysis because this station was the largest distance from the events epicentres (largest minimum distance of 20 km) (Fig. 1). Finally, the ground motion catalogue (GM catalogue) contained 934 records of ground motion parameters produced by 145 events measured at particular stations. The range of PGV was 0.0001 to 0.66 cm/s, and that of PGA was 0.00004 to 0.94 m/s 2 . The majority of the ground motion parameters were weak, i.e., PGA of less than 0.03 m/s 2 and PGV of less than 0.02 cm/s, accounting for 85% and 83% of the data, respectively. A significant part of these data is at the noise level (PGA of <0.001 m/s 2 and PGV of <0.0012 cm/s), that is, around 24%. However, all these data were used to invert PGA and PGV since the visual inspection of the data confirmed their good quality. This ensures the appropriate number of recordings of individual events, which is essential for the inversion of these data.

GMPE.
To perform the inversion, we needed the residual values of the peak motions. These residuals are the difference between the observed and estimated peak motions calculated using the GMPE. A completely new model of the GMPE was estimated because there was no existing GMPE model for the ST2 area. The following model of the GMPE for PGA and PGV was established: log ( , ) = α + β − log √ 2 + ℎ 2 + w , where is the station's ordinal number; is a ground motion parameter at the ℎ station (PGA, PGV); is the moment magnitude; is the source-station epicentre distance; ℎ is the common depth factor introduced to account for the non-linearity of the ground vibration propagation at short distances from the epicentre (ℎ is estimated to ensure the minimum standard error of the estimation); is a free value, which depends on the station's location (relative amplification factor at station relative to station ) and can also be interpreted as the logarithm of the relative amplification. , , , and are all parameters estimated via regression analysis, and they give a set of GMPEs dedicated to every station with common , , coefficients. The model considers an increase in the GM parameters with increasing source size, non-linear geometric scattering, and the influence of the site effect at the station. This complex model is based on the solutions proposed by Joyner & Boore (113), Field (114), and Lasocki (115). Amplification factors are usually included in GMPEs as site categories. To prepare the residuals according to the method of Boatwright (81), corrections to the stations needed to be made. Therefore, we chose the model dealing with the amplification factor relative to the station position (115). The GMPEs were estimated using the GMPE application of the IS-EPOS platform (116). The coefficients of the GMPE were estimated using the ordinary one-stage regression method in three steps, including the statistical and physical consistency of the model, the estimation of the fixed depth, and analysis of the residuals (the entire routine is explained in the IS-EPOS platform User Guide, 117). Table 1 presents the coefficients of the final models. Inversion. We estimated the peak values, event correction, and residuals of the peak values (PGA and PGD). These parameters are included in the GM Catalogue. The rupture velocity and direction were estimated using the inversion method proposed by Boatwright (81). The inversion results for 143 events were calculated based on the required parameters: the shear wave velocity at the source ( = 3.42 / at depths of 0-17 km), the propagation direction of the shear wave at each station, the azimuth of each station, the PGA and PGV at each station, and the residual peak motions. For some cases, the inversion mechanism did not give a convergent result, so the iteration stopped in the initial stage and the result was close to zero. Therefore, we tested the stability and correctness of the results for each event. The inversion of the peak ground values of each event was calculated using the jack-knife resampling technique, i.e., systematically leaving out one observation at a particular station. Next, the outliers of the inversion results were evaluated for each event. Some of the incorrect inversion results were improved by replacing them with inversion results obtained using resampling and by omitting data from one station. Additionally, the standard deviation of the correct inversion results without outliers was calculated. The stability and correctness of the results were marked by the attributes using the following rules. (i) Attribute 1: stable and correct inversion results, that is, the results obtained for all the stations were correct without outliers and the standard deviation was less than 0.1. (ii) Attribute 11: correct inversion results, that is, the results computed for all the stations were incorrect, but the results were replaced by the outcome obtained by omitting one of the stations, and the standard deviation was less than 0.1. (iii) Attribute 2: unstable inversion results, that is, correct inversion results were obtained for all the stations or for the selected stations, but the standard deviation was greater than 0.1. (iv) Attribute 0: there were not any correct results even for the selected stations.

Table 1 Final model of the GMPE for the peak ground motion parameters
/ model. The / model was investigated using the local P-wave tomography. The earthquake catalogue and 1D velocity model (Fig. 8A) were downloaded from the IS-

EPOS platform of the European Plate Observing System's (EPOS) Thematic Core Service
Anthropogenic Hazards (https://tcs.ah-epos.eu). This catalogue contains the parameters of more than 7000 events that occurred from August 24, 2013, to June 2, 2017. We selected 4433 events based on the following strict data selection criteria. (i) The event location uncertainties (both horizontal and vertical) were less than 500 m. (ii) The RMS value was less than 0.2 s. (iii) The azimuthal gap was less than 180°. (iv) At least 6 Pg arrivals per event were used. Then, a linearized, iterative, damped least-squares computer program (SIMULPS) (118) was applied to solve the non-linear / structure problem. Using the initial 1D velocity model, we parameterized the study area with an even geographic grid spacing of 1 km (in both the north-south and east-west directions), and vertically, with multiple layers with a constant thickness (1 km) at depths of 0-5 km. Using the final RMS, which should reach a minimum value, the optimum number of inversions and regularization parameters (e.g., damping, etc.) for the iteration were set based on the standard L-curve (Fig. 8B) (119). The quality of the results provide us with the areas with fairly good resolution in any seismic tomography study. Thus, the checkerboard resolution was tested to assess the quality of the inversion results. Therefore, the 1D initial velocity model was alternately perturbed with a maximum error of 10% at similar actual grid points. The horizontal slices of the tomographic results (left maps in Fig. 9) and the checkerboard resolution tests (right maps in Fig. 9) are shown in Fig. 9.  Competing interests: The Authors declare that they have no competing interests.
Data availability: All data used in the analyses are available on the IS-EPOS platform of Thematic Core Service Anthropogenic Hazards: https://tcs.ahepos.eu/#episode:SONG_TRANH. All data needed to evaluate the conclusions in the paper are presented in the main text of the paper and Methods. Additional data related to this paper may be requested from the authors.   (43,44). The timing of the stiffness drops is shown in light orange. The dashed yellow and blue rectangles mark the two different parameter signatures of the process indicating potential different mechanisms responsible for the stiffness drop. The first parameter leads to the porosity loss resulting in the pore pressure increase (compression, shown by yellow rectangle), and the second leads to the volume expansion resulting in the pore pressure decrease (dilatation, shown by blue rectangle). All the recorded seismicity is marked by the grey circles.     7 | Model of the interaction between the slow slip and seismic events triggered by the water reservoir impoundment. VW, VN, and VS are the velocity weakening, velocity neutral, and velocity strengthening regions, respectively, and they are represented by the light yellow to dark brown colours in the seismogenic zone. The contacting rupturing interfaces may possess different elastic properties, which promotes the pulse-like mode of the rupture (e.g., 79). The observed rupture directions are up-and downdip.

Fig. 8 | (A) 1D (red) and (blue) models used as the input for the 3D inversion. (B)
Standard Lcurve: model roughness versus data misfit, which was used to obtain the optimal damping regularization parameter. The black circle is the optimal damping value of 20.   ST2. The inversions were performed using the instability criterion in a moving window of 15 earthquakes with a 1 event-moving step. The dotted lines indicate the 95% confidence intervals. C) Time variations of the shear stress along the fault plane (orange line) based on the relative stress magnitude analysis; the stiffness changes (brown line) based on the stress drop and slip of the earthquakes; and the corresponding friction rate parameter (green line) calculated assuming a critical slip distance of 10 −2 (43,44). The timing of the stiffness drops is shown in light orange. The dashed yellow and blue rectangles mark the two different parameter signatures of the process indicating potential different mechanisms responsible for the stiffness drop. The first parameter leads to the porosity loss resulting in the pore pressure increase (compression, shown by yellow rectangle), and the second leads to the volume expansion resulting in the pore pressure decrease (dilatation, shown by blue rectangle).   Fig. 2. B) Slip weakening curve, with the area below the curve indicating the minimum breakdown work of the slow slip events (S-SE) over the characteristic distance − . Here, − is considered to be a proxy of the critical slip weakening distance of the S-SE. C) Stacking of the stiffness evolution in all of the identified seismic cycles in the analysed period. Colours indicate groups with the same − . D) Proxy of the critical slip weakening distance − vs. the pore fluid pressure factor (dots) and the maximum observed magnitude (stars) in the weakening cycle. The colours correspond to the groups in c, the unfilled symbols correspond to the particular weakening cycles and the larger filled symbols correspond to the mean and maximum values of and , respectively, for the identified groups in c. The dashed and solid grey lines are the logarithmic functions relating − to and , respectively. The approximated rise times and the 95% confidence intervals (blue circles) are shown, and they are related to the data from the mining induced seismicity in the Upper Silesian Coal Basin (red circle), the Legnica-Głogów Copper District in Poland (brown circle) and the geothermal energy production seismicity in the Geysers geothermal field in the USA (green circle). The large blue circles correspond to the events with a magnitude of greater than 3.0, and the colours linked to the depth and the size of the magnitude are the same as in Fig. 2.   Fig. 6. Scaling of fracture energy with slip attributable to thermal pressurization. The fracture energy was estimated for the seismic events (blue empty circles) and slow slip events (blue solid cycles) in the Song Trangh 2 region and earthquakes in the continental crust and subduction systems (grey empty circles for small continental events and grey solid cycles for large continental and subduction events) (compiled by 94). For the Cajon Pass data, the USCB, LGCD, and TG data, see Fig. 5).
Regarding the seismic events in the ST2 region, the fracture energies were approximated from the energetic quantity ′ (95), which was calculated from the seismic source parameters, , the static stress drop, ∆ , and the apparent stress, : ′ = (