Loading effects beneath the Gotvand-e Olya Reservoir (south-west of Iran) deduced from ambient noise tomography

Loading effects beneath the Gotvand-e Olya Reservoir (south-west of Iran) deduced from ambient... Abstract In order to define the precise shallow velocity structure beneath the second largest dam reservoir in Iran and to understand the loading effects on the underlying crust, the shear wave velocity of the shallow structure beneath the Gotvand-e Olya (hereinafter referred to as Gotvand) reservoir is determined through the inversion of group velocities obtained from seismic ambient noise tomography, using continuous data from 10 stations of a local network, installed to monitor the induced seismicity in the region surrounding the Gotvand and Masjed Soleyman dams for potential hazard. We obtained Rayleigh waves from cross-correlation of waveforms recorded 10 months before and the same duration after impoundment of the Gotvand reservoir and calculated the group velocity from dispersion analysis in the period range 2–8 s. The group velocity dispersion curves are used to produce 2-D group velocity tomographic maps. The resulting tomographic maps at short periods are well correlated with subsurface geological features and delineate distinct low- and high-velocity zones separated mainly by geological boundaries. The 3-D shear wave velocity structure provides detailed information about the crustal features underneath the reservoir. The results are consistent with the lithology of the region, and attest that ambient noise tomography (ANT) can be used for detailed studies of the velocity structure and lithology at shallow depths using continuous data from a dense local seismic network. An increase of shear wave velocity is observed at the deep parts (4–6 km) underneath the reservoir after impoundment of the dam, which could be caused by the changes in rocks properties after impoundment. However, at shallow depths (2–4 km), a decrease of Vs velocity is observed that can be associated to the penetration of water after the impoundment. Changes of material strength underneath the Gotvand reservoir, calculated using loading effects of the reservoir by considering the pore pressure are in good agreement with the observed shear wave velocity changes across the reservoir. Tomography14, Crustal imaging, Induced seismicity, Seismic noise 1 INTRODUCTION Building a dam causes severe stress load on the crust underneath the reservoir when the valley behind it is filled. The pore pressure of penetrated fluids can also change the physical properties of the Earth's crust via lowering the strength of rocks (Talwani & Acree 1985; Talwani 1997). Gotvand Dam with a total height of 180 m is the highest rock-fill dam with central clay core in Iran. The Gotvand reservoir with a total volume of 4500 Milion Cubic Meters and an area of ∼96.5 km2 at the level of 243 m from sea level is the second large reservoir in Iran after the Karkhe dam. Gotvand dam is located in the Zagros Mountains of western Iran, which is known as one of the most seismically active regions along the Alp-Himalayan belt. The Gotvand dam is located in the Zagros simply fold and thrust belt (SFB). Shortening at the surface in this area is accommodated by parallel trains of anticlines and synclines, while deformation in depth is modulated by parallel NW-SE high angle basement reverse faults and roughly N–S, right-lateral surface faults like Kazerun, Borazjan, Karebas, Sabzpushan and Sarvestan (Fig. 1). The N-dipping basement reverse faults have caused significant changes in the stratigraphic level across certain folds in the study region (Berberian 1995). Lahbari known also as Mountain Front Fault (MFF), which marks the southwestern limit of the surface exposure of the Asmari limestone is one of the most important N-dipping faults in the study region. Figure 1. View largeDownload slide Map of the Zagros Fold and Thrust belt with SRTM digital topography coloured to emphasize the high mountains (green) and the deep valleys (yellow) that penetrate the range. Seismicity (IIEES catalogue ML > 4.5) is shown by red circles. The white rectangle is the location of the study region. Faults from Hessami Azar et al. (2003). MFF, Mountain Front Fault; ZFF, Zagros Foredeep Faults; MRF, Main Recent Fault; BZ, Borazjan; KZ, Kazerun; KB, Karebas; SB, Sabzpushan. Figure 1. View largeDownload slide Map of the Zagros Fold and Thrust belt with SRTM digital topography coloured to emphasize the high mountains (green) and the deep valleys (yellow) that penetrate the range. Seismicity (IIEES catalogue ML > 4.5) is shown by red circles. The white rectangle is the location of the study region. Faults from Hessami Azar et al. (2003). MFF, Mountain Front Fault; ZFF, Zagros Foredeep Faults; MRF, Main Recent Fault; BZ, Borazjan; KZ, Kazerun; KB, Karebas; SB, Sabzpushan. Topographic analysis and seismotectonic studies show that activation of blind thrust faults in the basement have resulted in the evolution of asymmetric folds in the sedimentary cover (Sattarzadeh 1997; McQuarrie 2004; Mouthereau et al. 2006). Microseismicity of the region is dominated by reverse faulting, mainly confined at depths between 5 and 20 km. However the distribution of teleseismic centroid depths shows that most of the seismic strain is released at depths of 5–10 km (Nissen et al. 2011) mainly within the lower sedimentary cover. O’Brien (1950) and Dunnington (1968) compiled the Zagros stratigraphic column by using regional stratigraphy. They divided it into five distinct units, which from bottom to top, are the Precambrian metamorphic Basement Group; the Lower Mobile Group, which is considered as a decollement level and is formed by Late Proterozoic–Early Cambrian Hormuz salt, the Competent Group, which consists of Cambrian to Lower Miocene Formation; the Upper Mobile Group (Gachsaran Miocene salt); the Passive Group (mostly clastic sediments) at the top. In the study region (Fig. 2), Gachsaran evaporate is the main formation in the Upper Mobile Group. James & Wynd (1965) and Setudehnia (1972) described this formation as composed by seven distinct sections, consisting of anhydrite with bituminous shale at the bottom, salt and limestone at the middle and alternating marl at the top. Gachsaran Formation is incompetent and it is highly subject to extreme mobility as a result of differential pressure. Agha Jari and Bakhtyari formations are the main clastic deposits in the Passive Group at the top (Fig. 2). The principal lithology of Agha Jari formation consists of gypsum veins, calcareous sandstone, red marl and siltstone and the Bakhtiari Formation is terrigenous, clastic unit ranging in grain size from a silt grade to boulder conglomerate (Setudehnia 1972). Figure 2. View largeDownload slide Geologic map and cross-section view of the Gotvand and Masjed Soleyman Reservoirs (National Geoscience Database of Iran). Figure 2. View largeDownload slide Geologic map and cross-section view of the Gotvand and Masjed Soleyman Reservoirs (National Geoscience Database of Iran). In the last decade, the methodology called ambient noise tomography (ANT) has been developed to produce detailed seismic velocity models for the crust and upper mantle (Shapiro et al. 2005; Bensen et al. 2007). Cross-correlation of ambient seismic noise can produce short period surface waves between station pairs and thus can provide more accurate and detailed constraints on upper crustal structures than traditional surface wave tomography. It provides much higher resolution without need of using any earthquake surface-wave data (Yao et al. 2006; Guo et al. 2009). ANT is based on the estimation of empirical Green’s functions (EGF) between two receivers by cross-correlating the continuous ambient seismic noise recorded at the two receivers (Lobkis &Weaver 2001; Snieder 2004). Recently ambient noise tomography has been widely used from regional scale (Shapiro et al. 2005; Yao et al. 2006) to local scale such as volcanic regions (Brenguier et al. 2007; Koulakov et al. 2014; Droznin et al. 2015). This paper presents the results of ambient noise analysis of data collected about 10 months before and after the impoundment of the Gotvand dam by 10 seismological stations belonging to the dam seismic monitoring system. Gotvand dam has been continuously monitored by a local seismological network (Fig. 3) since 2006, about 5 yr before impoundment of the reservoir. However, the continuous seismic data was accessible only from 2011 March 20, about 10 months before starting the impoundment. The same duration of 10 months data set recorded after the impoundment and this until 2012 November 20 was used. The Gotvand reservoir impoundment actually was started on 2012 January 20 when a considerable increase in water level is observed. The impoundment of the Masjed Soleyman dam is started on 2001 December 12 and finished on 2002 June 25. There was no independent permanent seismological network monitoring the Masjed Soleyman reservoir prior to its impoundment. Figure 3. View largeDownload slide Topographic map of the study area including Gotvand and Masjed Soleyman dams reservoirs. Seismological stations of two Gotvand and Masjed Soleyman seismic networks are presented in black triangles. Cities are shown in black squares. Distribution of the reservoir lake at 243 m above the sea level is presented in blue. The black solid lines represent the location of vertical cross-sections. Yellow star represent the location of the gridpoint used to calculate sensitivity kernel (Fig. 7) for the fundamental Rayleigh wave group velocity. Figure 3. View largeDownload slide Topographic map of the study area including Gotvand and Masjed Soleyman dams reservoirs. Seismological stations of two Gotvand and Masjed Soleyman seismic networks are presented in black triangles. Cities are shown in black squares. Distribution of the reservoir lake at 243 m above the sea level is presented in blue. The black solid lines represent the location of vertical cross-sections. Yellow star represent the location of the gridpoint used to calculate sensitivity kernel (Fig. 7) for the fundamental Rayleigh wave group velocity. The purpose of this study is to determine the shallow structure beneath the Gotvand reservoir before and after impoundment of the dam and to characterize the possible changes in the shallow crust due to loading effect. We measure the group velocity of the cross-correlation Rayleigh waves at periods from 2 to 8 s. Then a surface wave tomography is performed to obtain high-resolution group velocity maps for the study area beneath the Gotvand reservoir. Group velocities from these maps will be inverted to produce high resolution 3-D model of shear wave velocity, which can provide important constraints for better understanding the loading effects due to reservoir impoundment on shallow structure characterization beneath the reservoir. 2 DATA PROCESSING We used 20 months of continuous vertical-component seismic data recorded by the Gotvand-Masjed Soleyman seismological network (Fig. 3) between 2011 March 20 and 2012 November 20. This seismic network consists of 10 medium band Trilium 40 seismometers equipped with 24-bit Taurus recorders both manufactured by Nanometrics company. All stations operate in continuous mode with a sampling rate of 80 sps. This seismological network belongs to Gotvand project of Iran Water and Power Resource Development Co. (IWPCO) and to Masjed Solyman dam of Khuzestan Water and Power Authority (KWPA). It has been designed and installed by International Institute of Earthquake Engineering and Seismology (IIEES) and is operated by IIEES since 2006. About 20 months of data, collected during 10 months before and after impoundment of the Gotvand dam are used to image shallow velocity structure of the area beneath the Gotvand reservoir from ambient noise tomography. The same procedure as described in Bensen et al. (2007) is applied for ambient noise data processing. Using only vertical components of ambient noise, we could only extract Rayleigh waves from waveform cross-correlations. First, instrument responses, mean and trend are removed from the continuous data that are then bandpass filtered in the period range 1–30 s. The contamination of ambient noise recordings by earthquakes and non-stationary noise sources near the stations is partly removed by applying the procedure of ‘temporal normalization’ which is called running-absolute-mean normalization (Bensen et al.2007). This normalization method computes the running average of absolute value of the seismic waveform in a normalization time window of fixed length and weights at the centre of the window by the inverse of this average (Bensen et al. 2007). This is followed by spectral whitening which is applied to flatten spectra over the entire period band (1–30 s). Then daily cross-correlations of continuous noise data in the period range of 1–30 s is performed and the results are stacked over the two periods of time, before and after impoundment of the Gotvand dam (Fig. 4). Data were selected according to signal-to-noise ratio (SNR) and interstation distance. According to previous studies, it is possible to obtain a sufficient randomization of the seismic noise to produce realistic Green's functions cross-correlating long time periods of data (e.g. more than 1 month, e.g. Shapiro et al. 2005; Bensen et al.2007). So the availability of 10 months continuous data, for both before and after impoundment, provides a suitable dataset to obtain coherent EGF with good SNR. For EGFs with SNR > 10 on both lags, we average the positive and negative lags to obtain the so called ‘symmetric component’ and therefore taking into account heterogeneities in the noise wavefield. For the remaining interstation measurements, we use only the lag with the higher SNR and retain the measurement. Fig. 4 shows cross-correlations of recorded signals for the time period before and after impoundment. Furthermore, to avoid spurious signals due to interference between the positive and negative parts of the correlation (Shapiro et al.2005; Lin et al.2007), we select only cross-correlation signals belonging to station pairs at a distance greater than two wavelengths. Figure 4. View largeDownload slide Stacks of daily cross-correlation pairs in the frequency bands of 1–30 s for before (a) and after (b) the impoundment with SNR > 10. They aligned depending on their interstation distance. Figure 4. View largeDownload slide Stacks of daily cross-correlation pairs in the frequency bands of 1–30 s for before (a) and after (b) the impoundment with SNR > 10. They aligned depending on their interstation distance. 3 GROUP VELOCITY MEASUREMENT AND SURFACE WAVE TOMOGRAPHY In order to construct group velocity distribution maps at 2- to 8-s periods using all the estimated dispersion measurements from Rayleigh wave EGFs, tomographic inversion technique of Ditmar & Yanovskaya (1987) and Yanovskaya & Ditmar (1990) was used. The methodology which has been applied in several papers (e.g. Ritzwoller & Levshin 1998; Guidarelli & Aoudia 2016) is a 2-D generalization of the classical 1-D method of Backus & Gilbert (1968). The study region was parametrized into grids with cell size of 0.05° × 0.05° and with proper regularization parameters to provide relatively smooth maps with small data misfits. The same regularizations parameters were used for producing the maps, before and after the impoundment. The resolution of the data set is controlled by the average path length, the density of paths and the azimuthal coverage. In this study, the spatial averaging kernels are estimated for each point (Backus & Gilbert 1968). In addition to group velocity maps, the corresponding resolution information will be provided at each period. The Yanovskaya's methodology (1997) is used to calculate the spatial resolution, which varies from 4 to about 10 km in the study region. This methodology allows calculating the spatial averaging kernels at each point and in different directions. A functional S(x, y) is calculated for different orientations of the coordinate system in order to determine sizes of the averaging area along different directions. The averaging area can be approximated by an ellipse, centred in a point of the studied region, with axes equal to the largest and the smallest values of S(x, y). The size of the averaging area along different directions are defined by S(x, y), which gives an idea of the resolution in each point (x, y). It is approximated by an ellipse with half-axes of the Smax(x, y) and Smin(x, y). The smallest [Smin(x, y)] and largest [Smax(x, y)] axes of the ellipse are calculated, and the resolution at each point is given by the mean size of the averaging area L = [Smin(x, y) + Smax(x, y)]/2. Applying the described tomographic method, the Rayleigh wave group velocity maps for periods ranging from 2 to 8 s are generated. Fig. 5 shows the group velocity distribution at 3, 5 and 8 s with the corresponding resolution maps and ray path coverage in the study region. Only the seismic velocities of the regions where the resolution is ≤10 km (Fig. 5) are taken into consideration for further discussion. Figure 5. View largeDownload slide Rayleigh wave group speed maps derived from ambient noise tomography at periods of 3, 5 and 8 s for before (a-c) and after (g-i) the impoundment. The resolution maps and ray path coverage at periods of 3, 5 and 8 s are presented for before (d-f) and after (j-l) the impoundment, respectively. Figure 5. View largeDownload slide Rayleigh wave group speed maps derived from ambient noise tomography at periods of 3, 5 and 8 s for before (a-c) and after (g-i) the impoundment. The resolution maps and ray path coverage at periods of 3, 5 and 8 s are presented for before (d-f) and after (j-l) the impoundment, respectively. The 3 s group velocity map before impoundment (Fig. 5) shows a distinct low velocity zone in the northern part of the reservoir, south of Lali, which extends as a NW-SE oriented narrow band up to Masjed Soleyman reservoir, located east of city with the same name. This trend of Low Velocity Zone (LVZ) is also visible in 5 and 8 s map, but with a little higher velocity (from 1.5 to 2.0 km s–1), respectively. Another low velocity zone is clear around the dam location at 3 and 5 s period map, but it vanishes at 8 s. A very narrow high velocity zone is in parallel present at the south of the LVZ. The area and velocity of this zone which is visible at all different periods, increases with period from 3 to 8 s. In the eastern part of the area, beneath the Masjed Soleyman reservoir a distinct high velocity zone is observed at all periods. After impoundment, distribution of the LVZ remains almost the same as before, mainly beneath the reservoir. A little increase is observed in velocity and area of the High Velocity Zone (HVZ), which is located around the LVZ, especially on southeastern margin of the study area at 5 and 8 s periods. Beneath, south and north of the Masjed Soleyman reservoir, a clear HVZ is observed at all periods, before and after the Gotvand reservoir impoundment. The impoundment of the Masjed Soleyman dam backs to 2001 about 5 yr before installation of the Gotvand and Masjed Soleyman seismological networks. 4 3-D SHEAR VELOCITY MODEL A 3-D shear wave velocity model was obtained extracting local group velocity dispersion curves for each node across the 0.05° × 0.05° gridpoints of the 2-D tomographic maps in the period range of 2–8 s. We adopt an iterative linearized least-square inversion scheme (Herrmann 2013) to invert the extracted dispersion curve for a 1-D shear wave profile. The linearized inversion technique shows some dependence on the initial model. The initial model is constructed using the available shear wave velocity models of the region obtained from 1-D inversion of the locally recorded earthquakes by the Gotvand seismic network (Tatar et al. 2009). Our initial model is a constant value two-layer velocity model composed of 10 thin layers (0–10 km) each one 1 km thick, overlying a half-space (black dashed line in Fig. 6b). The S-wave velocity in each layer is constant. In order to illustrate the inversion process, an example of Vs inversion for an extracted group velocity dispersion curve at a gridpoint, which its location is shown by a yellow star in Fig. 3, is presented in Fig. 6. Figure 6. View largeDownload slide (a) Predicted group velocity dispersion curve (grey lines) from the ensemble of Vs models shown in (b), and average of the ensemble dispersion curves (blue line). The black line is the observed (extracted) dispersion curve with error bars at individual periods. Result of inversion of the average of Vs models (red line in b) is presented as red line. (b) The resulting ensemble of acceptable Vs models (grey lines) and the average of the ensemble (red line). The initial model is presented as black dashed line, composed of 10 separated layers, each one 1 km thick. Figure 6. View largeDownload slide (a) Predicted group velocity dispersion curve (grey lines) from the ensemble of Vs models shown in (b), and average of the ensemble dispersion curves (blue line). The black line is the observed (extracted) dispersion curve with error bars at individual periods. Result of inversion of the average of Vs models (red line in b) is presented as red line. (b) The resulting ensemble of acceptable Vs models (grey lines) and the average of the ensemble (red line). The initial model is presented as black dashed line, composed of 10 separated layers, each one 1 km thick. To verify the resolution of the inversion result with depth, we plotted the quantity dU(T)/dVs (Fig. 7) which reflects the sensitivity of the fundamental mode Rayleigh-wave group velocity on the S-wave velocity versus depth for one examined cell located on profile AA’ (Fig. 3) for the periods of 2, 3, 5 and 8 s. This figure shows that Rayleigh group velocity at 2, 3 and 5 s are sensitive to variation of shear wave velocity at depths smaller than 2, 3 and 4 km, respectively. At 8 s, the shortest period of our data shows a peak at a depth of about 6 km. So S-wave velocities at depths greater than 6 km do not have a significant influence on the values of group velocities (very low absolute sensitivity), which indicate our results are valid until depth ∼6 km. Deeper anomalies are not reliable due to large errors and low resolution of results. Figure 7. View largeDownload slide Normalized sensitivity kernels for the fundamental mode Rayleigh-wave group velocity, calculated at different periods of 2, 3, 5 and 8 s for a grid located on profile AA’ (Fig. 3). Figure 7. View largeDownload slide Normalized sensitivity kernels for the fundamental mode Rayleigh-wave group velocity, calculated at different periods of 2, 3, 5 and 8 s for a grid located on profile AA’ (Fig. 3). In order to smooth out the Vs profiles and reduce the velocity contrasts over depths, during the inversion each Vs profile was vertically damped to minimize the differences at neighbouring layers. All constructed 1-D Vs profiles for each gridpoint then are assembled to form the final 3-D Vs model. Horizontal sections of our 3-D shear wave velocity model at depths of 3, 6 and 9 km for before and after impoundment are shown in Fig. 8. Figure 8. View largeDownload slide Shear wave velocity maps at depths of 3, 6 and 9 km before (a–c) and after (d–f) impoundment of the Gotvand dam. Figure 8. View largeDownload slide Shear wave velocity maps at depths of 3, 6 and 9 km before (a–c) and after (d–f) impoundment of the Gotvand dam. Vs model at different depths (Fig. 8) reveal almost the same features as the group velocity maps at different periods. Before impoundment of the reservoir, the Vs model at 3 km depth shows a LVZ just in north of the central part of the reservoir, south of Lali, where a distinct low velocity region is observed. A narrow HVZ with NW-SE orientation is visible in southern margin of the reservoir, especially at depth of 6 km. Comparison of the Vs models for before and after impoundment (Fig. 8) suggest a relatively velocity reduction beneath the reservoir at shallow depths. In order to estimate the uncertainty in our shear-wave velocity model we used a bootstrap analysis (Efron & Gong 1983; Efron & Tibshirani 1993); according to this methodology error estimation is based on the inversion of different subsets of data from all produced sets of solutions. In this analysis, we produced 1000 data sets, called bootstrap samples, from the full set of interstation group velocity measurements. Each set produced by randomly selecting 75 per cent of the individual group velocity measurements; then each set is processed by the 2-D tomographic inversion and 1-D shear wave velocity inversion. Therefore, the bootstrap analysis propagated errors through both the 2-D tomographic inversion and 1-D shear wave velocity inversion. Fig. 9 shows the standard deviation of the model ensemble for our shear wave velocity model at shallow depths (<6 km) vary in the range of ∼0.03–0.10 km s–1 for the inversion before (Figs 9a–c) and after (Figs 9d–f) impoundment. Generally larger standard deviation values (up to ∼0.15 km s–1) are found at greater depths (>6 km), which is possibly due to decreasing the number of available dispersion measurements at the longer periods and also relatively poor depth resolution of the Rayleigh waves for longer periods compare to short periods (Fig. 7). Figure 9. View largeDownload slide Uncertainty in the shear wave velocity estimates from 1000 bootstrap resampled data sets at 3, 6 and 9 km depths, before (a–c) and after (d–f) impoundment of the dam. Figure 9. View largeDownload slide Uncertainty in the shear wave velocity estimates from 1000 bootstrap resampled data sets at 3, 6 and 9 km depths, before (a–c) and after (d–f) impoundment of the dam. Two vertical cross-sections across the Gotvand reservoir (Fig. 3), plotted in Fig. 10, reveal the velocity variations in depth before and after impoundment. Corresponding standard deviation in these two vertical cross-sections are also plotted in Fig. 11, showing that uncertainty at shallow parts (2–6 km) almost vary in the range of ∼0.04–0.07 km s–1, which is much better than deeper parts (>6 km)that reaches up to ∼0.10 km s–1. The location of Profile BB’ is the same as geological map cross-section (Fig. 2). There is very good consistency between the geological features and the determined velocity structure. Figure 10. View largeDownload slide Vertical cross-sections of shear wave velocity for before (a and c) and after (b and d) impoundment of the Gotvand dam along two profiles delineated by the black lines in Figs 3 and 8(a). Topography is over plotted above individual cross sections. The location of the lake along the profile is presented as solid grey line above the topography. Figure 10. View largeDownload slide Vertical cross-sections of shear wave velocity for before (a and c) and after (b and d) impoundment of the Gotvand dam along two profiles delineated by the black lines in Figs 3 and 8(a). Topography is over plotted above individual cross sections. The location of the lake along the profile is presented as solid grey line above the topography. Figure 11. View largeDownload slide Uncertainty in shear wave estimated from 1000 bootstrap resampled data sets for before (a and c) and after (b and d) impoundment of the Gotvand dam along two profiles delineated by the black lines in Figs 3 and 8(a). Topography is over plotted above individual cross sections. The location of the lake along the profile is presented as solid grey line above the topography. Figure 11. View largeDownload slide Uncertainty in shear wave estimated from 1000 bootstrap resampled data sets for before (a and c) and after (b and d) impoundment of the Gotvand dam along two profiles delineated by the black lines in Figs 3 and 8(a). Topography is over plotted above individual cross sections. The location of the lake along the profile is presented as solid grey line above the topography. As it is inferred from the cross-section view of the geological map (Fig. 2), Gachsaran is the most dominant formation underneath the reservoir. This sedimentary formation caused by an evaporate environment at the end of the Upper Miocene. As faulting, folding, and flow are the controlling factors for the thickness variations in the Gachsaran formation, so strong differences in seismic images can be observed due to variations in the lithology and thickness caused by an uplifted rock or faulting along this formation (Alavi 1994). The Gachsaran Formation consists of the high velocity anhydrite, salt, and limestone, and low velocity rocks of marl and shale (Abdollahie Fard et al. 2011). For the time period before the impoundment of the dam, a very low velocity zone is observed in the northern part of the reservoir along the lines AA’ (Fig. 10a) and BB’ (Fig. 10c). Comparing with the geological cross section along the BB’ profile, this low velocity zone is related to the Gachsaran formation which is the dominant formation in the area located in northern parts of the reservoir. This gipsy formation that made mainly from low velocity evaporate deposits is spread over a wide area at surface and in depth. The shape, size and S-wave velocity (∼2 km s–1) of observed anomaly north of the lake on BB’ shear wave velocity profile (Fig. 10c) matches perfectly with the specification of Gachsaran formation on geological cross section, which is located north of the lake (Fig. 2). The high velocity zone observed southwest of the low velocity zone of Gachsaran formation, at a depth of ∼2–8 km (Figs 10c and d) shows a very good consistency with the folded and uplifted limestone formations (Fig. 2) of the Asmari-Pabedeh (Vs ∼ 2800 m s–1), Ilam-Sarvak and Kazhdum (Vs ∼ 3300 m s–1). These cross sections (A and C on Fig. 10), reveal that the ANT technique is a power tool with enough good resolution for detecting lithologically different formations even at small scale. It means that very dense small-spacing seismological networks can be used for detecting of hydrocarbon traps applying ANT technique. The effect of the reservoir impoundment is observed in Figs 10(b) and (d). It seems that both low and high velocity anomalies related to the Gachsaran and Asmari, Pabdeh, Ilam-Sarvak group respectively, have been compacted due to loading of the above reservoir water column. A very distinct high velocity anomaly is observed just south of the Gotvand lake (Fig. 10d). Based on field investigation, this high velocity anomaly can be related to the existence of a salt dome which is responsible for very high salinity level of the reservoir water. Due to misestimating of the size and capacity of this salt dome which its northern flank has direct contact with the reservoir water, the reservoir water that became hypersaline (very salty) is almost unusable for agricultural purposes. 5 DISCUSSION In order to better define the impoundment effects, the difference between the shear-wave velocity after and before the impoundment of the Gotvand reservoir was calculated along the AA’ and BB’ cross section (Fig. 12). This Figure shows clearly that along both cross sections we observe a decrease in S-wave velocity (red colour) up to ∼0.3 ± 0.05 km s–1 at shallow depths (2–4 km) just beneath the lake (reservoir). This can be due to the seepage of water at shallow depth, which causes reduction of shear wave velocity. But at deeper parts (4–6 km) beneath the reservoir, an increase of ∼0.3 ± 0.09 km s–1 (blue colour) is observed in S-wave velocity that can be due to compaction effects of the reservoir water. However, due to the large uncertainty in S-wave velocity based on computed resolution kernel curves (Fig. 7) large velocity anomalies at deeper parts are not enough reliable to be taken into consideration for further interpretation. Figure 12. View largeDownload slide Difference in shear wave velocity observed after and before the impoundment (Vs after –Vs before) for drawn vertical cross-sections along AA’ (a) and BB’ (b) profiles (Fig. 8). Solid grey lines above the topography show the location of the lake along the profile. Figure 12. View largeDownload slide Difference in shear wave velocity observed after and before the impoundment (Vs after –Vs before) for drawn vertical cross-sections along AA’ (a) and BB’ (b) profiles (Fig. 8). Solid grey lines above the topography show the location of the lake along the profile. In order to better understand the variations of velocities before and after impoundment due to the effect of water load on the crust beneath the dam, we modelled the incremental strength underneath the dam along the profiles shown in Fig. 3. According to coulombs law, this incremental strength due to reservoir impoundment is given by (Bell & Nur 1978):   $$\left\{ \begin{array}{@{}l@{}} \Delta S = \mu \left( {\Delta {\sigma _n} - \Delta P} \right){\rm{ }} - \Delta {\tau _n},\ \ \ \ \ \ \ \ \\ \Delta P = \Delta {P_u} + \Delta {P_{diff}} \end{array} \right.$$ (1)  $$\Delta {P_u} = B\Delta {\sigma _m}$$ (2)where Δσn and Δτn are the incremental normal and shear stresses on the fault plane due to the water loading, respectively, μ is the coefficient of friction, B is the Skempton's coefficient, Δσm is the mean stress, which is obtained from summation over the diagonal elements of the stress tensor, ΔPu is the undrained pore pressure, occurs instantly in response to undrained loading of water, and after a while pore pressure, ΔPdiff, is diffused from the reservoir into the crust beneath it. The change in strength due to loading can be calculated by knowing the filling history of the dam, location, and geometry of the fault. Then, the Δσn, Δτn, ΔPu and ΔS can be found from the eqs (1) and (2). In order to calculate the strength beneath the AA’ and BB’ lines, we assumed a vertical fault beneath them, while its western block move upward. The normal and the shear stresses can be given by (Jaeger & Cook 1979):   $$\left\{ \begin{array}{@{}l@{}} \Delta {\sigma _n} = l\Delta {\sigma _x}(T) + m\Delta {\sigma _y}(T) + n\Delta {\sigma _z}(T),\\ \Delta {\tau _n} = \pm \sqrt {(\Delta {\sigma _x}{{(T)}^2} + \Delta {\sigma _y}{{(T)}^2} + \Delta {\sigma _z}{{(T)}^2}) - \Delta \sigma _n^2} , \end{array} \right.$$ (3)where l, m, n are the direction cosines of the normal to the fault plane and Δτn is positive when the shear stress favours faulting. We also considered B = 0.7 (Talwani 1997) and μ = 0.75 (Byerlee1978). A negative value for ΔS may imply a weakening or destabilization and a positive value means strengthening or stabilization (Bell & Nur 1978) in the sense of causing the stress state to move toward or away from failure. Load of a reservoir can cause compressive stresses underneath the reservoir, which can lead to stabilizing the underlying crust and produce tensile stresses in horizontal direction on the surrounding area of the reservoir (Liu & Xu 2005). These tensile stresses can create small shallow fractures and facilitate the penetration of water into deeper rocks and change the properties of the crust (Liu et al. 2011). Assuming a homogeneous medium, Boussinesq solutions (Jaeger & Cook 1979) were used to calculate the elastic stress changes due to reservoir loading. The surface of the reservoir is divided into small grid blocks with 100 m × 100 m dimension. Summing all the stress changes in each grid block the total stress change at each point underneath the reservoir is calculated. The final results for incremental strength at total volume of impoundment 2.1 × 109 m3, and lake level of 197 m (Fig. 13) indicate that all parts underneath the lake have a positive value which means stabilization or strengthening. Away from the lake, the results shows negative values which leads to destabilization or weakening. Figure 13. View largeDownload slide Strength changes after impoundment of the Gotvand lake (total volume of impoundment ∼2.1 × 109 m3, and lake level of 197 m) along (a) Profile AA’ and (b) Profile BB’. Solid grey line shows the location of the lake. Figure 13. View largeDownload slide Strength changes after impoundment of the Gotvand lake (total volume of impoundment ∼2.1 × 109 m3, and lake level of 197 m) along (a) Profile AA’ and (b) Profile BB’. Solid grey line shows the location of the lake. Both elastic loading stress and pore pressure diffusion (ΔPdiff) can affect the area underneath the reservoir (Jaeger & Cook 1979). As the largest volume of the water is behind the crest of the Gotvand dam, so it is expected that the effect of loading stress being more effective in this part. This effect can be seen in Fig. 13. Beneath the profile AA’ which crosses the deepest part of the reservoir the amount of incremental strength is higher than what we observed beneath the profile BB’. Compressive stress due to loading the dam can change lithology of the rocks and enhance small fractures in the periphery of the reservoir and increase the pore pressure. Estimated analytical strength due to loading, leads to stabilizing the crust underneath the reservoir, which is in agreement with earlier reports (Roeloffs 1988; Liu & Xu 2005; Santoyo et al. 2010). Shear wave velocity difference for before and after impoundment (Fig. 12) indicates on considerable decrease of Vs (red colour) up to ∼0.4 ± 0.1 km s–1 at very deep parts (depths greater than 8 and 6 km on profile AA’ and BB’, respectively) beneath and outside of the lake where the compaction effects (analytical strength) is less than 0.2 bar. It seems that deep parts of the crust underneath the reservoir (5–8 km depth) is stabilized due to the compressional stress loading, while outside of the reservoir we observe weakening due to extend of the microcracking and fractures as a result of horizontal tensile stress (Liu et al. 2011). It should be taken into consideration that the presence of incompetent lithology in some parts can facilitate the permeability and migration of water in to deeper parts (Talwani et al. 2007). Finally, the loading effect of the Gotvand reservoir on underneath strength changes is calculated along the profile AA’ by incorporating the effect of pore pressure (ΔPdiff) in lowering the strength of rocks (Fig. 14a). Figure 14. View largeDownload slide Vertical cross-section 3-D view of (a) distribution of pore pressure and (b) strength changes along the profile AA’ after the impoundment of the lake (total volume of impoundment ∼2.1 × 109 m3, and lake level of 197 m) by adding the pore pressure effect (coordinate system in UTM). Figure 14. View largeDownload slide Vertical cross-section 3-D view of (a) distribution of pore pressure and (b) strength changes along the profile AA’ after the impoundment of the lake (total volume of impoundment ∼2.1 × 109 m3, and lake level of 197 m) by adding the pore pressure effect (coordinate system in UTM). Diffused pore pressure (ΔPdiff) can be calculated by solving the following equation at a time t and depth z:   \begin{eqnarray} \Delta {P_{{\rm{diff}}}}(z,t) &=& (1 - \alpha ){p_0}\,{\rm{erfc}}\,\left[ {\frac{z}{{{{(4ct)}^{1/2}}}}} \right]\, + \,\alpha H(t){p_0}\\ \alpha &=& \,B(1 + {v_u})/[3(1 - {v_u})], \end{eqnarray} (4)where the erfc is the complementary error function, H(t) is the Heaviside unit step function, B is the Skempton's coefficient, vu is the undrained Poisson's ratio, z is the depth beneath the reservoir, p0 is the pressure at the bottom of the reservoir, c is the hydraulic diffusivity and t is time. The hydraulic diffusivity constant in eq. (4) was estimated to be c = 0.1 m2 s–1, using, c = r2/4Δt, where Δt is the time lag between the start of the impoundment and the time of estimation and r is the hypocentral distance from the dam (Talwani et al. 2007). Fig. 14(b) shows that due to penetration of water or pore pressure diffusion, which affects shallow depths, decrease in rock strength value is expected. However, in deeper parts, due only to loading effect, an increase in rock strength is observed. These observations are in good agreement with shear wave velocity decrease at shallow depths and its increase at deeper parts beneath the reservoir as it can be seen in Fig. 12 after impoundment of the dam. 6 CONCLUSION We obtained in this study a detailed image of the seismic velocity structures in the top ∼6 km of the region beneath the second largest dam reservoir in south of Iran using Rayleigh wave group velocity measurements from cross-correlation of ambient noise data in the period range of 2–8 s. Group velocity maps at 3, 5 and 8 s periods are obtained by surface wave tomography. The observed velocity structures primarily correlate well with geological features. 3-D shear wave velocity model was obtained by constructing a local group velocity dispersion curves for each node of the 0.05° × 0.05° grid of the 2-D tomographic maps in the period range of 2–8 s. The obtained images for shear wave velocities reveal clear velocity contrasts across the study region. A very low velocity zone is observed beneath the northern part of the reservoir related to the Gachsaran formation which is the dominant formation in the area located in northern parts of the reservoir. The high velocity zone observed southwest of the low velocity zone of Gachsaran formation, at a depth of ∼4–8 km consistent with the folded and uplifted limestone formations of the Asmari-Pabedeh, Ilam-Sarvak and Kazhdum. Calculated difference between the shear wave velocity after and before the impoundment of the Gotvand reservoir shows a clear decrease in S-wave velocity up to ∼0.3 ± 0.05 km s–1 at shallow depths (2–4 km) just beneath the lake (reservoir) which can be related to the seepage of water in shallow depths. At deeper parts (4–6 km) beneath the reservoir, an increase of ∼0.3 ± 0.1 km s–1 is observed in S-wave velocity that can be associated to compaction due to loading effects of the reservoir water. Calculation of loading effects of the Gotvand reservoir on underneath strength changes, considering the effect of pore pressure, show good agreement with the observed shear-wave velocity changes across the reservoir. The results of this research show that building a dam can make severe changes in the crust properties underneath it. Considering the lithology, volume of the water and location of the dam, these changes could be different at distinct distances and depths. Acknowledgements This research was supported and has been done as part of the contract between International Institute of earthquake Engineering and Seismology (IIEES), the Iranian Water and Power Resources Development Co. (IWPCO)-Upper Gotvand Project and the Khuzestan Water and Power Authority (KWPA) for seismic monitoring of the Gotvand and Masjed Soleyman reservoirs. The data used in this work were recorded by the Gotvand-Masjed Soleyamn permanent Seismic Network operated by IIEES. We are grateful to the colleagues of IIEES, IWPCO-Upper Gotvand Projects, and KWPA Masjed Solymand dam for helps and supports during installation, gathering, and preparing of the data. Special thanks to Dr Dariush Mahjoub from IWPCO and Dr Yaghob Arab from KWPA for their continuous supports. We thank Regione Friuli-Venezia Giulia (Cooperazione Internazionale e allo Sviluppo) for support through the project “Sostegno all Éducazione e alla Ricerca nellámbito della prevenzione del rischio sismico in Iran”. We gratefully thank and acknowledge two anonymous reviewers for constructive reviews. Some figures were made using the public domain Generic Mapping Tools (Wessel & Smith 1998). REFERENCES Abdollahie Fard I., Sepehr S., Sherkati S. 2011. Neogene salt in SW Iran and its interaction with Zagros folding, Geol. Mag. , 148, 854– 867. Google Scholar CrossRef Search ADS   Alavi M., 1994. Tectonics of the zagros orogenic belt of iran; new data and interpretations, Tectonophysics , 229, 211– 238. Google Scholar CrossRef Search ADS   Backus G., Gilbert F., 1968. The resolving power of gross Earth data, Geophys. J. R. astr. Soc. , 16, 169– 205. Google Scholar CrossRef Search ADS   Bell M.L., Nur A., 1978. Strength changes due to reservoir-induced pore pressure and stress and application to Lake Oroville, J. geophys. Res. , 83, 4469– 4483. Google Scholar CrossRef Search ADS   Bensen G.D., Ritzwoller M.H., Barmin M.P., 2007. Processing seismic ambient noise data to obtain reliable broad-band surface wave dispersion measurements, Geophys. J. Int. , 169, 1239– 1260. Google Scholar CrossRef Search ADS   Berberian M., 1995. Master blind thrust faults hidden under the Zagros folds: active basement tectonics and surface morphotectonics, Tectonophysics , 241, 193– 224. Google Scholar CrossRef Search ADS   Brenguier F., Shapiro N.M., Campillo M., Nercessian A., Ferrazzini V., 2007. 3D surface wave tomography of the Piton de la Fournaise volcano using seismic noise correlations, Geophys. Res. Lett. , 34, L02305, 1– 5. Google Scholar CrossRef Search ADS   Byerlee J., 1978. Friction of rocks, Pure appl. Geophys. , 122, 947– 965. Ditmar P.G., Yanovskaya T.B., 1987. Generalization of the Backus-Gilbert method for estimation of the horizontal variations of surface wave velocity, Izv. Acad. Sci. USSR Phys. Solid Earth. (English Transl.) , 23, 470– 477. Droznin D.V., Shapiro N.M., Droznina S.Y.A., Senyukov S.L., Chebrov V.N., Gordeev E.I., 2015. Detecting and locating volcanic tremors on the Klyuchevskoy group of volcanoes (Kamchatka) based on correlations of continuous seismic records, Geophys. J. Int. , 203, 1001– 1010. Google Scholar CrossRef Search ADS   Dunnington H.V., 1968. Salt-tectonic features of northern Iraq. In: MATTOX, R. B. (ed.) A symposium based on paper from the International Conference on Saline Deposits, Geol. Soc. Am., Special Papers , 88, 183– 227. Google Scholar CrossRef Search ADS   Efron B., Gong G., 1983. A leisurely look at the Bootstrap, the Jackknife, and cross validation, Am. Stat. , 37, 36– 48. Efron B., Tibshirani R.J., 1993. An Introduction to the Bootstrap , Chapman and Hall/CRC. Google Scholar CrossRef Search ADS   Guidarelli M., Aoudia A., 2016. Ambient noise tomography of the Cameroon Volcanic Line and Northern Congo craton: new constraints on the structure of the lithosphere, Geophys. J. Int. , 204, 1756– 1765. Google Scholar CrossRef Search ADS   Guo Z., Gao X., Yao H.J., Li J., Wang W.M., 2009. Midcrustal low-velocity layer beneath the central Himalaya and southern Tibet revealed by ambient noise array tomography, Geochem. Geophys. Geosyst. , 10( 5), Q05007, 1–12. Google Scholar CrossRef Search ADS   Herrmann R.B., 2013. Computer programs in seismology: an evolving tool for instruction and research, Seism. Res. Lett. , 84, 1081– 1088. Google Scholar CrossRef Search ADS   Jaeger J.C., Cook N.G.W., 1979. Fundamentals of Rock Mechanics , 3rd edn. Chapman and Hall. James G.A., Wynd J.G., 1965. Stratigraphic nomenclature of Iranian oil consortium agreement area, AAPG Bull ., 49, 2182– 2245. Hessami Azar K., Jamali F., Tabasi H., 2003. Major Active Faults of Iran , IIEES Publication. Koulakov I.et al.  , 2014. Asymmetric caldera-related structures in the area of the Avacha group of volcanoes in Kamchatka as revealed by ambient noise tomography and deep seismic sounding, J. Volc. Geotherm. Res. , 285, 36– 46. Google Scholar CrossRef Search ADS   Lin F.C., Ritzwollerl M.H., Townend J., Bannister S., Savage M.K., 2007. Ambient noise Rayleigh wave tomography of new Zealand, Geophys. J. Int. , 170, 649– 666. Google Scholar CrossRef Search ADS   Liu S.M., Xu L.H., 2005. Finite-element simulation of hydraulic-pressure stress field in Danjiangkou Reservoir area, J. Hydraul. Eng. , 36, 863– 869. Liu S.M., Xu L.H., Talwani P., 2011. Reservoir-induced seismicity in the Danjiangkou Reservoir: a quantitative analysis, Geophys. J. Int. , 185, 514– 528. Google Scholar CrossRef Search ADS   Lobkis O., Weaver R.L., 2001. On the emergence of the Green's function in the correlations of a diffuse field, J. acoust. Soc. Am. , 110( 6), 3011– 3017. Google Scholar CrossRef Search ADS   Mcquarrie N., 2004. Crustal scale geometry of the zagros fold–thrust belt, Iran, J. Struct. Geol. , 26, 519– 535. Google Scholar CrossRef Search ADS   Mouthereau F., Lacombe O., Meyer B., 2006. The Zagros folded belt (Fars, Iran): constraints from topography and critical wedge modelling, Geophys. J Int. , 165, 336– 356. Google Scholar CrossRef Search ADS   Nissen E., Tatar M., Jackson J.A., Allen M.B., 2011. New views on earthquake faulting in the Zagros fold-and-thrust belt of Iran, Geophys. J. Int. , 186( 3), 928– 944. Google Scholar CrossRef Search ADS   O’brien C.A.E., 1950. Tectonic problems of the oil field belt of southwest Iran, Proceedings of 18th International Geological Congress, part 6 , Great Britain, pp. 45– 58. Ritzwoller M.H., Levshin A.L., 1998. Eurasian surface wave tomography: group velocity, J. geophys. Res. , 103, 4839– 4878. Google Scholar CrossRef Search ADS   Roeloffs E.A., 1988. Fault stability changes induced beneath a reservoir with cyclic variations in water level, J. geophys. Res. , 93, 2107– 2124. Google Scholar CrossRef Search ADS   Santoyo M., García-Jerez A., Luzón F., 2010. A subsurface stress analysis and its possible relation with seismicity near the Itoiz Reservoir, Navarra, Northern Spain, Tectonophysics , 482, 205– 215. Google Scholar CrossRef Search ADS   Sattarzadeh Y. 1997. Active tectonics in the Zagros Mountains, Iran , PhD thesis, Imperial College, London. Setudehnia A., 1972. Iran du Sud-Ouest: Lexique Strat. Internat, Centre Nat. Rech. Scientifique , Paris, III, Asie, Fasc. 9b, p. 289– 376. Shapiro N.M., Campillo M., Stehly L., Ritzwoller M.H., 2005. High-resolution surface wave tomography from ambient seismic noise, Science , 307, 1615– 1618. Google Scholar CrossRef Search ADS PubMed  Snieder R., 2004. Extracting the Green's function from the correlation of coda waves: a derivation based on stationary phase, Phys. Rev. E. , 69, 046610. Google Scholar CrossRef Search ADS   Talwani P., 1997. On the nature of reservoir-induced seismicity, Pure appl. Geophys. , 150, 473– 492. Google Scholar CrossRef Search ADS   Talwani P., Acree S., 1984/1985. Pore pressure diffusion and the mechanism of reservoir-induced seismicity, Pure appl. Geophys. , 122, 947– 965. Google Scholar CrossRef Search ADS   Talwani P., Chen L., Gahalaut K., 2007. Seismogenic permeability, ks, J. geophys. Res. , 112, B07309, 1–18. Google Scholar CrossRef Search ADS   Tatar M., Yamini Fard F., Dezvareh M. 2009. Seismicity and seismotectonic study of the Gotvand-e Olya dam region using local earthquakes, annual report, IIEES, pp. 120. Wessel P., Smith W.H.F., 1998. New, improved version of the Generic Mapping Tools released, EOS, Trans. Am. geophys. Un. , 79, 579. Google Scholar CrossRef Search ADS   Yanovskaya T. B., 1997. Resolution estimation in the problems of seismic ray tomography, Izv. Phys. Solid Earth , 33( 9), 762– 765. Yanovskaya T.B., Ditmar P.G., 1990. Smoothness criteria in surface wave tomography, Geophys. J. Int. , 102, 63– 72. Google Scholar CrossRef Search ADS   Yao H., Van der Hilst R.D., De Hoop M.V., 2006. Surface-wave tomography in SE Tibet from ambient seismic noise and two-station analysis—I. Phase velocity maps, Geophys. J. Int. , 166, 732– 744. Google Scholar CrossRef Search ADS   © The Authors 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Geophysical Journal International Oxford University Press

Loading effects beneath the Gotvand-e Olya Reservoir (south-west of Iran) deduced from ambient noise tomography

, Volume 212 (1) – Jan 1, 2018
15 pages

Publisher
Oxford University Press
ISSN
0956-540X
eISSN
1365-246X
D.O.I.
10.1093/gji/ggx411
Publisher site
See Article on Publisher Site

Abstract

Abstract In order to define the precise shallow velocity structure beneath the second largest dam reservoir in Iran and to understand the loading effects on the underlying crust, the shear wave velocity of the shallow structure beneath the Gotvand-e Olya (hereinafter referred to as Gotvand) reservoir is determined through the inversion of group velocities obtained from seismic ambient noise tomography, using continuous data from 10 stations of a local network, installed to monitor the induced seismicity in the region surrounding the Gotvand and Masjed Soleyman dams for potential hazard. We obtained Rayleigh waves from cross-correlation of waveforms recorded 10 months before and the same duration after impoundment of the Gotvand reservoir and calculated the group velocity from dispersion analysis in the period range 2–8 s. The group velocity dispersion curves are used to produce 2-D group velocity tomographic maps. The resulting tomographic maps at short periods are well correlated with subsurface geological features and delineate distinct low- and high-velocity zones separated mainly by geological boundaries. The 3-D shear wave velocity structure provides detailed information about the crustal features underneath the reservoir. The results are consistent with the lithology of the region, and attest that ambient noise tomography (ANT) can be used for detailed studies of the velocity structure and lithology at shallow depths using continuous data from a dense local seismic network. An increase of shear wave velocity is observed at the deep parts (4–6 km) underneath the reservoir after impoundment of the dam, which could be caused by the changes in rocks properties after impoundment. However, at shallow depths (2–4 km), a decrease of Vs velocity is observed that can be associated to the penetration of water after the impoundment. Changes of material strength underneath the Gotvand reservoir, calculated using loading effects of the reservoir by considering the pore pressure are in good agreement with the observed shear wave velocity changes across the reservoir. Tomography14, Crustal imaging, Induced seismicity, Seismic noise 1 INTRODUCTION Building a dam causes severe stress load on the crust underneath the reservoir when the valley behind it is filled. The pore pressure of penetrated fluids can also change the physical properties of the Earth's crust via lowering the strength of rocks (Talwani & Acree 1985; Talwani 1997). Gotvand Dam with a total height of 180 m is the highest rock-fill dam with central clay core in Iran. The Gotvand reservoir with a total volume of 4500 Milion Cubic Meters and an area of ∼96.5 km2 at the level of 243 m from sea level is the second large reservoir in Iran after the Karkhe dam. Gotvand dam is located in the Zagros Mountains of western Iran, which is known as one of the most seismically active regions along the Alp-Himalayan belt. The Gotvand dam is located in the Zagros simply fold and thrust belt (SFB). Shortening at the surface in this area is accommodated by parallel trains of anticlines and synclines, while deformation in depth is modulated by parallel NW-SE high angle basement reverse faults and roughly N–S, right-lateral surface faults like Kazerun, Borazjan, Karebas, Sabzpushan and Sarvestan (Fig. 1). The N-dipping basement reverse faults have caused significant changes in the stratigraphic level across certain folds in the study region (Berberian 1995). Lahbari known also as Mountain Front Fault (MFF), which marks the southwestern limit of the surface exposure of the Asmari limestone is one of the most important N-dipping faults in the study region. Figure 1. View largeDownload slide Map of the Zagros Fold and Thrust belt with SRTM digital topography coloured to emphasize the high mountains (green) and the deep valleys (yellow) that penetrate the range. Seismicity (IIEES catalogue ML > 4.5) is shown by red circles. The white rectangle is the location of the study region. Faults from Hessami Azar et al. (2003). MFF, Mountain Front Fault; ZFF, Zagros Foredeep Faults; MRF, Main Recent Fault; BZ, Borazjan; KZ, Kazerun; KB, Karebas; SB, Sabzpushan. Figure 1. View largeDownload slide Map of the Zagros Fold and Thrust belt with SRTM digital topography coloured to emphasize the high mountains (green) and the deep valleys (yellow) that penetrate the range. Seismicity (IIEES catalogue ML > 4.5) is shown by red circles. The white rectangle is the location of the study region. Faults from Hessami Azar et al. (2003). MFF, Mountain Front Fault; ZFF, Zagros Foredeep Faults; MRF, Main Recent Fault; BZ, Borazjan; KZ, Kazerun; KB, Karebas; SB, Sabzpushan. Topographic analysis and seismotectonic studies show that activation of blind thrust faults in the basement have resulted in the evolution of asymmetric folds in the sedimentary cover (Sattarzadeh 1997; McQuarrie 2004; Mouthereau et al. 2006). Microseismicity of the region is dominated by reverse faulting, mainly confined at depths between 5 and 20 km. However the distribution of teleseismic centroid depths shows that most of the seismic strain is released at depths of 5–10 km (Nissen et al. 2011) mainly within the lower sedimentary cover. O’Brien (1950) and Dunnington (1968) compiled the Zagros stratigraphic column by using regional stratigraphy. They divided it into five distinct units, which from bottom to top, are the Precambrian metamorphic Basement Group; the Lower Mobile Group, which is considered as a decollement level and is formed by Late Proterozoic–Early Cambrian Hormuz salt, the Competent Group, which consists of Cambrian to Lower Miocene Formation; the Upper Mobile Group (Gachsaran Miocene salt); the Passive Group (mostly clastic sediments) at the top. In the study region (Fig. 2), Gachsaran evaporate is the main formation in the Upper Mobile Group. James & Wynd (1965) and Setudehnia (1972) described this formation as composed by seven distinct sections, consisting of anhydrite with bituminous shale at the bottom, salt and limestone at the middle and alternating marl at the top. Gachsaran Formation is incompetent and it is highly subject to extreme mobility as a result of differential pressure. Agha Jari and Bakhtyari formations are the main clastic deposits in the Passive Group at the top (Fig. 2). The principal lithology of Agha Jari formation consists of gypsum veins, calcareous sandstone, red marl and siltstone and the Bakhtiari Formation is terrigenous, clastic unit ranging in grain size from a silt grade to boulder conglomerate (Setudehnia 1972). Figure 2. View largeDownload slide Geologic map and cross-section view of the Gotvand and Masjed Soleyman Reservoirs (National Geoscience Database of Iran). Figure 2. View largeDownload slide Geologic map and cross-section view of the Gotvand and Masjed Soleyman Reservoirs (National Geoscience Database of Iran). In the last decade, the methodology called ambient noise tomography (ANT) has been developed to produce detailed seismic velocity models for the crust and upper mantle (Shapiro et al. 2005; Bensen et al. 2007). Cross-correlation of ambient seismic noise can produce short period surface waves between station pairs and thus can provide more accurate and detailed constraints on upper crustal structures than traditional surface wave tomography. It provides much higher resolution without need of using any earthquake surface-wave data (Yao et al. 2006; Guo et al. 2009). ANT is based on the estimation of empirical Green’s functions (EGF) between two receivers by cross-correlating the continuous ambient seismic noise recorded at the two receivers (Lobkis &Weaver 2001; Snieder 2004). Recently ambient noise tomography has been widely used from regional scale (Shapiro et al. 2005; Yao et al. 2006) to local scale such as volcanic regions (Brenguier et al. 2007; Koulakov et al. 2014; Droznin et al. 2015). This paper presents the results of ambient noise analysis of data collected about 10 months before and after the impoundment of the Gotvand dam by 10 seismological stations belonging to the dam seismic monitoring system. Gotvand dam has been continuously monitored by a local seismological network (Fig. 3) since 2006, about 5 yr before impoundment of the reservoir. However, the continuous seismic data was accessible only from 2011 March 20, about 10 months before starting the impoundment. The same duration of 10 months data set recorded after the impoundment and this until 2012 November 20 was used. The Gotvand reservoir impoundment actually was started on 2012 January 20 when a considerable increase in water level is observed. The impoundment of the Masjed Soleyman dam is started on 2001 December 12 and finished on 2002 June 25. There was no independent permanent seismological network monitoring the Masjed Soleyman reservoir prior to its impoundment. Figure 3. View largeDownload slide Topographic map of the study area including Gotvand and Masjed Soleyman dams reservoirs. Seismological stations of two Gotvand and Masjed Soleyman seismic networks are presented in black triangles. Cities are shown in black squares. Distribution of the reservoir lake at 243 m above the sea level is presented in blue. The black solid lines represent the location of vertical cross-sections. Yellow star represent the location of the gridpoint used to calculate sensitivity kernel (Fig. 7) for the fundamental Rayleigh wave group velocity. Figure 3. View largeDownload slide Topographic map of the study area including Gotvand and Masjed Soleyman dams reservoirs. Seismological stations of two Gotvand and Masjed Soleyman seismic networks are presented in black triangles. Cities are shown in black squares. Distribution of the reservoir lake at 243 m above the sea level is presented in blue. The black solid lines represent the location of vertical cross-sections. Yellow star represent the location of the gridpoint used to calculate sensitivity kernel (Fig. 7) for the fundamental Rayleigh wave group velocity. The purpose of this study is to determine the shallow structure beneath the Gotvand reservoir before and after impoundment of the dam and to characterize the possible changes in the shallow crust due to loading effect. We measure the group velocity of the cross-correlation Rayleigh waves at periods from 2 to 8 s. Then a surface wave tomography is performed to obtain high-resolution group velocity maps for the study area beneath the Gotvand reservoir. Group velocities from these maps will be inverted to produce high resolution 3-D model of shear wave velocity, which can provide important constraints for better understanding the loading effects due to reservoir impoundment on shallow structure characterization beneath the reservoir. 2 DATA PROCESSING We used 20 months of continuous vertical-component seismic data recorded by the Gotvand-Masjed Soleyman seismological network (Fig. 3) between 2011 March 20 and 2012 November 20. This seismic network consists of 10 medium band Trilium 40 seismometers equipped with 24-bit Taurus recorders both manufactured by Nanometrics company. All stations operate in continuous mode with a sampling rate of 80 sps. This seismological network belongs to Gotvand project of Iran Water and Power Resource Development Co. (IWPCO) and to Masjed Solyman dam of Khuzestan Water and Power Authority (KWPA). It has been designed and installed by International Institute of Earthquake Engineering and Seismology (IIEES) and is operated by IIEES since 2006. About 20 months of data, collected during 10 months before and after impoundment of the Gotvand dam are used to image shallow velocity structure of the area beneath the Gotvand reservoir from ambient noise tomography. The same procedure as described in Bensen et al. (2007) is applied for ambient noise data processing. Using only vertical components of ambient noise, we could only extract Rayleigh waves from waveform cross-correlations. First, instrument responses, mean and trend are removed from the continuous data that are then bandpass filtered in the period range 1–30 s. The contamination of ambient noise recordings by earthquakes and non-stationary noise sources near the stations is partly removed by applying the procedure of ‘temporal normalization’ which is called running-absolute-mean normalization (Bensen et al.2007). This normalization method computes the running average of absolute value of the seismic waveform in a normalization time window of fixed length and weights at the centre of the window by the inverse of this average (Bensen et al. 2007). This is followed by spectral whitening which is applied to flatten spectra over the entire period band (1–30 s). Then daily cross-correlations of continuous noise data in the period range of 1–30 s is performed and the results are stacked over the two periods of time, before and after impoundment of the Gotvand dam (Fig. 4). Data were selected according to signal-to-noise ratio (SNR) and interstation distance. According to previous studies, it is possible to obtain a sufficient randomization of the seismic noise to produce realistic Green's functions cross-correlating long time periods of data (e.g. more than 1 month, e.g. Shapiro et al. 2005; Bensen et al.2007). So the availability of 10 months continuous data, for both before and after impoundment, provides a suitable dataset to obtain coherent EGF with good SNR. For EGFs with SNR > 10 on both lags, we average the positive and negative lags to obtain the so called ‘symmetric component’ and therefore taking into account heterogeneities in the noise wavefield. For the remaining interstation measurements, we use only the lag with the higher SNR and retain the measurement. Fig. 4 shows cross-correlations of recorded signals for the time period before and after impoundment. Furthermore, to avoid spurious signals due to interference between the positive and negative parts of the correlation (Shapiro et al.2005; Lin et al.2007), we select only cross-correlation signals belonging to station pairs at a distance greater than two wavelengths. Figure 4. View largeDownload slide Stacks of daily cross-correlation pairs in the frequency bands of 1–30 s for before (a) and after (b) the impoundment with SNR > 10. They aligned depending on their interstation distance. Figure 4. View largeDownload slide Stacks of daily cross-correlation pairs in the frequency bands of 1–30 s for before (a) and after (b) the impoundment with SNR > 10. They aligned depending on their interstation distance. 3 GROUP VELOCITY MEASUREMENT AND SURFACE WAVE TOMOGRAPHY In order to construct group velocity distribution maps at 2- to 8-s periods using all the estimated dispersion measurements from Rayleigh wave EGFs, tomographic inversion technique of Ditmar & Yanovskaya (1987) and Yanovskaya & Ditmar (1990) was used. The methodology which has been applied in several papers (e.g. Ritzwoller & Levshin 1998; Guidarelli & Aoudia 2016) is a 2-D generalization of the classical 1-D method of Backus & Gilbert (1968). The study region was parametrized into grids with cell size of 0.05° × 0.05° and with proper regularization parameters to provide relatively smooth maps with small data misfits. The same regularizations parameters were used for producing the maps, before and after the impoundment. The resolution of the data set is controlled by the average path length, the density of paths and the azimuthal coverage. In this study, the spatial averaging kernels are estimated for each point (Backus & Gilbert 1968). In addition to group velocity maps, the corresponding resolution information will be provided at each period. The Yanovskaya's methodology (1997) is used to calculate the spatial resolution, which varies from 4 to about 10 km in the study region. This methodology allows calculating the spatial averaging kernels at each point and in different directions. A functional S(x, y) is calculated for different orientations of the coordinate system in order to determine sizes of the averaging area along different directions. The averaging area can be approximated by an ellipse, centred in a point of the studied region, with axes equal to the largest and the smallest values of S(x, y). The size of the averaging area along different directions are defined by S(x, y), which gives an idea of the resolution in each point (x, y). It is approximated by an ellipse with half-axes of the Smax(x, y) and Smin(x, y). The smallest [Smin(x, y)] and largest [Smax(x, y)] axes of the ellipse are calculated, and the resolution at each point is given by the mean size of the averaging area L = [Smin(x, y) + Smax(x, y)]/2. Applying the described tomographic method, the Rayleigh wave group velocity maps for periods ranging from 2 to 8 s are generated. Fig. 5 shows the group velocity distribution at 3, 5 and 8 s with the corresponding resolution maps and ray path coverage in the study region. Only the seismic velocities of the regions where the resolution is ≤10 km (Fig. 5) are taken into consideration for further discussion. Figure 5. View largeDownload slide Rayleigh wave group speed maps derived from ambient noise tomography at periods of 3, 5 and 8 s for before (a-c) and after (g-i) the impoundment. The resolution maps and ray path coverage at periods of 3, 5 and 8 s are presented for before (d-f) and after (j-l) the impoundment, respectively. Figure 5. View largeDownload slide Rayleigh wave group speed maps derived from ambient noise tomography at periods of 3, 5 and 8 s for before (a-c) and after (g-i) the impoundment. The resolution maps and ray path coverage at periods of 3, 5 and 8 s are presented for before (d-f) and after (j-l) the impoundment, respectively. The 3 s group velocity map before impoundment (Fig. 5) shows a distinct low velocity zone in the northern part of the reservoir, south of Lali, which extends as a NW-SE oriented narrow band up to Masjed Soleyman reservoir, located east of city with the same name. This trend of Low Velocity Zone (LVZ) is also visible in 5 and 8 s map, but with a little higher velocity (from 1.5 to 2.0 km s–1), respectively. Another low velocity zone is clear around the dam location at 3 and 5 s period map, but it vanishes at 8 s. A very narrow high velocity zone is in parallel present at the south of the LVZ. The area and velocity of this zone which is visible at all different periods, increases with period from 3 to 8 s. In the eastern part of the area, beneath the Masjed Soleyman reservoir a distinct high velocity zone is observed at all periods. After impoundment, distribution of the LVZ remains almost the same as before, mainly beneath the reservoir. A little increase is observed in velocity and area of the High Velocity Zone (HVZ), which is located around the LVZ, especially on southeastern margin of the study area at 5 and 8 s periods. Beneath, south and north of the Masjed Soleyman reservoir, a clear HVZ is observed at all periods, before and after the Gotvand reservoir impoundment. The impoundment of the Masjed Soleyman dam backs to 2001 about 5 yr before installation of the Gotvand and Masjed Soleyman seismological networks. 4 3-D SHEAR VELOCITY MODEL A 3-D shear wave velocity model was obtained extracting local group velocity dispersion curves for each node across the 0.05° × 0.05° gridpoints of the 2-D tomographic maps in the period range of 2–8 s. We adopt an iterative linearized least-square inversion scheme (Herrmann 2013) to invert the extracted dispersion curve for a 1-D shear wave profile. The linearized inversion technique shows some dependence on the initial model. The initial model is constructed using the available shear wave velocity models of the region obtained from 1-D inversion of the locally recorded earthquakes by the Gotvand seismic network (Tatar et al. 2009). Our initial model is a constant value two-layer velocity model composed of 10 thin layers (0–10 km) each one 1 km thick, overlying a half-space (black dashed line in Fig. 6b). The S-wave velocity in each layer is constant. In order to illustrate the inversion process, an example of Vs inversion for an extracted group velocity dispersion curve at a gridpoint, which its location is shown by a yellow star in Fig. 3, is presented in Fig. 6. Figure 6. View largeDownload slide (a) Predicted group velocity dispersion curve (grey lines) from the ensemble of Vs models shown in (b), and average of the ensemble dispersion curves (blue line). The black line is the observed (extracted) dispersion curve with error bars at individual periods. Result of inversion of the average of Vs models (red line in b) is presented as red line. (b) The resulting ensemble of acceptable Vs models (grey lines) and the average of the ensemble (red line). The initial model is presented as black dashed line, composed of 10 separated layers, each one 1 km thick. Figure 6. View largeDownload slide (a) Predicted group velocity dispersion curve (grey lines) from the ensemble of Vs models shown in (b), and average of the ensemble dispersion curves (blue line). The black line is the observed (extracted) dispersion curve with error bars at individual periods. Result of inversion of the average of Vs models (red line in b) is presented as red line. (b) The resulting ensemble of acceptable Vs models (grey lines) and the average of the ensemble (red line). The initial model is presented as black dashed line, composed of 10 separated layers, each one 1 km thick. To verify the resolution of the inversion result with depth, we plotted the quantity dU(T)/dVs (Fig. 7) which reflects the sensitivity of the fundamental mode Rayleigh-wave group velocity on the S-wave velocity versus depth for one examined cell located on profile AA’ (Fig. 3) for the periods of 2, 3, 5 and 8 s. This figure shows that Rayleigh group velocity at 2, 3 and 5 s are sensitive to variation of shear wave velocity at depths smaller than 2, 3 and 4 km, respectively. At 8 s, the shortest period of our data shows a peak at a depth of about 6 km. So S-wave velocities at depths greater than 6 km do not have a significant influence on the values of group velocities (very low absolute sensitivity), which indicate our results are valid until depth ∼6 km. Deeper anomalies are not reliable due to large errors and low resolution of results. Figure 7. View largeDownload slide Normalized sensitivity kernels for the fundamental mode Rayleigh-wave group velocity, calculated at different periods of 2, 3, 5 and 8 s for a grid located on profile AA’ (Fig. 3). Figure 7. View largeDownload slide Normalized sensitivity kernels for the fundamental mode Rayleigh-wave group velocity, calculated at different periods of 2, 3, 5 and 8 s for a grid located on profile AA’ (Fig. 3). In order to smooth out the Vs profiles and reduce the velocity contrasts over depths, during the inversion each Vs profile was vertically damped to minimize the differences at neighbouring layers. All constructed 1-D Vs profiles for each gridpoint then are assembled to form the final 3-D Vs model. Horizontal sections of our 3-D shear wave velocity model at depths of 3, 6 and 9 km for before and after impoundment are shown in Fig. 8. Figure 8. View largeDownload slide Shear wave velocity maps at depths of 3, 6 and 9 km before (a–c) and after (d–f) impoundment of the Gotvand dam. Figure 8. View largeDownload slide Shear wave velocity maps at depths of 3, 6 and 9 km before (a–c) and after (d–f) impoundment of the Gotvand dam. Vs model at different depths (Fig. 8) reveal almost the same features as the group velocity maps at different periods. Before impoundment of the reservoir, the Vs model at 3 km depth shows a LVZ just in north of the central part of the reservoir, south of Lali, where a distinct low velocity region is observed. A narrow HVZ with NW-SE orientation is visible in southern margin of the reservoir, especially at depth of 6 km. Comparison of the Vs models for before and after impoundment (Fig. 8) suggest a relatively velocity reduction beneath the reservoir at shallow depths. In order to estimate the uncertainty in our shear-wave velocity model we used a bootstrap analysis (Efron & Gong 1983; Efron & Tibshirani 1993); according to this methodology error estimation is based on the inversion of different subsets of data from all produced sets of solutions. In this analysis, we produced 1000 data sets, called bootstrap samples, from the full set of interstation group velocity measurements. Each set produced by randomly selecting 75 per cent of the individual group velocity measurements; then each set is processed by the 2-D tomographic inversion and 1-D shear wave velocity inversion. Therefore, the bootstrap analysis propagated errors through both the 2-D tomographic inversion and 1-D shear wave velocity inversion. Fig. 9 shows the standard deviation of the model ensemble for our shear wave velocity model at shallow depths (<6 km) vary in the range of ∼0.03–0.10 km s–1 for the inversion before (Figs 9a–c) and after (Figs 9d–f) impoundment. Generally larger standard deviation values (up to ∼0.15 km s–1) are found at greater depths (>6 km), which is possibly due to decreasing the number of available dispersion measurements at the longer periods and also relatively poor depth resolution of the Rayleigh waves for longer periods compare to short periods (Fig. 7). Figure 9. View largeDownload slide Uncertainty in the shear wave velocity estimates from 1000 bootstrap resampled data sets at 3, 6 and 9 km depths, before (a–c) and after (d–f) impoundment of the dam. Figure 9. View largeDownload slide Uncertainty in the shear wave velocity estimates from 1000 bootstrap resampled data sets at 3, 6 and 9 km depths, before (a–c) and after (d–f) impoundment of the dam. Two vertical cross-sections across the Gotvand reservoir (Fig. 3), plotted in Fig. 10, reveal the velocity variations in depth before and after impoundment. Corresponding standard deviation in these two vertical cross-sections are also plotted in Fig. 11, showing that uncertainty at shallow parts (2–6 km) almost vary in the range of ∼0.04–0.07 km s–1, which is much better than deeper parts (>6 km)that reaches up to ∼0.10 km s–1. The location of Profile BB’ is the same as geological map cross-section (Fig. 2). There is very good consistency between the geological features and the determined velocity structure. Figure 10. View largeDownload slide Vertical cross-sections of shear wave velocity for before (a and c) and after (b and d) impoundment of the Gotvand dam along two profiles delineated by the black lines in Figs 3 and 8(a). Topography is over plotted above individual cross sections. The location of the lake along the profile is presented as solid grey line above the topography. Figure 10. View largeDownload slide Vertical cross-sections of shear wave velocity for before (a and c) and after (b and d) impoundment of the Gotvand dam along two profiles delineated by the black lines in Figs 3 and 8(a). Topography is over plotted above individual cross sections. The location of the lake along the profile is presented as solid grey line above the topography. Figure 11. View largeDownload slide Uncertainty in shear wave estimated from 1000 bootstrap resampled data sets for before (a and c) and after (b and d) impoundment of the Gotvand dam along two profiles delineated by the black lines in Figs 3 and 8(a). Topography is over plotted above individual cross sections. The location of the lake along the profile is presented as solid grey line above the topography. Figure 11. View largeDownload slide Uncertainty in shear wave estimated from 1000 bootstrap resampled data sets for before (a and c) and after (b and d) impoundment of the Gotvand dam along two profiles delineated by the black lines in Figs 3 and 8(a). Topography is over plotted above individual cross sections. The location of the lake along the profile is presented as solid grey line above the topography. As it is inferred from the cross-section view of the geological map (Fig. 2), Gachsaran is the most dominant formation underneath the reservoir. This sedimentary formation caused by an evaporate environment at the end of the Upper Miocene. As faulting, folding, and flow are the controlling factors for the thickness variations in the Gachsaran formation, so strong differences in seismic images can be observed due to variations in the lithology and thickness caused by an uplifted rock or faulting along this formation (Alavi 1994). The Gachsaran Formation consists of the high velocity anhydrite, salt, and limestone, and low velocity rocks of marl and shale (Abdollahie Fard et al. 2011). For the time period before the impoundment of the dam, a very low velocity zone is observed in the northern part of the reservoir along the lines AA’ (Fig. 10a) and BB’ (Fig. 10c). Comparing with the geological cross section along the BB’ profile, this low velocity zone is related to the Gachsaran formation which is the dominant formation in the area located in northern parts of the reservoir. This gipsy formation that made mainly from low velocity evaporate deposits is spread over a wide area at surface and in depth. The shape, size and S-wave velocity (∼2 km s–1) of observed anomaly north of the lake on BB’ shear wave velocity profile (Fig. 10c) matches perfectly with the specification of Gachsaran formation on geological cross section, which is located north of the lake (Fig. 2). The high velocity zone observed southwest of the low velocity zone of Gachsaran formation, at a depth of ∼2–8 km (Figs 10c and d) shows a very good consistency with the folded and uplifted limestone formations (Fig. 2) of the Asmari-Pabedeh (Vs ∼ 2800 m s–1), Ilam-Sarvak and Kazhdum (Vs ∼ 3300 m s–1). These cross sections (A and C on Fig. 10), reveal that the ANT technique is a power tool with enough good resolution for detecting lithologically different formations even at small scale. It means that very dense small-spacing seismological networks can be used for detecting of hydrocarbon traps applying ANT technique. The effect of the reservoir impoundment is observed in Figs 10(b) and (d). It seems that both low and high velocity anomalies related to the Gachsaran and Asmari, Pabdeh, Ilam-Sarvak group respectively, have been compacted due to loading of the above reservoir water column. A very distinct high velocity anomaly is observed just south of the Gotvand lake (Fig. 10d). Based on field investigation, this high velocity anomaly can be related to the existence of a salt dome which is responsible for very high salinity level of the reservoir water. Due to misestimating of the size and capacity of this salt dome which its northern flank has direct contact with the reservoir water, the reservoir water that became hypersaline (very salty) is almost unusable for agricultural purposes. 5 DISCUSSION In order to better define the impoundment effects, the difference between the shear-wave velocity after and before the impoundment of the Gotvand reservoir was calculated along the AA’ and BB’ cross section (Fig. 12). This Figure shows clearly that along both cross sections we observe a decrease in S-wave velocity (red colour) up to ∼0.3 ± 0.05 km s–1 at shallow depths (2–4 km) just beneath the lake (reservoir). This can be due to the seepage of water at shallow depth, which causes reduction of shear wave velocity. But at deeper parts (4–6 km) beneath the reservoir, an increase of ∼0.3 ± 0.09 km s–1 (blue colour) is observed in S-wave velocity that can be due to compaction effects of the reservoir water. However, due to the large uncertainty in S-wave velocity based on computed resolution kernel curves (Fig. 7) large velocity anomalies at deeper parts are not enough reliable to be taken into consideration for further interpretation. Figure 12. View largeDownload slide Difference in shear wave velocity observed after and before the impoundment (Vs after –Vs before) for drawn vertical cross-sections along AA’ (a) and BB’ (b) profiles (Fig. 8). Solid grey lines above the topography show the location of the lake along the profile. Figure 12. View largeDownload slide Difference in shear wave velocity observed after and before the impoundment (Vs after –Vs before) for drawn vertical cross-sections along AA’ (a) and BB’ (b) profiles (Fig. 8). Solid grey lines above the topography show the location of the lake along the profile. In order to better understand the variations of velocities before and after impoundment due to the effect of water load on the crust beneath the dam, we modelled the incremental strength underneath the dam along the profiles shown in Fig. 3. According to coulombs law, this incremental strength due to reservoir impoundment is given by (Bell & Nur 1978):   $$\left\{ \begin{array}{@{}l@{}} \Delta S = \mu \left( {\Delta {\sigma _n} - \Delta P} \right){\rm{ }} - \Delta {\tau _n},\ \ \ \ \ \ \ \ \\ \Delta P = \Delta {P_u} + \Delta {P_{diff}} \end{array} \right.$$ (1)  $$\Delta {P_u} = B\Delta {\sigma _m}$$ (2)where Δσn and Δτn are the incremental normal and shear stresses on the fault plane due to the water loading, respectively, μ is the coefficient of friction, B is the Skempton's coefficient, Δσm is the mean stress, which is obtained from summation over the diagonal elements of the stress tensor, ΔPu is the undrained pore pressure, occurs instantly in response to undrained loading of water, and after a while pore pressure, ΔPdiff, is diffused from the reservoir into the crust beneath it. The change in strength due to loading can be calculated by knowing the filling history of the dam, location, and geometry of the fault. Then, the Δσn, Δτn, ΔPu and ΔS can be found from the eqs (1) and (2). In order to calculate the strength beneath the AA’ and BB’ lines, we assumed a vertical fault beneath them, while its western block move upward. The normal and the shear stresses can be given by (Jaeger & Cook 1979):   $$\left\{ \begin{array}{@{}l@{}} \Delta {\sigma _n} = l\Delta {\sigma _x}(T) + m\Delta {\sigma _y}(T) + n\Delta {\sigma _z}(T),\\ \Delta {\tau _n} = \pm \sqrt {(\Delta {\sigma _x}{{(T)}^2} + \Delta {\sigma _y}{{(T)}^2} + \Delta {\sigma _z}{{(T)}^2}) - \Delta \sigma _n^2} , \end{array} \right.$$ (3)where l, m, n are the direction cosines of the normal to the fault plane and Δτn is positive when the shear stress favours faulting. We also considered B = 0.7 (Talwani 1997) and μ = 0.75 (Byerlee1978). A negative value for ΔS may imply a weakening or destabilization and a positive value means strengthening or stabilization (Bell & Nur 1978) in the sense of causing the stress state to move toward or away from failure. Load of a reservoir can cause compressive stresses underneath the reservoir, which can lead to stabilizing the underlying crust and produce tensile stresses in horizontal direction on the surrounding area of the reservoir (Liu & Xu 2005). These tensile stresses can create small shallow fractures and facilitate the penetration of water into deeper rocks and change the properties of the crust (Liu et al. 2011). Assuming a homogeneous medium, Boussinesq solutions (Jaeger & Cook 1979) were used to calculate the elastic stress changes due to reservoir loading. The surface of the reservoir is divided into small grid blocks with 100 m × 100 m dimension. Summing all the stress changes in each grid block the total stress change at each point underneath the reservoir is calculated. The final results for incremental strength at total volume of impoundment 2.1 × 109 m3, and lake level of 197 m (Fig. 13) indicate that all parts underneath the lake have a positive value which means stabilization or strengthening. Away from the lake, the results shows negative values which leads to destabilization or weakening. Figure 13. View largeDownload slide Strength changes after impoundment of the Gotvand lake (total volume of impoundment ∼2.1 × 109 m3, and lake level of 197 m) along (a) Profile AA’ and (b) Profile BB’. Solid grey line shows the location of the lake. Figure 13. View largeDownload slide Strength changes after impoundment of the Gotvand lake (total volume of impoundment ∼2.1 × 109 m3, and lake level of 197 m) along (a) Profile AA’ and (b) Profile BB’. Solid grey line shows the location of the lake. Both elastic loading stress and pore pressure diffusion (ΔPdiff) can affect the area underneath the reservoir (Jaeger & Cook 1979). As the largest volume of the water is behind the crest of the Gotvand dam, so it is expected that the effect of loading stress being more effective in this part. This effect can be seen in Fig. 13. Beneath the profile AA’ which crosses the deepest part of the reservoir the amount of incremental strength is higher than what we observed beneath the profile BB’. Compressive stress due to loading the dam can change lithology of the rocks and enhance small fractures in the periphery of the reservoir and increase the pore pressure. Estimated analytical strength due to loading, leads to stabilizing the crust underneath the reservoir, which is in agreement with earlier reports (Roeloffs 1988; Liu & Xu 2005; Santoyo et al. 2010). Shear wave velocity difference for before and after impoundment (Fig. 12) indicates on considerable decrease of Vs (red colour) up to ∼0.4 ± 0.1 km s–1 at very deep parts (depths greater than 8 and 6 km on profile AA’ and BB’, respectively) beneath and outside of the lake where the compaction effects (analytical strength) is less than 0.2 bar. It seems that deep parts of the crust underneath the reservoir (5–8 km depth) is stabilized due to the compressional stress loading, while outside of the reservoir we observe weakening due to extend of the microcracking and fractures as a result of horizontal tensile stress (Liu et al. 2011). It should be taken into consideration that the presence of incompetent lithology in some parts can facilitate the permeability and migration of water in to deeper parts (Talwani et al. 2007). Finally, the loading effect of the Gotvand reservoir on underneath strength changes is calculated along the profile AA’ by incorporating the effect of pore pressure (ΔPdiff) in lowering the strength of rocks (Fig. 14a). Figure 14. View largeDownload slide Vertical cross-section 3-D view of (a) distribution of pore pressure and (b) strength changes along the profile AA’ after the impoundment of the lake (total volume of impoundment ∼2.1 × 109 m3, and lake level of 197 m) by adding the pore pressure effect (coordinate system in UTM). Figure 14. View largeDownload slide Vertical cross-section 3-D view of (a) distribution of pore pressure and (b) strength changes along the profile AA’ after the impoundment of the lake (total volume of impoundment ∼2.1 × 109 m3, and lake level of 197 m) by adding the pore pressure effect (coordinate system in UTM). Diffused pore pressure (ΔPdiff) can be calculated by solving the following equation at a time t and depth z:   \begin{eqnarray} \Delta {P_{{\rm{diff}}}}(z,t) &=& (1 - \alpha ){p_0}\,{\rm{erfc}}\,\left[ {\frac{z}{{{{(4ct)}^{1/2}}}}} \right]\, + \,\alpha H(t){p_0}\\ \alpha &=& \,B(1 + {v_u})/[3(1 - {v_u})], \end{eqnarray} (4)where the erfc is the complementary error function, H(t) is the Heaviside unit step function, B is the Skempton's coefficient, vu is the undrained Poisson's ratio, z is the depth beneath the reservoir, p0 is the pressure at the bottom of the reservoir, c is the hydraulic diffusivity and t is time. The hydraulic diffusivity constant in eq. (4) was estimated to be c = 0.1 m2 s–1, using, c = r2/4Δt, where Δt is the time lag between the start of the impoundment and the time of estimation and r is the hypocentral distance from the dam (Talwani et al. 2007). Fig. 14(b) shows that due to penetration of water or pore pressure diffusion, which affects shallow depths, decrease in rock strength value is expected. However, in deeper parts, due only to loading effect, an increase in rock strength is observed. These observations are in good agreement with shear wave velocity decrease at shallow depths and its increase at deeper parts beneath the reservoir as it can be seen in Fig. 12 after impoundment of the dam. 6 CONCLUSION We obtained in this study a detailed image of the seismic velocity structures in the top ∼6 km of the region beneath the second largest dam reservoir in south of Iran using Rayleigh wave group velocity measurements from cross-correlation of ambient noise data in the period range of 2–8 s. Group velocity maps at 3, 5 and 8 s periods are obtained by surface wave tomography. The observed velocity structures primarily correlate well with geological features. 3-D shear wave velocity model was obtained by constructing a local group velocity dispersion curves for each node of the 0.05° × 0.05° grid of the 2-D tomographic maps in the period range of 2–8 s. The obtained images for shear wave velocities reveal clear velocity contrasts across the study region. A very low velocity zone is observed beneath the northern part of the reservoir related to the Gachsaran formation which is the dominant formation in the area located in northern parts of the reservoir. The high velocity zone observed southwest of the low velocity zone of Gachsaran formation, at a depth of ∼4–8 km consistent with the folded and uplifted limestone formations of the Asmari-Pabedeh, Ilam-Sarvak and Kazhdum. Calculated difference between the shear wave velocity after and before the impoundment of the Gotvand reservoir shows a clear decrease in S-wave velocity up to ∼0.3 ± 0.05 km s–1 at shallow depths (2–4 km) just beneath the lake (reservoir) which can be related to the seepage of water in shallow depths. At deeper parts (4–6 km) beneath the reservoir, an increase of ∼0.3 ± 0.1 km s–1 is observed in S-wave velocity that can be associated to compaction due to loading effects of the reservoir water. Calculation of loading effects of the Gotvand reservoir on underneath strength changes, considering the effect of pore pressure, show good agreement with the observed shear-wave velocity changes across the reservoir. The results of this research show that building a dam can make severe changes in the crust properties underneath it. Considering the lithology, volume of the water and location of the dam, these changes could be different at distinct distances and depths. Acknowledgements This research was supported and has been done as part of the contract between International Institute of earthquake Engineering and Seismology (IIEES), the Iranian Water and Power Resources Development Co. (IWPCO)-Upper Gotvand Project and the Khuzestan Water and Power Authority (KWPA) for seismic monitoring of the Gotvand and Masjed Soleyman reservoirs. The data used in this work were recorded by the Gotvand-Masjed Soleyamn permanent Seismic Network operated by IIEES. We are grateful to the colleagues of IIEES, IWPCO-Upper Gotvand Projects, and KWPA Masjed Solymand dam for helps and supports during installation, gathering, and preparing of the data. Special thanks to Dr Dariush Mahjoub from IWPCO and Dr Yaghob Arab from KWPA for their continuous supports. We thank Regione Friuli-Venezia Giulia (Cooperazione Internazionale e allo Sviluppo) for support through the project “Sostegno all Éducazione e alla Ricerca nellámbito della prevenzione del rischio sismico in Iran”. We gratefully thank and acknowledge two anonymous reviewers for constructive reviews. Some figures were made using the public domain Generic Mapping Tools (Wessel & Smith 1998). REFERENCES Abdollahie Fard I., Sepehr S., Sherkati S. 2011. Neogene salt in SW Iran and its interaction with Zagros folding, Geol. Mag. , 148, 854– 867. Google Scholar CrossRef Search ADS   Alavi M., 1994. Tectonics of the zagros orogenic belt of iran; new data and interpretations, Tectonophysics , 229, 211– 238. Google Scholar CrossRef Search ADS   Backus G., Gilbert F., 1968. The resolving power of gross Earth data, Geophys. J. R. astr. Soc. , 16, 169– 205. Google Scholar CrossRef Search ADS   Bell M.L., Nur A., 1978. Strength changes due to reservoir-induced pore pressure and stress and application to Lake Oroville, J. geophys. Res. , 83, 4469– 4483. Google Scholar CrossRef Search ADS   Bensen G.D., Ritzwoller M.H., Barmin M.P., 2007. Processing seismic ambient noise data to obtain reliable broad-band surface wave dispersion measurements, Geophys. J. Int. , 169, 1239– 1260. Google Scholar CrossRef Search ADS   Berberian M., 1995. Master blind thrust faults hidden under the Zagros folds: active basement tectonics and surface morphotectonics, Tectonophysics , 241, 193– 224. Google Scholar CrossRef Search ADS   Brenguier F., Shapiro N.M., Campillo M., Nercessian A., Ferrazzini V., 2007. 3D surface wave tomography of the Piton de la Fournaise volcano using seismic noise correlations, Geophys. Res. Lett. , 34, L02305, 1– 5. Google Scholar CrossRef Search ADS   Byerlee J., 1978. Friction of rocks, Pure appl. Geophys. , 122, 947– 965. Ditmar P.G., Yanovskaya T.B., 1987. Generalization of the Backus-Gilbert method for estimation of the horizontal variations of surface wave velocity, Izv. Acad. Sci. USSR Phys. Solid Earth. (English Transl.) , 23, 470– 477. Droznin D.V., Shapiro N.M., Droznina S.Y.A., Senyukov S.L., Chebrov V.N., Gordeev E.I., 2015. Detecting and locating volcanic tremors on the Klyuchevskoy group of volcanoes (Kamchatka) based on correlations of continuous seismic records, Geophys. J. Int. , 203, 1001– 1010. Google Scholar CrossRef Search ADS   Dunnington H.V., 1968. Salt-tectonic features of northern Iraq. In: MATTOX, R. B. (ed.) A symposium based on paper from the International Conference on Saline Deposits, Geol. Soc. Am., Special Papers , 88, 183– 227. Google Scholar CrossRef Search ADS   Efron B., Gong G., 1983. A leisurely look at the Bootstrap, the Jackknife, and cross validation, Am. Stat. , 37, 36– 48. Efron B., Tibshirani R.J., 1993. An Introduction to the Bootstrap , Chapman and Hall/CRC. Google Scholar CrossRef Search ADS   Guidarelli M., Aoudia A., 2016. Ambient noise tomography of the Cameroon Volcanic Line and Northern Congo craton: new constraints on the structure of the lithosphere, Geophys. J. Int. , 204, 1756– 1765. Google Scholar CrossRef Search ADS   Guo Z., Gao X., Yao H.J., Li J., Wang W.M., 2009. Midcrustal low-velocity layer beneath the central Himalaya and southern Tibet revealed by ambient noise array tomography, Geochem. Geophys. Geosyst. , 10( 5), Q05007, 1–12. Google Scholar CrossRef Search ADS   Herrmann R.B., 2013. Computer programs in seismology: an evolving tool for instruction and research, Seism. Res. Lett. , 84, 1081– 1088. Google Scholar CrossRef Search ADS   Jaeger J.C., Cook N.G.W., 1979. Fundamentals of Rock Mechanics , 3rd edn. Chapman and Hall. James G.A., Wynd J.G., 1965. Stratigraphic nomenclature of Iranian oil consortium agreement area, AAPG Bull ., 49, 2182– 2245. Hessami Azar K., Jamali F., Tabasi H., 2003. Major Active Faults of Iran , IIEES Publication. Koulakov I.et al.  , 2014. Asymmetric caldera-related structures in the area of the Avacha group of volcanoes in Kamchatka as revealed by ambient noise tomography and deep seismic sounding, J. Volc. Geotherm. Res. , 285, 36– 46. Google Scholar CrossRef Search ADS   Lin F.C., Ritzwollerl M.H., Townend J., Bannister S., Savage M.K., 2007. Ambient noise Rayleigh wave tomography of new Zealand, Geophys. J. Int. , 170, 649– 666. Google Scholar CrossRef Search ADS   Liu S.M., Xu L.H., 2005. Finite-element simulation of hydraulic-pressure stress field in Danjiangkou Reservoir area, J. Hydraul. Eng. , 36, 863– 869. Liu S.M., Xu L.H., Talwani P., 2011. Reservoir-induced seismicity in the Danjiangkou Reservoir: a quantitative analysis, Geophys. J. Int. , 185, 514– 528. Google Scholar CrossRef Search ADS   Lobkis O., Weaver R.L., 2001. On the emergence of the Green's function in the correlations of a diffuse field, J. acoust. Soc. Am. , 110( 6), 3011– 3017. Google Scholar CrossRef Search ADS   Mcquarrie N., 2004. Crustal scale geometry of the zagros fold–thrust belt, Iran, J. Struct. Geol. , 26, 519– 535. Google Scholar CrossRef Search ADS   Mouthereau F., Lacombe O., Meyer B., 2006. The Zagros folded belt (Fars, Iran): constraints from topography and critical wedge modelling, Geophys. J Int. , 165, 336– 356. Google Scholar CrossRef Search ADS   Nissen E., Tatar M., Jackson J.A., Allen M.B., 2011. New views on earthquake faulting in the Zagros fold-and-thrust belt of Iran, Geophys. J. Int. , 186( 3), 928– 944. Google Scholar CrossRef Search ADS   O’brien C.A.E., 1950. Tectonic problems of the oil field belt of southwest Iran, Proceedings of 18th International Geological Congress, part 6 , Great Britain, pp. 45– 58. Ritzwoller M.H., Levshin A.L., 1998. Eurasian surface wave tomography: group velocity, J. geophys. Res. , 103, 4839– 4878. Google Scholar CrossRef Search ADS   Roeloffs E.A., 1988. Fault stability changes induced beneath a reservoir with cyclic variations in water level, J. geophys. Res. , 93, 2107– 2124. Google Scholar CrossRef Search ADS   Santoyo M., García-Jerez A., Luzón F., 2010. A subsurface stress analysis and its possible relation with seismicity near the Itoiz Reservoir, Navarra, Northern Spain, Tectonophysics , 482, 205– 215. Google Scholar CrossRef Search ADS   Sattarzadeh Y. 1997. Active tectonics in the Zagros Mountains, Iran , PhD thesis, Imperial College, London. Setudehnia A., 1972. Iran du Sud-Ouest: Lexique Strat. Internat, Centre Nat. Rech. Scientifique , Paris, III, Asie, Fasc. 9b, p. 289– 376. Shapiro N.M., Campillo M., Stehly L., Ritzwoller M.H., 2005. High-resolution surface wave tomography from ambient seismic noise, Science , 307, 1615– 1618. Google Scholar CrossRef Search ADS PubMed  Snieder R., 2004. Extracting the Green's function from the correlation of coda waves: a derivation based on stationary phase, Phys. Rev. E. , 69, 046610. Google Scholar CrossRef Search ADS   Talwani P., 1997. On the nature of reservoir-induced seismicity, Pure appl. Geophys. , 150, 473– 492. Google Scholar CrossRef Search ADS   Talwani P., Acree S., 1984/1985. Pore pressure diffusion and the mechanism of reservoir-induced seismicity, Pure appl. Geophys. , 122, 947– 965. Google Scholar CrossRef Search ADS   Talwani P., Chen L., Gahalaut K., 2007. Seismogenic permeability, ks, J. geophys. Res. , 112, B07309, 1–18. Google Scholar CrossRef Search ADS   Tatar M., Yamini Fard F., Dezvareh M. 2009. Seismicity and seismotectonic study of the Gotvand-e Olya dam region using local earthquakes, annual report, IIEES, pp. 120. Wessel P., Smith W.H.F., 1998. New, improved version of the Generic Mapping Tools released, EOS, Trans. Am. geophys. Un. , 79, 579. Google Scholar CrossRef Search ADS   Yanovskaya T. B., 1997. Resolution estimation in the problems of seismic ray tomography, Izv. Phys. Solid Earth , 33( 9), 762– 765. Yanovskaya T.B., Ditmar P.G., 1990. Smoothness criteria in surface wave tomography, Geophys. J. Int. , 102, 63– 72. Google Scholar CrossRef Search ADS   Yao H., Van der Hilst R.D., De Hoop M.V., 2006. Surface-wave tomography in SE Tibet from ambient seismic noise and two-station analysis—I. Phase velocity maps, Geophys. J. Int. , 166, 732– 744. Google Scholar CrossRef Search ADS   © The Authors 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society.

Journal

Geophysical Journal InternationalOxford University Press

Published: Jan 1, 2018

DeepDyve is your personal research library

It’s your single place to instantly
that matters to you.

over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month Explore the DeepDyve Library Search Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly Organize Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place. Access Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals. Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. DeepDyve Freelancer DeepDyve Pro Price FREE$49/month
\$360/year

Save searches from
PubMed

Create lists to

Export lists, citations