Investigating reservoir-triggered seismicity in the Koyna–Warna region through nodal ambient noise tomography

Abstract

The Koyna–Warna region in southern India is a unique case of dual Reservoir-Triggered Seismicity (RTS), surrounding impoundments on Koyna and Warna rivers. In this study, we perform ambient noise surface wave tomography using data from a temporary dense deployment of three-component geophones in order to understand the origin of RTS. To achieve this objective, surface wave dispersion curves from noise cross-correlations in rotated coordinates were used to create group velocity maps in 0.5–5 s periods. The group velocity dispersions from map bins were inverted for shear wave velocity () and radial anisotropy () maps down to 5 km depth. The and till 1 km depth are correlated with topographical features in agreement with the basement depth estimates. Large earthquakes in the region are marked by negative with negative anomaly or anomaly contrast, suggesting shear strength could be an important control on their genesis. A large negative anomaly southwest of the Koyna reservoir is linked with reservoir water migration to greater depths, part of a larger interconnected system of pathways extending across the reservoirs. These results, set decades post reservoir impoundment, can help understand the evolution of RTS around dam sites elsewhere.

Introduction

The RTS has been evidenced in different parts of the world such as Lake Mead (USA), Xinfengjiang reservoir (China), Kariba reservoir (Zambia–Zimbabwe Border), Kremasta (Greece), Koyna (India) etc.¹. Among these, the Koyna–Warna region is characterized by one of the largest RTS event ever recorded ( 6.3) around the world. In addition, unlike the remaining RTS sites, the Koyna–Warna is a multi-reservoir system, which provides a unique opportunity for understanding the impact of multiple impoundments on the aquifers in a broad area. The insights from the seismic velocity structure in such a complex RTS region would help scientists and engineers in the planning of future dam sites globally.

On 10th December 1967, an 6.3 earthquake occurred south-west of the Koyna Dam, which is located on the Koyna river in Maharashtra, India (Fig. 1a,b). It caused widespread damage in the nearby town of Koynanagar, with houses completely destroyed and over 200 people killed². Such a large earthquake was unprecedented in this region with no record of seismicity prior to the start of the dam construction in 1956. After the associated reservoir was filled in 1962, seismicity started in the Koyna region, which could be recorded at the seismological station in Pune, 120 km away³. The minor seismicity culminated in the occurrence of the Koyna earthquakes of 1967^{4, 5–6}. Later when the neighbor Warna reservoir was filled after the creation of Warna dam in 1985, seismicity started to appear around Warna reservoir as well. Over the last 60 years, more than 20 earthquakes with M > 5 and 200 earthquakes with M > 4 earthquakes have occured in the Koyna–Warna region⁷.

[See PDF for image]

Fig. 1

Maps of the Koyna–Warna region, showing topography, shallow large earthquakes (M > 5, depth 5 km; red asterisk), seismicity color scaled by depth and compiled from various sources^{8, 9–10}, and the stations used in this study (green triangles). In (a) broad geographical features are labeled, while (b) is a zoomed-in map (area under blue rectangle in (a)) highlighting smaller features. Locations of nine cored boreholes are shown in grey solid circles (KBH1 borehole marked). The pilot deep borehole KFD1 is shown using black solid star. KRFZ stands for Koyna River Fault Zone. Profiles A to E along which results are analyzed are also shown in (b) in brown. The maps were generated using Python-v3.8.18 scripts with pygmt-v0.9¹¹ and combined using Inkscape-v1.2.2¹², please see Data Availability section for more details.

The origin of the seismicity pattern in the Koyna–Warna region is still debated, however some understanding has been established. The emergence of minor seismicity around the time of the first impoundment of Koyna and Warna dams was later investigated, and a high correlation between water level changes in the reservoirs and seismicity rate was found^13,14. However, water level changes are not the only factors determining the seismicity pattern, as stress transfer between local faults can explain some of the focal mechanisms, suggesting that the region is divided into multiple tectonic blocks^{15, 16–17}. Currently, the consensus is that the seismicity in this region is “triggered” along a network of local faults by the diffusion of pore pressure, due to reservoir water filling cycle^1,8,14,18

The Koyna–Warna region is part of the Western ghats (stepped terrain) in the Deccan Volcanic Province (DVP), which formed due to the passage of the Indian plate over the Reunion hotspot about 65 Ma ago¹⁹. In this region, the DVP is composed of nearly horizontal simple tholeiitic lava flows with thickness up to 1.2 km, underlain by a granitic-granitoid basement²⁰. As part of the scientific deep drilling project, nine cored boreholes were drilled in the Koyna–Warna region (Fig. 1a,b), which shed new light on the stratigraphy²¹. The basement granite-granitic gneiss is inter-layered with migmatitic gneiss, in addition, evidences of fluid-rock interaction that has led to mineralization along micro-fractures has been found^22,23. Using results from these cored boreholes, a deep research borehole was drilled recently down to 3 km depth with the purpose of well logging²⁴.

The Koyna–Warna region hosts multiple NE-SW faults along with NW-SE structural lineaments. The 1967 Koyna earthquake (M 6.3) occurred along a NNE-SSW fault and generated surficial fissures up to 4 km SSE of Donichawadi (Fig. 1b). Focal mechanism solutions of the mainshock and the aftershocks revealed a left-lateral strike-slip fault extending SSE of the Koyna reservoir^4,17,25. Aeromagnetic maps also indicate a NNE-SSE magnetic lineament side-by-side this fault zone²⁶. This NNE-SSW lineament along with the Donichawadi fault is referred to as the Koyna River Fault Zone (KRFZ) or Koyna Seismogenic Zone (KSZ) in literature^16,27. By combining aeromagnetic, SAR interferometry and bouger gravity anomaly maps, the NNE-SSW lineament has been found to be a shallow feature and seems to be connected to deeper NW-SE lineaments²⁸. The majority of the seismicity near the Koyna reservoir originates in 7–13 km depth range of mainly strike-slip mechanisms, in contrast to the Warna region where majority of the seismicity is shallow (< 7 km) with normal faulting mechanisms^8,9,17,29.

The subsurface structure in the Koyna–Warna region has been imaged using multiple geophysical tools such as seismic, gravity, magnetotelluric methods etc. Deep seismic sounding along two roughly E–W profiles bounding the Koyna reservoir in the north and the south, revealed a low-velocity layer between 6 and 11.5 km depth³⁰. The 2D gravity models indicate a gravity low at Koyna, which is attributed to isostatic adjustment that leads to thickened crust³¹, or an additional low density layer at 5–15 km depth³². The 2D magnetotelluric image along a WNW-ESE profile across the Koyna Seismogenic Zone reveals a shallow (< 10 km) conductive anomaly interspersed between very high resistive blocks that are likely Deccan basalts³³. The first detailed P and S velocity maps of the Warna region was performed using 343 local earthquakes around the Warna reservoir³⁴. Zones of high / and low were observed beneath the Warna reservoir at less than 1.5 km depth and extending down to 7 km depth in lateral sections.

Despite decades of research in the Koyna–Warna region, the linkage between the genesis of large earthquakes and the shallow structure is not well understood, partly because it has only been studied in a small area of the Koyna–Warna region. Due to non-uniform distribution of earthquakes, past travel time tomography has yielded resolution limited to the Warna reservoir^34,35. Ambient noise tomography is based on measuring phase or group velocity changes of the surface wave arrivals that emerge after cross-correlating ambient noise waveforms from pairs of stations³⁶. Because ambient noise sources are abundant around the globe, this technique allows better spatial resolution in regions with non-uniform or sparse earthquake distribution³⁷. It has also been shown that at shorter periods (< 10 s), stacking of 1 month of cross-correlations provides high enough SNR for ambient noise tomography and the required number of stacked days reduces with decreasing period^38,39. Ambient noise cross-correlations of vertical–vertical and transverse–transverse components are used to compute Rayleigh wave and Love wave group velocities respectively, from which the shear wave vertical velocity () and shear wave horizontal velocity () as well as shear wave radial anisotropy () are derived. Conventionally, is negative for vertically aligned structures such as dykes, faults etc. and positive for horizontally aligned structures such as lava flows, sedimentary layering etc.^{40, 41–42}. So far ambient noise cross-correlations have only been used to derive 1D shear wave velocity model for the whole Koyna-Warna region⁴³.

In this study, we use waveform data from ninety seven 4.5 Hz three-component geophone nodes installed in the Koyna–Warna region between January 2010 and May 2010 in order to compute cross-correlations of seismic noise waveforms recorded on the vertical as well as transverse components. We measure the the fundamental mode Rayleigh wave group velocity dispersion curves from the stacked vertical–vertical (ZZ) cross-correlations and the fundamental mode Love wave dispersion curves from the stacked transverse-transverse (TT) cross-correlations using all station pairs. After that, we invert the Rayleigh and Love-wave group velocity dispersion measurements to obtain 2-D group velocity maps. Finally, by separately inverting the group velocity maps from ZZ and TT cross-correlations, the shear wave radial anisotropy and the isotropic shear wave velocity maps are also generated.

Results

The isotropic shear wave velocity () in the upper 1 km of the Koyna–Wanra region indicates strong correlation with the topographical features (Fig. 2a,b). In 0 km depth lateral section, towards the west of the Western ghats escarpment, negative anomalies are observed that seem to follow the Western ghats escarpment. Similarly, negative anomaly beneath the Warna reservoir and the southern part of the Koyna reservoir is observed. In addition, moderately negative anomalies are also observed beneath the Morana reservoir. On the other hand, positive anomalies are observed in the eastern part of the Koyna-Warna region. In 1 km depth lateral section, similar correlation between anomalies and the topographical features are observed, although the low beneath the Warna reservoir shifts towards the Warna dam in the east and the Western ghats escarpment in the west (Fig. S1).

[See PDF for image]

Fig. 2

Isotropic shear wave velocity anomalies in lateral sections (a–d) from surface down to 5 km depth. The white contours are plotted at every 2.5 % anomaly. The outline of the Koyna, Warna, and Morana rivers as well as respective associated reservoirs are shown (black). The western ghats escarpment is shown with a ridged line. The stations are shown as red triangles in (a). The seismicity within 0.5 km depth of each lateral section is shown by green solid circles. Focal mechanisms of large earthquakes ( > 5)^{4, 5–6,8} within 0.5 km depth of the lateral sections are also shown using beach balls. The locations of the cored borehole KBH1 and the pilot deep borehole KFD1 are shown with hollow circle and star respectively and the dashed black line marks the Donichawadi fissure same as in Fig. 1. The maps in this figure were generated using Python-v3.8.18 scripts with MSNoise-v1.3 and MSNoise-Tomo⁴⁴, swprepost-v2.0.0⁴⁵, Geopsy-v3.4.0⁴⁶, and pygmt-v0.9¹¹, and combined using Inkscape-v1.2.2¹², please see Data Availability section for more details.

From 2 km depth onwards, a negative anomaly is observed beneath the central part of the Western ghats escarpment down to 5 km depth. One of the earthquakes in the 1967 Koyna earthquake sequence also occurred in this region (Fig. 2c). Apart from that, negative anomaly is only observed directly beneath the north-western part of Warna reservoir down to 2 km depth. At greater depths, moderately positive anomalies are observed directly beneath the reservoir, however these are surrounded by negative anomalies. The majority of the seismicity in 0–5 km depth also occurs outside the boundaries of the Warna reservoir, mainly towards the south-west.

[See PDF for image]

Fig. 3

Shear wave radial anisotropy () in lateral sections (a–d) from surface down to 5 km depth. See Fig. 2 caption for more details. The maps in this figure were generated using Python-v3.8.18 scripts with MSNoise-v1.3 and MSNoise-Tomo⁴⁴, swprepost-v2.0.0⁴⁵, Geopsy-v3.4.0⁴⁶, and pygmt-v0.9¹¹, and combined using Inkscape-v1.2.2¹², please see Data Availability section for more details.

The radial anisotropy () is mostly positive down to 1 km depth (Fig. 3a, Fig. S1), which means at these depths is higher than . At greater depths however, a negative is observed overall. At depths 2 km, negative is localized around the 1967 Koyna earthquake epicenters, the northern part of Warna reservoir, and between Warna and Morana reservoirs.

The cross-sections in the central part of the Koyna–Warna region reveal evolution of seismicity in the region across several and anomaly zones (Fig. 4a–g). The depth resolution is not uniform across the cross-sections. It varies between 5 and 6.5 km after topographic adjustment depending on the wavelength content of the dispersion data (see Methods section).

estimated locally around the Warna reservoir using local earthquake tomography³⁴ agree well with the results in the present study. At 1.5 km depth, low beneath the northern part of the Warna reservoir, found previously, correlates with the negative anomaly at 2 km depth (Fig. 2b) and 1.5 km depth (Fig. S1) in the present study. In addition, it was observed that most of the of the seismicity occurred in low velocity zones that were present outside the reservoir boundaries³⁴, similar to present study. Furthermore, they also observed that the low zones close to the Warna reservoir migrated northwards at greater depths. This pattern could be related to the negative anomaly that we observe in the central part of the western ghats escarpment in lateral sections at 2 km depth and deeper (Fig. 2b–d).

The shallow structure in the Koyna–Warna region is dominated by Deccan basalts, which are found to be 800–1200 m thick across the region based on borehole data²² and magnetotelluric experiments²⁷. These formed as a result of large scale melting when the Indian plate passed over the Reunion hotspot¹⁹. Most of the lava flow in this region is characterized by rubbly pahoehoe type that consists of a thick sheet like smooth flow base up and a thin rubbly top with slabs of lava embedded in it⁴⁷. This sheet like geometry of the Deccan lava flow can explain the positive observed in the top 1 km depth (Fig. 3a, Fig. S1)

[See PDF for image]

Fig. 4

Isotropic shear wave velocity anomalies (a,c,e) and shear wave radial anisotropy (b,d,f) in lateral sections from surface down to 5 km depth and the outline of cross sections in a spatial map (g). radial aniso. In (b,d,f) refers to shear wave radial anisotropy. The Koyna reservoir location is denoted with “K” and the warna reservoir locations are denoted with “W” on the top of the cross-sections. The seismicity in the area in 1967–2017 period^{8, 9–10} within 1.25 km of the cross-sections has been plotted on top and color coded by their date of occurrence. The contour in (a) highlights linked negative anomalies from surface down to greater depths connecting the two reservoirs. The focal mechanisms have been projected based on the orientations of the vertical planes. The maps in this figure were generated using Python-v3.8.18 scripts with MSNoise-v1.3 and MSNoise-Tomo⁴⁴, swprepost-v2.0.0⁴⁵, Geopsy-v3.4.0⁴⁶, and pygmt-v0.9¹¹, and combined using Inkscape-v1.2.2¹².

Large earthquakes in Koyna–Warna and the strain budget

Large earthquakes globally exhibit two major patterns in terms of seismic velocities of the source zones: either they occur at the transition between high and low velocity lithology or they occur in a low velocity zone indicating a fractured medium allowing migration of fluids from a deeper source^48,49. The occurrence of these large earthquakes mostly control the strain budget in a region.

Two main theories have been proposed to explain the strain budget in the Koyna–Warna region. First theory states that the geometrical orientations of faults allow stress transfer such that the release of stress due to a large earthquake on one of the faults destabilizes another fault. The stress triggering due to earthquakes in the Koyna region generates earthquakes in the Warna region and has resulted in continued seismicity until now¹⁶. This theory may explain the occurrence of a normal faulting earthquake ( 5.6) following the strike-slip earthquake ( 6.3) a few days apart, in the 1967 sequence^{4, 5–6}. The second theory advocates for alternating phases of strike-slip and normal faulting earthquake predominance in the region due to periodic changes in the regional stress pattern⁸. However, most of the large earthquakes following the 1967 sequence have been normal faulting type with the fewer strike slip earthquakes also having a significant normal component. Therefore, stress transfer due to geometric orientation of faults as well as periodic variation in stress levels may not fully explain the seismicity pattern.

The source zones of the first two earthquakes in the 1967 sequence are characterized by transition from mostly positive to slightly negative anomaly (Fig. 4a), while highly negative (Fig. 4b). The highly negative is correlated with the predominantly strike-slip mechanism of the two earthquakes on steeply dipping fault planes. Steeply dipping fault planes result in lower compared to , hence negative . The anomaly pattern suggests that the earthquakes occurred at a structural boundary. The third earthquake in the sequence is characterized by a negative anomaly and slightly positive to negligible (Fig. 4c,d). The negative anomaly indicates a weakened crust before or after the earthquake. A similar pattern can be observed near the Warna reservoir, where some of the earliest recorded large earthquakes happened in a neutral zone (1974 event; Fig. 2d) and at a structural boundary marked by transition between high and low anomaly (1976 event; Fig. 2c). The source zones of both events lie in a negative zone (Fig. 3c,d). The subsequent large earthquakes occurred in negative anomaly zones with moderately negative (Figs. 3c,d, 4e,f). An exception is the 1975, Mw 5.0 event that occurred in a negative zone with highly positive at 0.1 km depth (Figs. 2a, 3a). The unusual depth of this event i.e. above the basement, along with its position close to the western ghats escarpment, away from the reservoirs, suggests either poor event location or unknown extremely shallow fault activity.

Theoretical models of earthquake nucleation on both slip-weakening and rate-and-state faults indicate that critical nucleation length after which the quasi-static slip develops into a dynamic rupture is proportional to the shear wave modulus or shear strength, and it is particularly significant for slip-weakening faults^50,51. In addition, numerical experiments on earthquake nucleation in presence of a precursory drop in shear wave velocity suggests that the velocity anomaly can lead to reduction in nucleation size as well as earlier nucleation. Thus, such precursory velocity drop can lead to earlier onset of large earthquakes⁵².

The distribution of smaller earthquakes across the cross-sections suggests a net southward migration of seismicity from Koyna reservoir towards Warna reservoir in the last 60 years (Fig. 4). The recent seismicity (year > 2000) at shallow depth (< 5–6 km) is concentrated around the negative anomaly zones particularly in the central and southern part of Koyna-Warna region (Fig. 4a,e,f). Measurements of geodetic strain rates indicate low annual tectonic deformation of 0.2–0.3 0.04 mm/year in the region, much lower than that of Indian plate⁵³. Even with such low deformation rates, seismicity has continued in this region till date. In recent years, it has concentrated in the negative anomaly zones likely due to their low shear strength. This also means that the accumulated stress is more frequently released in the south-central part compared to the northern part that is close to the Koyna reservoir. Therefore, there is a probability of higher stress accumulation and hence a greater chance that the next large earthquake may occur close to the Koyna reservoir. On the other hand, seismicity near the Warna reservoir would likely continue but become smaller with further decrease in shear strength.

Influence of reservoir impoundment on the structure

The immediate effect (< few years) of reservoir impoundment on the structure has been observed in the form of compaction beneath the reservoir, leading to higher velocities⁵⁴. However, over time, water infiltration to greater depths has also been evidenced by increased seismicity levels along with lower shear wave velocity⁵⁵. In this section, we study after-effects of reservoir impoundment on recent shear wave velocity structure of the Koyna–Warna region.

A large negative anomaly is observed in the spatial maps as well as the cross-sections SW of the Koyna reservoir and extending westwards across the western ghats escarpment (Figs. 2b–d, 4a,g). This anomaly continues from 2 km depth down to the 5 km depth and possibly beyond (Fig. 4a) and it is marked by moderately negative to negligible (Figs. 3b–d, 4b). The isotropic nature of this velocity anomaly suggests that the structure underneath exhibits multi-directional deformation. At shallow depths, this anomaly seems to be linked to the Koyna reservoir (Fig. 4a; “K”). Fluid migration pathways from shallow levels down to greater depths along weakened zones may explain the lateral and depth continuity of this anomaly.

Increased pore pressure diffusion to greater depths due to annual filling and emptying cycle of the Koyna and Warna reservoirs, and release after the occurrence of a large earthquake, has been identified as a possible mechanism of seismicity in the region^1,14,18. Well logs in a 3 km pilot research borehole (KFD1) at a nearby site, indicate higher water content in 1.5–3 km depth that suggest deeper penetration of reservoir water⁵⁶. Besides, mineralogical and geochemical analyses of the cores from KBH1 borehole near the pilot borehole have identified precipitation of ferruginous and siliceous minerals along fractures at different depth levels (800–1200 m), which indicate percolation of water to greater depths^23,57. Evidence of fluid migration SW of the Koyna reservoir has also been found in the form of a westward dipping conductive zone in an E-W MT profile passing south of the 1967 Koyna earthquake ( 6.3) epicenter²⁷. Shallow parts of the dipping conductive structure roughly coincides with the negative anomaly zone especially around 4–5 km depth. Therefore, the conductive zone could be a downward continuation of the negative anomaly.

The spatio-temporal distribution of seismicity reveals that the earliest earthquakes in the region occurred very close to the Koyna reservoir (Fig. 4a) likely due to reservoir loading and subsurface compaction. The smaller fractures developed due to this seismicity may have facilitated reservoir water migration. The left lateral slip along with a small normal faulting component during the 1967 Koyna earthquake ( 6.3) might have weakened the structure SW of the Koyna reservoir that allowed reservoir water to migrate downwards and southwestwards over several decades aided by smaller earthquakes (Fig. 4e,f). The downward migrating water under higher pressure and temperature would have induced chemical alterations in the structure which is evidenced as negative anomaly.

Tectonic link between Koyna and Warna reservoirs

The anomalies provide insights into how the RTS in the Koyna–Warna region may have evolved. The first two 1967 Koyna earthquakes ( 5.9 and 6.3) produced left-lateral strike slip movement leading to the transfer of stress southwards and westwards (Fig. 4a). The regions surrounding the two epicenters were weakened due to the reservoir water migration, evidenced by negative anomalies with low degree of anisotropy (Fig. 4a,b). The stress transfer may have pushed the already weakened zones beyond the critical limit, allowing small fractures to develop into large earthquakes⁵⁸, the first of which occurred 2 days after the 1967 main shock with a normal faulting mechanism (Figs. 2c, 3c, 4c,d). Over time the stress may have transferred further southwards, leading to the 1974–1976 Warna earthquakes in a region of relative higher shear strength, evidenced by moderately positive anomaly (Figs. 2c, 3c). After the Warna reservoir got filled in late 80s–early 90s, the reservoir water migration likely started occurring near Warna reservoir, similar to the Koyna reservoir case. This phenomena produced weakened shear strength in the surrounding regions (Figs. 4e, 5) and coupled with stress transfer, likely generated the 2009 Warna earthquake sequence with normal faulting earthquakes (Figs. 4e,f, 5). Thus, Koyna–Warna region may represent an interconnected aquifer system, where reservoir water migrates from shallow to deeper levels as well as laterally, generating zones of weakened shear strength (Fig. 5).

[See PDF for image]

Fig. 5

Isotropic shear wave velocity anomalies in two cross-sections intersecting koyna, warna, and morana reservoirs, and a spatial map showing the outline of cross sections. The locations of cross-section profiles are shown in a lateral map in Fig. 1b. The maps in this figure were generated using Python-v3.8.18 scripts with MSNoise-v1.3 and MSNoise-Tomo⁴⁴, swprepost-v2.0.0⁴⁵, Geopsy-v3.4.0⁴⁶, and pygmt-v0.9¹¹, and combined using Inkscape-v1.2.2¹², please see Data Availability section for more details.

Conclusion

We performed ambient noise tomography of the Koyna–Warna region using data from 97 three-component geophone nodes installed for a 5 month period between Jan 2010 and May 2010. Examining the spectrogram and power spectral density of the waveforms as well as manual picking of the group velocity dispersion curves, provided reliable measurements between 0.5 and 5 s period. Inversion of the dispersion curves in this period range allowed estimation of 1D model down to 5 km depth across the region. The 1D models were interpolated to generate 2D lateral and cross-section maps of as well as . The anomalies correlate with topographical features down to 1 km depth and the basement depth below the Deccan basalts is identified using values, in agreement with the previous results. The source zones of large earthquakes in the region exhibit anomaly contrast or negative anomaly, consistent with global occurrences as well as theoretical and numerical models, suggesting structural shear strength changes as an important factor in assessing potential for future large earthquakes. Negative anomaly SW of the Koyna reservoir close to the location of KFD1 and KBH1 boreholes provides seismic evidence of reservoir water percolation to 5 km depth and beyond. Extended anomalies in depth and laterally, indicate an interconnected system of water migration pathways between Koyna, Warna, and Morana reservoirs, and their associated rivers. Further studies on the temporal variation of can enhance our understanding of the evolution of these water migration pathways and their implications for seismic hazard.

Methods

Preprocessing and computation of cross-correlations

A network of ninety-seven 3-component geophones (Geospace GS-11D) of 4.5 Hz natural frequency coupled with Taurus dataloggers was deployed in the Koyna–Warna region from January 2010 to May 2010 to record the local seismicity³⁴. Spectrogram of waveforms indicate high noise levels below 0.1 s period and above 2 s period, however variations in the noise are detectable at shorter periods as well as longer periods (Fig. 6a). The power spectral density of the waveforms shows that the noise amplitudes lie between the global noise models up to 7 s period (0.14 Hz) and down to 0.1 s period or shorter (> 10 Hz) (Fig. 6b).

[See PDF for image]

Fig. 6

(a) Spectrogram for waveforms from all stations in the network. The spectrogram amplitudes have been converted to decibels. (b) Power spectral density of waveforms from all the stations. Grey curves indicate the global New High Noise Model (NHNM) and New Low Noise Models (NLNM)⁵⁹. (c) Heat map of all the Rayleigh wave dispersion curves using ZZ cross-correlations. (d) Heatmap of all the Love wave dispersion curves using TT cross-correlations. The plots in this figure were generated using Python-v3.8.18 scripts with Obspy-v1.4.0⁶⁰ and Matplotlib-v3.7.3⁶¹, and combined using Inkscape-v1.2.2¹², please see Data Availability section for more details.

Considering this analysis, the waveforms were deconvolved using their instrument responses and then bandpass filtered between 0.16 and 12 Hz. Next, the waveforms were decimated from 100 Hz sampling frequency to 40 Hz. After that, the daily waveforms were divided into 2 h time windows and spectrally whitened to flatten spurious events in the waveforms. Finally, the waveforms were cross-correlated with max lag of 1 h and all cross-correlations were linearly stacked for the entire duration of the installation. Using this process, cross-component and same component cross-correlations: North-North (NN), East-East (EE), East-North (EN), North-East(NE), and Vertical-Vertical(ZZ) were computed. The waveform preprocessing and computation of cross-correlations was performed using MSNoise version-1.3⁴⁴. All the cross-correlations of horizontal components were rotated to radial and transverse coordinates: Radial-Radial(RR), Transverse-Transverse (TT), Radial-Transverse (RT), and Transverse-Radial (TR)⁶². In total, 4050 ZZ cross-correlations and 3989 TT as well as RR cross-correlations were obtained.

Computation of dispersion curves

Group velocity dispersion curves were computed using Frequency-Time ANalysis (FTAN) of the cross-correlations that involves simultaneously analyzing the cross-correlations in both frequency and time domains⁶³. The Rayleigh wave dispersion curves were computed using the ZZ cross-correlations, while the Love wave dispersion curves were computed using TT cross-correlations. After an automatic computation of dispersion curves, these were manually re-picked one by one in order to remove discrepancies and the FTAN normalized amplitudes of the picks were recorded to determine the quality of picks (see Fig. S2 for FTAN illustration). At short distances, the dispersion measurements of long periods are not reliable. Therefore, measurements were only made for periods corresponding to inter-station distances that are larger than one wavelength. At this stage, 2912 dispersion curves from TT cross-correlations and 3066 dispersion curves from ZZ cross-correlations were obtained.

By simultaneously plotting dispersion curves from all the cross-correlations in the form of heat map, their average trend and the scatter in the measurements can be observed (Fig. 6c,d). In general, Love wave group velocities are higher than Rayleigh wave group velocities. For periods less than 1 s, larger scatter in the dispersion curves is observed, which could be due to the incoherent noise, possibly human-induced, interfering with the coherent noise sources⁶⁴. On the other hand, measurements at periods longer than 3 s also show large scatter, however this is likely due to fewer measurements made at those periods. In order to purge measurements that may potentially be erroneous, we only consider period range of 0.5 s onward for further analysis. In addition, the standard deviations of the measurements obtained using all the station pairs will be used to quantify relative uncertainty during depth inversion (see “Inversion for shear wave velocity” section).

Ambient noise group velocity maps and resolution tests

The data from dispersion curves were transformed to group velocities in a 2D grid using damped least squares inversion. The group velocities at each period were forward modeled based on ray theory and superposition principle, where the arrival time of a surface wave group in a cross-correlation waveform is the sum of its travel times in 2D bins along the station pair path. For each path, the traveltime perturbation with respect to the reference traveltime in mean model depends linearly on the path length and the relative group velocity perturbation with respect to mean group velocity in each 2D bin. Then the L2 norm of the misfit between observed and modeled group velocities are inverted linearly using damped least squares. The forward modeling and the inversion methodology is well described and implemented in a plugin of MSNoise (MSNoise-Tomo)^44,65. The 2D grid for inversion was defined by dividing the Koyna-Warna region into 2.5 km 2.5 km bins, and the weights for damping and smoothing were calculated for each period using l-curve⁶⁶ (Fig. S3). The data covariance matrix was diagonally populated with FTAN amplitudes of the dispersion picks calculated (see “Computation of dispersion curves” section), in order to assign weights to the measurements based on the quality of the picks. The parameters for the generation of group velocity maps using MSNoise-Tomo are provided in Tables Fig. S1, Fig. S2, and Fig. S3.

[See PDF for image]

Fig. 7

Spike resolution tests for Rayleigh wave group velocity maps (a,c) and Love wave group velocity maps (b,d) in 0.5 s and 5.0 s periods. Rayleigh wave group velocity maps (e,g) and Love wave group velocity maps (f,h) calculated from ZZ and TT cross-correlations respectively in 0.5 s and 5.0 s periods. The bins with less than three station-pair are masked. The map has been overlayed on the topography of the Koyna–Warna region⁶⁷. A common color scale has been used for both Rayleigh and Love wave group velocities at each period. The outline of the Koyna and Warna rivers as well as their respective reservoirs are also overlayed on the map (black). The stations used are shown in red triangles in each sub-figure. The maps in this figure were generated using Python-v3.8.18 scripts with MSNoise-v1.3, MSNoise-Tomo⁴⁴, GMT-6⁶⁸, and pygmt-v0.9¹¹, and combined using Inkscape-v1.2.2¹², please see Data Availability section for more details.

In order to test the resolution of the group velocity maps, we performed spike resolution tests at each period using spike anomalies of size 5 km 5 km each, separated by 5 km. The spike resolution tests or sparse checkerboard allow clear identification of the distortion in anomaly shapes compared to conventional checkerboards⁶⁹. Smearing of the anomalies is observed near the boundaries of the region at all periods for both Rayleigh and Love wave maps (Fig. 7a–d). The well-recovered region gradually shrinks in size with increasing period, such that at 5.0 s period only the central part of the Koyna-Warna region comprising of the Warna reservoir and southern Koyna reservoir is well-recovered. The resolution is roughly similar for both Rayleigh and Love wave maps. The amplitude of the anomalies also decreases with increasing period. This pattern is likely due to fewer available station-pair paths at longer periods (Fig. S4).

Based on the resolution test, the poorly recovered regions at all periods have been masked during the generation of group velocity maps (Fig. 7e–h). In addition, regions that show abnormally high group velocities for shorter periods (> 3.75 km/s for Love wave and > 3.50 km/s for Rayleigh wave) that are unrealistic for crustal scale from regional studies^70,71 have been masked. In general, the Love wave velocities are higher than Rayleigh wave velocities across all periods. In most parts of the region, the group velocities increase with increasing period because longer period waves sample greater depths and the subsurface seismic velocities tend to generally increase with depth.

Inversion for shear wave velocity

The dispersion curve obtained from each bin of the group velocity maps was used to invert for 1D shear wave velocity in depth. In total, 667 bins of 2.5 km 2.5 km size from Rayleigh wave maps and 664 bins from Love wave maps maps were used. The depth inversion was performed using the nearest neighbor algorithm⁷², which is implemented in the dinver program of the Geopsy software package for ambient noise processing⁴⁶. The dispersion curves obtained from Rayleigh wave maps and Love wave maps were separately inverted for shear wave velocity in vertical direction () and shear wave velocity in horizontal direction () respectively. The relative standard deviations of the dispersion curves (Fig. 6c,d) across different periods were used as standard errors during the inversion.

As with other inversion problems in earth sciences, the inversion of surface wave dispersion curves for shear wave velocity is ill-posed and suffers from non-uniqueness of the solution. In order to overcome these challenges to a certain extent, a workflow wherein resampling of the dispersion curves based on wavelength and experimenting with multiple parametrization of the 1D model has been suggested⁷³. Following this workflow, each dispersion curve was resampled in terms of wavelength () and six different parameterizations (three based on fixed number of layers and three based on fixed thickness ratio of successive layers) for 1D model were defined using swprepost package⁴⁵. The resampled dispersion curve was inverted in geopsy-dinver by applying the nearest neighbor algorithm to a starting model generated following each parametrization. At the end of 6 inversions for each dispersion curve, 20 lowest misfit models from the 6 types of parameterizations were weighted-averaged with respect to their misfits to obtain the best-fit 1D model. This process was followed for both 1D and inversions. The full list of parameters related to the 1D model inversions are provided in Table S4.

After best-fit and 1D models were obtained for each lateral bin, the isotropic model () and the radial anisotropy () were calculated as follows:

1

2

where and are the misfit errors of and model estimations respectively.

The models and the for all the bins were smoothed in the lateral direction using a 2D Gaussian smoothing with kernel width of 2 bins (5 km), which is equivalent to the spatial resolution found laterally (Fig. 7a-d). In the vertical direction, a 1D Gaussian smoothing with a width of 300 m was applied. This width is less than the minimum resolvable layer thickness of the 1D model ()⁷⁴, which also corresponds to the vertical resolution. The fundamental mode surface waves are not sensitive to depths larger than ⁷⁵, therefore, this value was used as the maximum resolvable depth. Finally, a topographical correction based on Scripps’ SRTM15+V2.6 grid⁷⁶ was applied to the and results using pygmt. The model volume was sliced laterally and vertically to obtain the horizontal sections and the vertical cross-sections respectively. In the previous sections, has been referred as just .

Acknowledgements

The digital seismic waveform data, collected using the 3-Component Taurus Standalone recording system, belongs to the CSIR-National Geophysical Research Institute (NGRI), and permission for its use has been granted vide reference no. Lib/2025/Pub-42. We acknowledge the contributions of the Controlled Source Seismic (CSS) field crew in data acquisition. The financial assistance provided by the Ministry of Earth Sciences (MoES), New Delhi, is gratefully acknowledged. The authors extend their special thanks to Dr. H. K. Gupta for his encouragement and kind support in facilitating this study as a pilot project and Dr. M.M. Dixit for successfully executing the pilot project. We also thank Dr. Rahul Biswas for his valuable comments, which helped improve the manuscript. Finally, we are thankful to the two anonymous reviewers for their constructive comments that helped improve the manuscript.

Author contributions

P.R. conducted the analysis and wrote the initial manuscript draft as well as revised the manuscript. P.K. conceived idea for the analysis and provided computational resources. B.M. and P.K prepared the data, reviewed the results, and prepared some of the figures. All authors revised the initial manuscript draft.

Data availability

The Supplementary Material for this article includes four additional figures, four tables and a description of the contents. It also contains description of all the datasets and scripts necessary to generate the figures, which are uploaded to figshare.org at https://figshare.com/s/e739c54537ba304ee0f3.

Declarations

Competing interests

The authors declare no competing interests.

Supplementary Information

The online version contains supplementary material available at https://doi.org/10.1038/s41598-025-12099-z.

Publisher’s note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

References

1. Gupta, HK. Artificial water reservoir-triggered seismicity (rts): Most prominent anthropogenic seismicity. Surv. Geophys.; 2022; 43, pp. 619-659.2022SGeo..43.619G [DOI: https://dx.doi.org/10.1007/s10712-021-09675-z]

2. Narain, H; Gupta, H. Koyna earthquake. Nature; 1968; 217, pp. 1138-1139.1968Natur.217.1138N [DOI: https://dx.doi.org/10.1038/2171138a0]

3. Guha, S; Gosavi, P; Padale, J; Marwadi, S. An earthquake cluster at Koyna. Bull. Seismol. Soc. Am.; 1971; 61, pp. 297-315. [DOI: https://dx.doi.org/10.1785/BSSA0610020297]

4. Langston, CA. A body wave inversion of the Koyna, India, earthquake of December 10, 1967, and some implications for body wave focal mechanisms. J. Geophys. Res.; 1976; 1896–1977, 81 pp. 2517-2529.1976JGR..81.2517L [DOI: https://dx.doi.org/10.1029/JB081i014p02517]

5. Langston, CA. Source inversion of seismic waveforms: The Koyna, India, earthquakes of 13 September 1967. Bull. Seismol. Soc. Am.; 1981; 71, pp. 1-24. [DOI: https://dx.doi.org/10.1785/BSSA0710010001]

6. Langston, CA; Franco-Spera, M. Modeling of the Koyna, India, aftershock of 12 December 1967. Bull. Seismol. Soc. Am.; 1985; 75, pp. 651-660. [DOI: https://dx.doi.org/10.1785/BSSA0750030651]

7. Gupta, HK. Review: Reservoir triggered seismicity (rts) at Koyna, India, over the past 50 yrs. Bull. Seismol. Soc. Am.; 2018; 108, pp. 2907-2918. [DOI: https://dx.doi.org/10.1785/0120180019]

8. Rao, NP; Shashidhar, D. Periodic variation of stress field in the Koyna–Warna reservoir triggered seismic zone inferred from focal mechanism studies. Tectonophysics; 2016; 679, pp. 29-40.2016Tectp.679..29R [DOI: https://dx.doi.org/10.1016/j.tecto.2016.04.036]

9. Shashidhar, D et al. A catalogue of earthquakes in the Koyna–Warna region, western India (2005–2017). J. Geol. Soc. India; 2019; 93, pp. 7-24. [DOI: https://dx.doi.org/10.1007/s12594-019-1115-y]

10. Guha, S. et al. Catalogue of earthquakes (>= m 3.0) in peninsular India. Atomic Energy Regulatory Board—Technical Reports. https://inis.iaea.org/records/04xwa-pdq75 (1993).

11. Uieda, L. et al. PyGMT: A Python Interface for the Generic Mapping Tools. https://doi.org/10.5281/zenodo.7772533 (2023).

12. Inkscape Project. Inkscape Vector Graphics Editor, Version 1.2.2. https://inkscape.org (2022).

13. Gupta, HK; Rastogi, BK; Narain, H. Common features of the reservoir-associated seismic activities. Bull. Seismol. Soc. Am.; 1972; 62, pp. 481-492. [DOI: https://dx.doi.org/10.1785/BSSA0620020481]

14. Pandey, AP; Chadha, R. Surface loading and triggered earthquakes in the Koyna–Warna region, western India. Phys. Earth Planet. Inter.; 2003; 139, pp. 207-223.2003PEPI.139.207P [DOI: https://dx.doi.org/10.1016/j.pepi.2003.08.003]

15. Catchings, RD; Dixit, MM; Goldman, MR; Kumar, S. Structure of the Koyna–Warna seismic zone, Maharashtra, India: A possible model for large induced earthquakes elsewhere. J. Geophys. Res. Solid Earth; 2015; 120, pp. 3479-3506.2015JGRB.120.3479C [DOI: https://dx.doi.org/10.1002/2014JB011695]

16. Gahalaut, VK; Kalpna, A; Singh, SK. Fault interaction and earthquake triggering in the Koyna–Warna region, India. Geophys. Res. Lett.; 2004; 31, 818. [DOI: https://dx.doi.org/10.1029/2004GL019818]

17. Talwani, P. Seismotectonics of the Koyna–Warna area, India. Pure Appl. Geophys.; 1997; 150, pp. 511-550.1997PApGe.150.511T [DOI: https://dx.doi.org/10.1007/s000240050091]

18. Durá-Gómez, I; Talwani, P. Hydromechanics of the Koyna–Warna region, India. Pure Appl. Geophys.; 2010; 167, pp. 183-213.2010PApGe.167.183D [DOI: https://dx.doi.org/10.1007/s00024-009-0012-5]

19. Courtillot, V et al. Deccan flood basalts at the cretaceous/tertiary boundary?. Earth Planet. Sci. Lett.; 1986; 80, pp. 361-374.1986E&PSL.80.361C1:CAS:528:DyaL2sXitFekug%3D%3D [DOI: https://dx.doi.org/10.1016/0012-821X(86)90118-4]

20. Patro, BPK; Sarma, SVS. Trap thickness and the subtrappean structures related to mode of eruption in the deccan plateau of India: Results from magnetotellurics. Earth Planets Space; 2007; 59, pp. 75-81.2007EP&S..59..75P [DOI: https://dx.doi.org/10.1186/BF03352679]

21. Gupta, H et al. Investigations related to scientific deep drilling to study reservoir-triggered earthquakes at Koyna, India. Int. J. Earth Sci.; 2015; 104, pp. 1511-1522.1:CAS:528:DC%2BC2cXitFKltrfN [DOI: https://dx.doi.org/10.1007/s00531-014-1128-0]

22. Misra, S et al. Granite-gneiss basement below deccan traps in the Koyna region, western India: Outcome from scientific drilling. J. Geol. Soc. India; 2017; 90, pp. 776-782.2017cody.book...M1:CAS:528:DC%2BC1cXjs1ejsw%3D%3D [DOI: https://dx.doi.org/10.1007/s12594-017-0790-9]

23. Halder, P; Shukla, MK; Kumar, K; Sharma, A. Understanding the mechanism of shallow crustal fluid-rock interaction in the deccan trap basement rocks and its significance in the koyna-warna seismogenic region, western india. Phys. Chem. Earth A/B/C; 2025; [DOI: https://dx.doi.org/10.1016/j.pce.2025.103870]

24. Goswami, D; Hazarika, P; Roy, S. In situ stress orientation from 3 km borehole image logs in the koyna seismogenic zone, western India: Implications for transitional faulting environment. Tectonics; 2020; 39, e2019005647.2020Tecto.3905647G [DOI: https://dx.doi.org/10.1029/2019TC005647]

25. Modak, K; Rohilla, S; Podugu, N; Goswami, D; Roy, S. Fault associated with the 1967 m 6.3 koyna earthquake, India: A review of recent studies and perspectives for further probing. J. Asian Earth Sci. X; 2022; 8, 100123. [DOI: https://dx.doi.org/10.1016/j.jaesx.2022.100123]

26. Agrawal, P; Pandey, O; Chetty, T. Aeromagnetic anomalies, lineaments and seismicity in Koyna–Warna region. J. Indian Geophys. Union; 2004; 8, pp. 229-242.

27. Borah, UK et al. Role of fluid on seismicity of an intra-plate earthquake zone in western India: An electrical fingerprint from magnetotelluric study. Earth Planets Space; 2023; 75, 149.2023EP&S..75.149B [DOI: https://dx.doi.org/10.1186/s40623-023-01905-5]

28. Arora, K et al. Lineaments in deccan basalts: The basement connection in the Koyna–Warna rts region. Bull. Seismol. Soc. Am.; 2018; 108, pp. 2919-2932. [DOI: https://dx.doi.org/10.1785/0120180011]

29. Gupta, HK et al. Location of the pilot borehole for investigations of reservoir triggered seismicity at Koyna, India. Gondwana Res.; 2017; 42, pp. 133-139.2017GondR.42.133G [DOI: https://dx.doi.org/10.1016/j.gr.2016.10.014]

30. Krishna, V. G., Kaila, K. L. & Reddy, P. R. Synthetic seismogram modeling of crustal seismic record sections from the Koyna dss profiles in the western India. In Properties and Processes of Earth’s Lower Crust 143–157. https://doi.org/10.1029/GM051p0143 (American Geophysical Union (AGU), 1989).

31. Tiwari, V; Vyaghreswara Rao, M; Mishra, D. Density inhomogeneities beneath deccan volcanic province, India as derived from gravity data. J. Geodyn.; 2001; 31, pp. 1-17. [DOI: https://dx.doi.org/10.1016/S0264-3707(00)00015-6]

32. Vasanthi, A; SatishKumar, K. Understanding conspicuous gravity low over the Koyna–Warna seismogenic region (Maharashtra, India) and earthquake nucleation: A paradigm shift. Pure Appl. Geophys.; 2016; 173, pp. 1933-1948.2016PApGe.173.1933V [DOI: https://dx.doi.org/10.1007/s00024-015-1237-0]

33. Patro, PK; Borah, U; Babu, GA; Veeraiah, B; Sarma, SVS. Ground electrical and electromagnetic studies in Koyna–Warna region, India. J. Geol. Soc. India; 2017; 90, pp. 711-719. [DOI: https://dx.doi.org/10.1007/s12594-017-0780-y]

34. Dixit, MM et al. Seismicity, faulting, and structure of the Koyna–Warna seismic region, western India from local earthquake tomography and hypocenter locations. J. Geophys. Res. Solid Earth; 2014; 119, pp. 6372-6398.2014JGRB.119.6372D [DOI: https://dx.doi.org/10.1002/2014JB010950]

35. Kumar, S; Kumar, P. One-dimensional velocity model, station correction and earthquake relocation of local earthquakes in the koyna-warna region, india. Pure Appl. Geophys.; 2019; 176, pp. 4761-4782.2019PApGe.176.4761K [DOI: https://dx.doi.org/10.1007/s00024-019-02264-7]

36. Shapiro, NM; Campillo, M. Emergence of broadband Rayleigh waves from correlations of the ambient seismic noise. Geophys. Res. Lett.; 2004; 31, 491. [DOI: https://dx.doi.org/10.1029/2004GL019491]

37. Ranjan, P; Stehly, L. Estimation of seismic attenuation from ambient noise coda waves: Application to the hellenic subduction zone. Bull. Seismol. Soc. Am.; 2024; 114, pp. 2065-2082. [DOI: https://dx.doi.org/10.1785/0120230265]

38. Lehujeur, M et al. Reservoir imaging using ambient noise correlation from a dense seismic network. J. Geophys. Res. Solid Earth; 2018; 123, pp. 6671-6686.2018JGRB.123.6671L [DOI: https://dx.doi.org/10.1029/2018JB015440]

39. Wang, Y; Lin, F-C; Schmandt, B; Farrell, J. Ambient noise tomography across mount st. helens using a dense seismic array. J. Geophys. Res. Solid Earth; 2017; 122, pp. 4492-4508.2017JGRB.122.4492W [DOI: https://dx.doi.org/10.1002/2016JB013769]

40. Esteve, C et al. The seismic signature and geothermal potential of the Schwechat depression in the Vienna basin, Austria, from ambient noise tomography. Geothermics; 2025; 127, 103211. [DOI: https://dx.doi.org/10.1016/j.geothermics.2024.103211]

41. Jiang, C; Denolle, MA. Pronounced seismic anisotropy in kanto sedimentary basin: A case study of using dense arrays, ambient noise seismology, and multi-modal surface-wave imaging. J. Geophys. Res. Solid Earth; 2022; 127, e2022024613.2022JGRB.12724613J [DOI: https://dx.doi.org/10.1029/2022JB024613]

42. Wu, S-M et al. Crustal characterization of the hengill geothermal fields: Insights from isotropic and anisotropic seismic noise imaging using a 500-node array. J. Geophys. Res. Solid Earth; 2024; 129, e2024028915. [DOI: https://dx.doi.org/10.1029/2024JB028915]

43. Rohilla, S; Rao, NP; Gerstoft, P; Shashidhar, D; Satyanarayana, HVS. Shear-wave velocity structure of the Koyna–Warna region in western India using ambient noise correlation and surface-wave dispersion. Bull. Seismol. Soc. Am.; 2015; 105, pp. 473-479. [DOI: https://dx.doi.org/10.1785/0120140091]

44. Lecocq, T; Caudron, C; Brenguier, F. Msnoise, a python package for monitoring seismic velocity changes using ambient seismic noise. Seismol. Res. Lett.; 2014; 85, pp. 715-726. [DOI: https://dx.doi.org/10.1785/0220130073]

45. Vantassel, J. jpvantassel/swprepost: v2.0.0. https://doi.org/10.5281/zenodo.6558333 (2022).

46. Wathelet, M et al. Geopsy: A user-friendly open-source tool set for ambient vibration processing. Seismol. Res. Lett.; 2020; 91, pp. 1878-1889. [DOI: https://dx.doi.org/10.1785/0220190360]

47. Duraiswami, RA et al. Volcanology and lava flow morpho-types from the Koyna–Warna region, western deccan traps, India. J. Geol. Soc. India; 2017; 90, pp. 742-747.1:CAS:528:DC%2BC1cXjs1ejug%3D%3D [DOI: https://dx.doi.org/10.1007/s12594-017-0785-6]

48. Hua, Y; Zhao, D; Toyokuni, G; Xu, Y. Tomography of the source zone of the great 2011 Tohoku earthquake. Nat. Commun.; 2020; 11, 1163.2020NatCo.11.1163H [DOI: https://dx.doi.org/10.1038/s41467-020-14745-8] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/32127532][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC7054414]

49. Ranjan, P; Konstantinou, K. Local earthquake tomography of the aegean crust: Implications for active deformation, large earthquakes, and arc volcanism. Tectonophysics; 2024; 880, 230331. [DOI: https://dx.doi.org/10.1016/j.tecto.2024.230331]

50. Rubin, AM; Ampuero, J-P. Earthquake nucleation on (aging) rate and state faults. J. Geophys. Res. Solid Earth; 2005; 110, 686. [DOI: https://dx.doi.org/10.1029/2005JB003686]

51. Uenishi, K; Rice, JR. Universal nucleation length for slip-weakening rupture instability under nonuniform fault loading. J. Geophys. Res. Solid Earth; 2003; 108, 681. [DOI: https://dx.doi.org/10.1029/2001JB001681]

52. Thakur, P; Huang, Y. The effects of precursory velocity changes on earthquake nucleation and stress evolution in dynamic earthquake cycle simulations. Earth Planet. Sci. Lett.; 2024; 637, 118727.1:CAS:528:DC%2BB2cXpsFyqs7g%3D [DOI: https://dx.doi.org/10.1016/j.epsl.2024.118727]

53. Pandey, A; Gahalaut, VK; Yadav, RK; Catherine, JK; Naidu, MS. Conundrum of sustained seismicity and low deformation in the Koyna–Warna region, India. Bull. Seismol. Soc. Am.; 2025; [DOI: https://dx.doi.org/10.1785/0120240196]

54. Ebrahimi, M; Tatar, M; Aoudia, A; Guidarelli, M. Loading effects beneath the Gotvand-e olya reservoir (south-west of Iran) deduced from ambient noise tomography. Geophys. J. Int.; 2017; 212, pp. 229-243.2018GeoJI.212.229E [DOI: https://dx.doi.org/10.1093/gji/ggx411]

55. Dong, S et al. Seismic evidence for fluid-driven pore pressure increase and its links with induced seismicity in the Xinfengjiang reservoir, South China. J. Geophys. Res. Solid Earth; 2022; 127, e2021023548.2022JGRB.12723548D [DOI: https://dx.doi.org/10.1029/2021JB023548]

56. Podugu, N; Goswami, D; Akkiraju, VV; Roy, S. In-situ physical and elastic properties of archaean basement granitoids in the Koyna seismogenic zone, western India from 3 km downhole geophysical well logs: Implications for water percolation at depth. Tectonophysics; 2023; 848, 229725. [DOI: https://dx.doi.org/10.1016/j.tecto.2023.229725]

57. Modak, K. Evidence of water percolation in granitoid basement in koyna seismogenic zone: Implications for reservoir triggered seismicity. J. Earth Syst. Sci.; 2024; 133, 163. [DOI: https://dx.doi.org/10.1007/s12040-024-02380-6]

58. Kato, A; Ben-Zion, Y. The generation of large earthquakes. Nat. Rev. Earth Environ.; 2021; 2, pp. 26-39.2021NRvEE..2..26K [DOI: https://dx.doi.org/10.1038/s43017-020-00108-w]

59. Peterson, J. R. Observations and Modeling of Seismic Background Noise. https://doi.org/10.3133/ofr93322 (1993).

60. Krischer, L et al. ObsPy: A bridge for seismology into the scientific python ecosystem. Comput. Sci. Discov.; 2015; 8, 014003.

61. Hunter, JD. Matplotlib: A 2d graphics environment. Comput. Sci. Eng.; 2007; 9, pp. 90-95. [DOI: https://dx.doi.org/10.1109/MCSE.2007.55]

62. Lin, F-C; Moschetti, MP; Ritzwoller, MH. Surface wave tomography of the western united states from ambient seismic noise: Rayleigh and love wave phase velocity maps. Geophys. J. Int.; 2008; 173, pp. 281-298.2008GeoJI.173.281L [DOI: https://dx.doi.org/10.1111/j.1365-246X.2008.03720.x]

63. Levshin, A. et al.Recording, Identification, and Measurement of Surface Wave Parameters 131–182 (Springer, 1989).

64. Schwardt, M; Pilger, C; Gaebler, P; Hupe, P; Ceranna, L. Natural and anthropogenic sources of seismic, hydroacoustic, and infrasonic waves: Waveforms and spectral characteristics (and their applicability for sensor calibration). Surv. Geophys.; 2022; 43, pp. 1265-1361.2022SGeo..43.1265S [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/35911621][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC9309596]

65. Mordret, A et al. Near-surface study at the valhall oil field from ambient noise surface wave tomography. Geophys. J. Int.; 2013; 193, pp. 1627-1643.2013GeoJI.193.1627M [DOI: https://dx.doi.org/10.1093/gji/ggt061]

66. Calvetti, D; Morigi, S; Reichel, L; Sgallari, F. Tikhonov regularization and the l-curve for large discrete ill-posed problems. J. Comput. Appl. Math.; 2000; 123, pp. 423-446.2000JCoAm.123.423C1798536 [DOI: https://dx.doi.org/10.1016/S0377-0427(00)00414-3] 0977.65030

67. NOAA National Centers for Environmental Information. ETOPO 2022 15 Arc-Second Global Relief Model (2022).

68. Wessel, P et al. The generic mapping tools version 6. Geochem. Geophys. Geosyst.; 2019; 20, pp. 5556-5564.2019GGG..20.5556W [DOI: https://dx.doi.org/10.1029/2019GC008515]

69. Rawlinson, N; Spakman, W. On the use of sensitivity tests in seismic tomography. Geophys. J. Int.; 2016; 205, pp. 1221-1243.2016GeoJI.205.1221R [DOI: https://dx.doi.org/10.1093/gji/ggw084]

70. Dey, S; Ghosh, M; Mitra, S. 3d shear-wave velocity structure of the crust and upper mantle beneath India, Himalaya and Tibet. J. Geophys. Res. Solid Earth; 2024; 129, e2023028292. [DOI: https://dx.doi.org/10.1029/2023JB028292]

71. Ghosh, M; Dey, S; Mitra, S. Love-wave group-velocity tomography of India–Tibet: Insights into radially anisotropic s-wave velocity structure. Geophys. J. Int.; 2025; 241, pp. 901-918. [DOI: https://dx.doi.org/10.1093/gji/ggaf074]

72. Sambridge, M. Geophysical inversion with a neighbourhood algorithm—i. Searching a parameter space. Geophys. J. Int.; 1999; 138, pp. 479-494.1999GeoJI.138.479S [DOI: https://dx.doi.org/10.1046/j.1365-246X.1999.00876.x]

73. Vantassel, JP; Cox, BR. Swinvert: A workflow for performing rigorous 1-d surface wave inversions. Geophys. J. Int.; 2020; 224, pp. 1141-1156.2021GeoJI.224.1141V [DOI: https://dx.doi.org/10.1093/gji/ggaa426]

74. Cox, BR; Teague, DP. Layering ratios: A systematic approach to the inversion of surface wave data in the absence of a priori information. Geophys. J. Int.; 2016; 207, pp. 422-438.2016GeoJI.207.422C [DOI: https://dx.doi.org/10.1093/gji/ggw282]

75. Garofalo, F et al. Interpacific project: Comparison of invasive and non-invasive methods for seismic site characterization. Part i: Intra-comparison of surface wave methods. Soil Dyn. Earthq. Eng.; 2016; 82, pp. 222-240. [DOI: https://dx.doi.org/10.1016/j.soildyn.2015.12.010]

76. Tozer, B et al. Global bathymetry and topography at 15 arc sec: Srtm15+. Earth Space Sci.; 2019; 6, pp. 1847-1864.2019E&SS..6.1847T [DOI: https://dx.doi.org/10.1029/2019EA000658]