مجله ژئوفیزیک ایران
سال دهم شماره 1 (پیاپی 29، بهار 1395)

Rain is one of the most important climatic factors affecting human activities which has also an important role in the field of water resources management. This weather phenomenon is a complex atmospheric process, which is highly dependent on space and time and thus not easy to predict. The trends of change in rainfall with time is a non-stationary stochastic process with high uncertainty and it is subject to various random factors. There have been many attempts to find the most appropriate method for rainfall prediction using for example meteorological or satellite data with a numerical weather prediction model, or even applying several techniques such as the artificial neural network or fuzzy logic as a forecasting approach. Also some methods, such as the time sequence method, probability statistics method cannot fully reflect the characteristics of the rainfall phenomenon, and the prediction results cannot be satisfactory. In order to improve the accuracy of rainfall forecasts, it is necessary to use a new rainfall prediction model such as intelligent methods and meta-heuristic algorithms. In this study, the “imperialist competitive algorithm” (ICA for brevity) and the ICA combined with the fuzzy logic algorithm were used to evaluate and compare their performance and ability in forecasting the amount of daily rainfall in semi-arid climate of Kerman in the southeast of Iran. So, 30 years of daily data in Kerman’s synoptic station (1981–2010) and 10 years of daily data in Zarand and Rafsanjan’s synoptic stations (2001–2010) were used in the rainy season (7 months of the year). Therefore, based on the previous studies, five parameters including precipitation, wet temperature, dew point, relative humidity and cloudiness were used to forecast rainfall in futures days. Having surveyed the data, first the applied computer codes were written in Matlab 14. In the ICA with fuzzy logic, the ICA was used for determining the membership functions’ ranges and values of the weights instead of the trial and error usually used in application of the fuzzy logic. Three higher accurate outputs were identified for each station separately. Among these outputs, for each station, the best output was chosen and used for the final phase of optimization. Four more effective variables in Kerman’s station (precipitation, wet temperature, dew point, and cloudiness), two more effective variables in Rafsanjan’s station (precipitation and cloudiness) and three more effective variables in Zarand’s station (precipitation, wet temperature, and relative humidity) were identified after optimizing with five input variables. Results showed that the rainfall height’s prediction was accompanied with a significant error based on the mentioned methods, so that the coefficients of determination (R2values) were obtained 0.54, 0.44 and 0.40 in, respectively, Kerman, Rafsanjan and Zarand’s synoptic stations. On the other hand, the forecast of the occurrence and non-occurrence of the rainfall with the ICA indicated reasonable results and in the best results 61.4%, 51.9% and 51.2% of days were predicted correctly in, respectively, Kerman, Rafsanjan and Zarand’s synoptic stations. The accuracy of calculations was improved with the ICA combined with the fuzzy logic. Accordingly, 89.63%, 82.31% and 74.12% of days were predicted correctly in, respectively, Kerman, Rafsanjan and Zarand’s synoptic stations. The results of evaluating the performance showed that the ICA can produce a relatively appropriate simulation of the occurrence and non-occurrence of rainfall in future days, but falls short of ability to simulate the rainfall height properly. On the other hand, the combined ICA and fuzzy logic algorithm provides a better simulation of problems involving high uncertainty.

Ground Penetrating Radar (GPR) method has been extensively employed to map shallow subsurface targets. This method has been widely used to image faults, deformations and discontinuity in network inside rock in engineering and geological studies. The first aim of this study is to show applicability of GPR method in geotechnical study of conductive area (similar to the area around a dam) and the second goal of the study is to verify the capability of a new generation of GPR system (Loza) in imaging of deep targets. In addition, the main features of this system are described in this paper. GPR signal propagation is strongly controlled by water content. It has limited performance in fine grain soils such as clays, marl and silts, or in saline groundwater, all of which strongly attenuate signals. The main restriction of the method is the limited penetration especially in the areas with conductive materials. Unlike the common GPR systems, the Loza system can penetrate deep into the ground even in the conductive areas (up to 250 m). The GPR Loza is a portable, enhanced-power ground penetrating mono-pulse radar developed by VNIISMI Ltd. A distinctive feature of this instrument as compared to other commercial GPR systems is an increase in the transmitter peak power by a factor of 10000 to work in environments with high conductivity. The 10 KW high-power transmitter with 25 MHz unshielded antenna and 6 m length was applied for subsurface study to the depth of 10 m. The average velocity of subsurface was chosen 0.11 m(ns)−1according to the subsurface materials. With the Loza GPR system, the highpower transmission of radar waves in asynchronous mode are recorded with resistive loaded dipole receivers. Geophysical study by GPR method was carried out to find the geophysical properties around Khansar dam. Khansar dam is an earth dam with clay core and 5 million cubic meters reservoir capacity. The dam is of 770 m length, 38 m height, and 10 m crest width. It is located south of Khansar city in the east of Zagrous Chain Mountains. The main geological layers in the study area are limestone, schist and young alluvium. The objective of the study was to investigate contacts of the clay core of the dam with the bedrock and alluvium, groundwater level and the channel and cavities in the bedrock, alluvium and dam. A total of 9916 m parallel and perpendicular profiles was designed for achieving these purposes. The profiles over the dam consist of two groups of profiles. First, profiles that were carried out along the dam, over the rip rap and berms and second, the profiles that were carried out over the dam along the river. Three profiles have been selected for interpretation. Despite the conductive area of the dam, suitable data was received from the depth of 100 m. For each profile, the geological model is designed based on the interpretation and analysis of GPR data. Furthermore, the Krot software was applied for processing and data interpretation. Several anomalies have been detected based on the GPR processed data and geological information. Moreover, Geological layers and the bedrock (which is crashed along some profile) have been detected in the radargrams. In addition, a buried channel is distinguished in the profiles which is located over a crashed zone. The buried channel, weak and heterogeneous zones are interpreted in the GPR radargrams and plotted in the map of the dam. Separation of the geological structure to the depth of 100 m verifies the applicability of this system. Subsequent drilling results in the dam area approve the results of the GPR data.

Increasing the accuracy of numerical methods used for simulation of fluid dynamics problems, particularly the geophysical fluid dynamics problems (atmospheric and oceanic), has been the subject of many research works. Recently, due to the increasing computing power of computers, the advantage of high-resolution numerical methods for numerical simulation of the governing equations of fluid flow is further emphasized. The idea of compact finite difference methods goes back to some works conducted in 1920s and 1940s. However, the pioneering works of Kreiss and Oliger (1972), Hirsh (1975) and Lele (1992) made these methdos popular and showd that compact finite difference methdos can be used as a powerful tool for numerical simulation of fluid dynamics problems appearing in different branches of science. These methods have also been used in numerical simulation of geophysical fluid dynamics problems. Due to the promising performance of compact finite difference methods, application of these schemes to numerical simulations of atmospheric and oceanic flows has increased. The compact finite difference schemes have shown that are able to provide a simple way to reach one of the main objectives in the development of numerical algorithms, i.e., having in our disposal a low-cost and highly-accurate computational method. The compact methods have been used extensively for numerical solution of various fluid dynamics problems. These methods have also been applied to numerical solution of some prototype geophysical fluid dynamics problems (e.g., shallow water equations). Most of the compact finite difference methods are symmetric (usually with a 3- or 5-point stencil) and finding each derivative requires a matrix inversion. Recently, a new class of highly-accurate explicit MacCormack type methods has been introduced for computational fluid dynamic. The compact MacCormack type methods were developed by Hixon and Turkel (2000) to split the derivative operator of the central compact method into two one-sided forward and backward operators. This study is devoted to application of the fourth-order compact MacCormack method for numerical solution of the conservative form of two-dimensional non-hydrostatic and fully compressible Navier Stokes equations governing an inviscid and adiabatic atmosphere. Moreover, the second-order MacCormack method is used to compare the performance of the computations. This enables us to measure some aspects of the computational results (such as efficiency and accuracy). Various aspects of the computations such as discretization of the equations for the interior and boundary points, the details of implementation of boundary conditions for different boundary types (e.g., rigid and open boundaries), time step, grid resolution and dissipation are presented. Since, unlike the second-order MacCormack method, the forward operator in the fourth-order compact MacCormack method for approximation of the first derivative at an arbitrary grid point (e.g, j) is not equal to the backward operator at its adjacent point (i.e., j 1), the application of the fourthorder compact MacCormack method for spatial discretization of the source term in vertical momentum equation in non-hydrostatic models needs special treatment. In this work we have used the conventional second-order MacCormack method (MC2), the standard fourth-order compact MacCormack method (MC4) developed by Hixon and Turkel (2000) and a fourth-order compact MacCormack method with a four-stage Runge-Kutta for time advancing (MCRK4) in our numerical simulations. To evaluate the performance of these methods, two test cases including evolution of a warm bubble, and evolution of a cold bubble in a netural atmosphere were simulated. To simulate cold bubble, the test case presented by Straka et al. (1993) and for simulation of warm bubble, the test case of Mendea-Nunez and carrol (1993) are used. Qualitative and quantitative assessment of the results for different test cases showed the superiority of the MCRK4 and MC4 methods over the MC2 method.

Gradient methods have been extensively used in the interpretation of near-surface potential field anomalies because they increase the resolution of the edges of geologic sources and decrease the effect of the background or ‘regional’ field. Two gradient-based methods, namely the Euler and analytic signal, have proven particularly effective in the semi-automatic interpretation of magnetic anomaly data (i.e., determination of source outline and the depths of sources). Aeromagnetic data often contain anomalies with a large range in amplitude. Processed aeromagnetic images, such as horizontal and vertical derivatives, similarly contain features with large and small amplitudes. The smaller amplitude anomalies might be of considerable geological interest, but they can be hard to delineate among those of larger amplitude. Automatic Gain Control (AGC) filters divide the data at each point by an average value computed in a window centered on that point, thus produce a balanced image where all anomalies have similar amplitude. However, their output depends strongly on the window size, (Rajagopalan and Milligan, 1995), which makes them unsuitable for applications involving aeromagnetic data sets that in most cases contain anomalies of different size. In the last few years, there has been a surge of interest in the development of filters based on the derivatives of aeromagnetic data that produce a balanced output but avoid the window-size problem of the standard AGC filter. The first filter to be introduced was the tilt angle, (Miller and Singh, 1994), which is a balanced vertical derivative. Verduzco et al. (2004) suggest using the total horizontal derivative of the tilt angle as an interpretive tool. Wijns et al. (2005) developed the theta map, a balanced total horizontal derivative, Cooper and Cowan (2006) introduced a horizontal tilt angle, which is a balanced horizontal derivative. Another approach to edge detection based on derivatives was the use of the balanced windowed standard deviation. The analytic signal, which is a complex function and makes use of the Hilbert transform, has been shown to be effective in interpretation of the subsurface magnetic contacts (e.g., Nabighian, 1972, 1974). A geologic contact or fault with significant susceptibility contrast is detected by mapping the maxima of the simple analytic signal, which is composed of the two horizontal and one vertical gradient. The analytic signal response from the larger and shallow sources is clearly visible, but it is very faint for deeper bodies. To enhance this image and bring out the detail in smaller amplitude anomalies, in this paper a balanced analytic signal was computed as followed (Cooper, 2009):
where xH and yH are the Hilbert transforms of the analytic signal in the x- and y- directions, respectively. The width of anomalies in the balanced analytic signal amplitude image has increased slightly compared with those in the original analytic signal amplitude image, which might make accurate contact mapping more difficult. The curvature of geophysical data can be a useful attribute, and in the past it has been applied to potential-field data (Phillips et al., 2007) and reflection seismic data (Blumentritt et al., 2006). In this case, the profile curvature was used, defined as (Mitasova and Jarosalav, 1993):
where,p=(∂f/∂x)2(∂f/∂y)2 andq=p. In this paper we demonstrated the application of these filters to synthetic magnetic data and as well as real magnetic data from Soork Iron mine and the aeromagnetic data from Yilgran Craton in Australia. Application of this filter to the magnetic data from Soork Iron ore lead to detection of a new hidden body which is proved by exploratory drilling.

OjatAbad iron ore located in north east of Semnan city. Old indications of mining are evident in the area. Belich and Bragin (1993) introduced Semnan iron ores as hydrothermal deposits. Recent studies show that the iron ores are related to Oligo-Miocene magmatism. Most iron ores have magnetite which has high magnetic susceptibility; therefore, magnetic method is a conventional method for geophysical exploration of iron.
In order to identify and detect this deposit, a magnetic survey was carried out in OjatAbad area. The magnetic data corrected for diurnal change of magnetic field and then total magnetic field of the Earth has been reduced. To do this, a reduce to pole (RTP) filter was implemented on the grid for locating the anomalies and their sources. This method entails removing the dependence of magnetic data to the magnetic inclination, i.e., converting the data which were recorded in the inclined Earth’s magnetic field to what they would have been if the magnetic field had been vertical. This method simplified the interpretation because for sub-vertical prisms or sub-vertical contacts (including faults), it transforms their asymmetric responses to simpler symmetric and anti-symmetric forms. The symmetric “highs” are directly centered on the body, while the maximum gradient of the anti-symmetric dipolar anomalies coincides exactly with the body edges.
For depth estimation of anomalies, the upward continuation filter was implemented. This is a mathematical technique that projects the data taken at an elevation to a higher elevation. The effect is that the short-wavelength features are smoothed out because one is moving away from the anomaly. The upward continuation is a way of enhancing large scale (usually deep) features in the survey area. It attenuates the anomalies depending on their wavelengths; the shorter the wavelength, the greater the attenuation. Also upward continuation tends to accentuate the anomalies caused by deep sources at the expense of the anomalies caused by shallow sources. For 3D imaging of magnetic data, we chose the inversion method of Li and Oldenburg (1998) that minimizes a function composed by (1) the data-misfit function defined in the data space as the L2 norm of the difference between the observed and predicted data, and (2) the stabilizing function defined in the parameter model space as the L2 norm of the first-order derivative of the weighted density distribution in both vertical and horizontal directions. They introduced a depth-weighting function to counteract the spatial decay of the kernel function with depth by giving more weight to rectangular prisms as depth increases.
On the RTP magnetic map of study area we can see 8 magnetic anomalies (A, B, C, D, E, F, G and H) which are located at the northeast to the southwest trend. The H one has the lowest amplitude. The results of upward continuation filter show that the anomalies of F and G are shallower than the other anomalies. Also, H and B are only deeper than 100m. The inversion results recovered all of the 8 anomalous bodies and confirm the above results. They showed that the anomalous body H has lower magnetic susceptibility and is deeper than the others and it seems likely that it is an intrusive body. So, iron mineralization is probably happened in the other bodies. The anomalous body B is the deepest mineralized body which elongates to 100 meter.

The Madden–Julian oscillation (MJO) is a large-scale coupling between atmospheric circulation and tropical deep convection explaining a large part of the intra-seasonal (30–90 day) variability in the tropical and extra-tropical atmospheres. Rather than being a standing pattern like the El Niño–Southern Oscillation (ENSO), the MJO is an oceanic–atmospheric moving pattern that propagates eastward at approximately 4 to 8 ms−1, through the atmosphere above the warm parts of the Indian and Pacific oceans. Due to its quasi-cyclic pattern, the MJO is known as the 30–60 day oscillation, 30–60 day wave, or intra-seasonal oscillation. The oscillation is characterized by an eastward progression of large regions of both enhanced and suppressed tropical precipitation, observed mainly over equatorial parts of the Indian and Pacific Oceans. The anomalous convective precipitation is firstly triggered over the western Indian Ocean, and persists as it propagates over the warm ocean waters of the western and central tropical Pacific. Wheeler and Hendon (2004) categorized the whole cycle of the MJO into 8 phases namely phase 1 to 8. Phase 1 of the MJO that is the theme of this study, signifies the spells for which the convective activity is mainly centered over the western parts of the Indian Ocean equator. It has been previously reported that the occurrence of the MJO phase 1 improves the probability of precipitation event in southwestern Iran (Nazemosadat and Ghaedamini, 2010). In spite of these findings, the reasons of frequent clear sky and shiny days during this phase were not yet resolved. The aim of this study was to compare the oceanic and atmospheric features of the phase 1 for the spells that the incidence of this phase is concurrence with or without precipitation in the southwest of Iran. This analysis is beneficial for improving the MJO-based precipitation and climate forecast over this area. To accomplish this task, the MJO amplitudes during phase 1 were extracted for all days during the 1975–2012 period. Among the selected dates, those days with amplitudes greater than 1 were assigned as the strong MJO events. These MJO events were then divided into two different parts comprising the events with and without pervasive precipitation in southwestern Iran. According to the adopted definition, the pervasive precipitation occurred when at least five out of nine considered stations over the study area received more than 1.0 mm precipitation. For the pervasive precipitating dates, the MJO-precipitation composites were constructed for each individual station. Similar composites were also constructed for precipitation, vector wind, outgoing long-wave radiation (OLR) and vapor flux over the Middle Eastern region and tropical parts of the Indian and Pacific Oceans using the ESRL–NOAA composite facilities. Similar compositing procedure was also performed for the no-precipitating dates to investigate atmospheric condition during such spells.From about 8% to 15% of the November–April precipitation in the southwest of Iran was found to be associated with the periods that the MJO was centered in its phase1. It was concluded that, for improving the MJO-based precipitation forecast in southwest of Iran, not only the phase number and amplitude size of the MJO index, but the position and intensity of convective activities as well as the atmospheric circulations over the Indian and the Pacific Oceans should also be analyzed. In phase 1, precipitation event in southwest of Iran is usually associated with the state of convective activities and their relevant airflows over some areas extended from equatorial parts of east Africa up to the equatorial areas of the Indian Ocean around 90o E. In general, precipitation occurs over the study area if the core center of these activities locate around 60o E. For such situation precipitation anomaly over this area is greater than 1.5 mm/day and the OLR anomalies are less than −20Wm−2or lower. Coincidence with the Indian Ocean equatorial area, precipitation event over the study area in Iran is harmonized with the enhancement of convective activities over tropical parts of the Pacific Ocean for the areas between 160o W to 140o W and 10o S to 20o S. Compared to non-precipitating periods of phase 1, easterly or southerly wind enhances by about four times over the eastern parts of the Indian Ocean, tropical parts of North Africa and the Arabian Peninsula during the precipitating spells. The Persian Gulf was found to play an influential role for re-moisturizing the southerly airflows crossing this water body during the precipitating dates.

In this study¡ a modified “probabilistic seismic hazard assessment” (PSHA) method is used to estimate the level of the potential seismic ground motion in Firoozkouh. A problem that may be encountered in probabilistic studies of seismic hazard for a specific site¡ for engineering purposes¡ is the selection of design earthquakes corresponding to a given hazard value. In order to derive a seismic scenario consistent with the results of PSHA for a site and determine the relative contribution of events to the overall seismic hazard¡ the concept of disaggregation was introduced. The disaggregation of seismic hazard is an effective way to identify the scenario events that contribute to a selected seismic-hazard level. In other words¡ the disaggregation process separates the contributions to the mean annual rate of exceedance (MRE) of a specific ground motion value at a site due to scenarios of given magnitude M¡ distance R¡ and the ground motion error term¡ ε. Disaggregation results could change with the spectral ordinate and return period¡ thus more than one single event may dominate the hazard especially if multiple sources affect the hazard at the site. These results can provide useful information for better defining the design scenario and selecting corresponding time histories for seismic design. In most cases¡ as the probability decreases¡ the hazard sources closer to the site dominate. Larger¡ more distant earthquakes contribute more significantly to hazard for longer periods than shorter periods. In this study¡ the seismic hazard disaggregation process is performed to identify dominant scenarios in “peak ground acceleration” (PGA) and 5% damped 0.2 and 2.0 s spectral accelerations corresponding to mean return periods (MRPs) of 50 yr¡ and 475 yr (hazard levels of 63% and 10% probability of exceedance in 50 yr¡ respectively) in Firoozkouh city. In this regard¡ potential seismic sources and their seismicity parameters have been estimated based on the concept of spatial distribution function in 34º–37ºN latitudes and 52º–55º E longitudes in grid intervals of 0.1º. For each point using proper attenuation relationships¡ PGA¡ 0.2 and 2.0 s spectral acceleration values with 63% and 10% probabilities of exceedance in 50 yr have been calculated using the EZ-FRISK (version 7.43) code. The hazard can be simultaneously disaggregated in different types of bin. The result of seismic hazard disaggregation are presented in terms of 1-D M¡ R and ε bins and 2-D M-R bins. Bins of width 0.4 in magnitude¡ 10 km in distance¡ and 0.2 in ε are selected. The disaggregation results in terms of probability density function (PDF) are reported¡ which is obtained by dividing the probability mass function (PMF) contribution of each bin by the bin’s size¡ thus the PDF representation is independent of the bin’s amplitude. The results identify the distribution of the earthquake scenarios that contribute to exceedance of PGA and 5% damped 0.2 and 2. s spectral accelerations for 50 yr and 475 yr MRPs¡ in terms of magnitude and distance (M-R). Dominant scenarios are identified for interest hazard levels in Firoozkouh city.

Northwestern Iran is one of the seismically active regions with a high seismic risk in the world as confirmed by the historical background and instrumental earthquakes. The Ahar–Varzeghan double earthquakes of August 11th, 2012 with magnitudes higher than 6 are the examples of high seismic activities in this area. Our knowledge of stress state in a region is useful for a better understanding of different rupture mechanisms. The accumulation of stress is the main cause of earthquake, and its analysis is a crucial task considering the dense population of the region. The focal mechanism of an earthquake is one of the important source parameters which is practical for studying and analyzing the stress field. This parameter is affected by fault geometry and principal stress directions. One of the big advantages of the focal mechanism solutions is the ability to study the stress regime depp within the lithosphere. In this study, we determine the focal mechanisms of 15 earthquakes using the moment tensor inversion method of ISOLA program (Sokos and Zahradnik, 2008) in Ahar–Varzeghan region. This method was first proposed in order to calculate the source parameters at teleseismic distances (Kikuchi and Kanamori, 1991). It was developed later for regional and local distances by Zahradnik et al. (2005). In this method, Green’s functions are calculated by the discrete wavenumber method (Bouchon, 1981). The events used have moment magnitudes higher than 4 and encompass latitudes between 37° N and 40° N, and longitudes between 44° E and 49° E, during the period 2012–2014 for the focal mechanisms determined in this study and the period 1997–2014 for the Global Centroid Moment Tensor (GCMT) ones. We used the broadband stations of Iranian Seismological Center (IRSC), International Institute of Earthquake Engineering and Seismology (IIEES) and also stations of several other countries bordering northwestern of Iran. The focal mechanisms determined are often strike-slip and reverse which show a correspondence with the tectonic of this region. Then we analyze the stress state using these focal mechanisms and also by the focal mechanisms of other large and moderate earthquakes determined by GCMT in the region. We calculate the principal orientations of stresses by multiple inverse method that was originally proposed by Yamaji (2000). The method is a numerical technique to separate stresses from heterogeneous fault slip and focal mechanism data. Using the information of the acquired and GCMT focal mechanisms containing strike, dip and rake angles for the main shocks, we study the state of stress in the northwest of Iran. The result shows the average stress model with 1σ and 3σ equal to 141 and 50.2 degrees with a stress ratio of 0.6 in the region between Urmia Lake and Talesh in the northwest of Iran. This stress ratio shows that most of the motions are strike-slip. Using focal mechanisms of aftershocks, the same values are respectively 132.5, 42.4 degrees and 0.3 in Ahar–Varzeghan region. The value 0.3 shows that the motions are reverse in this part of the northwest of Iran. The small difference between these values in the northwest of Iran and Ahar–Varzeghan region shows that the values of the acquired principal stress directions using aftershocks are close to those of the main shocks. The faults directions of right-lateral strike-slip motion are in accordance with the stress direction determined.

Ground penetrating radar (GPR) method as a high-resolution non-destructive geophysical method acts based on the transmission of relatively high-frequency electromagnetic waves inside the ground and recording the reflected waves from the interfaces between the subsurface layers. As the method uses electromagnetic waves in the frequency ranges of 12.5 megahertz to 1 gigahertz (called GPR waves), it can only be used for shallow subsurface investigations. As a result of fast and dense data measurements in this method, continuous images of the reflections of GPR waves from the interfaces of subsurface media with different electrical or electromagnetic properties are obtained. These properties comprising of dielectric constant (or relative permittivity), electrical conductivity and magnetic permeability play key roles in GPR responses. The GPR equipment measures the travel time of the waves. Thus, the preliminary display of the acquired GPR data is in the form of a time section in which the vertical axis indicates the two-way time taken from the transmission of the GPR wave by the transmitter into the ground until its reflection and receipt by the receiver. GPR method has been successfully used in a variety of applications including hydrogeological investigations, mapping of bedrock surfaces, and detection of subsurface targets such as buried pipes, cavities, foundations, subsurface contaminations, waste deposits, water tables, soil horizons and other subsurface interfaces or targets. In this research, GPR method has been used to examine the subsurface sediments in southeast margin of the Caspian Sea. After acquiring the GPR data along a large number of relatively long survey lines in the study area, effort has been made to apply various processing techniques to the acquired GPR data in order to investigate the effect of each of processing techniques on the data enhancement. Considering the collection of vast GPR datasets along different long survey lines in the study area containing various subsurface targets with different depths and sizes, the results or performances of applying each of the processing techniques to the GPR data have not been similar. Due to the low distance between the GPR transmitter and receiver as well as the electrical properties, especially the conductivity of the ground, and also, to remove the unwanted low-frequency signals or reflections while preserving the highfrequency signals, the dewow filter has been applied before any other processing to all the GPR datasets. The short time intervals between the transmitted GPR pulses and the pulses received directly from the air-ground surface, and also, the existence of reflections from the shallow subsurface targets, cause signal saturation in the receiver. For this, the dewow filter is applied to the GPR data to correct for signal saturation or wow in the data. Different types of gain are also among the processing methods applied to the data to reduce the attenuating effect of the GPR waves as the depth increases. To demonstrate the effects of different gains and to select the optimum gain, we have applied different gains to the GPR data. To convert the trace from a wavelet with both positive and negative components (i.e., sine or cosine nature) to a monopulse wavelet with all positives, we have used the envelope filter. This process removes the oscillatory nature of the radar wavelet and shows the data in its true resolution, making it easier to interpret. To determine the locations of the GPR events, the GPR time section should be converted to its corresponding GPR depth section in which the vertical axis shows the depth. To do this, it is necessary to know the velocity of the GPR wave in the subsurface structures of the area under study. This research indicates that using the characteristics of GPR waves in the GPR sections, we can detect the subsurface targets and discriminate the coarse-grained sediments from the fine-grained sediments, and also determine the electrical properties of subsurface layers with high success. High resolution of the GPR data has enabled us to characterize most of the subsurface sediments. Furthermore, the shallow subsurface bedding can be easily observed in the GPR sections obtained. High moisture, salinity, clay and silt contents of the shallow subsurface sediments cause high conductivity of the ground in the area, and thus, cause the depth of penetration of the GPR waves to be mostly less than 1 meter.

For modeling the blast wave, its Green’s function must be solved first. Therefore, we use the Green’s function solution of explosive sources by Herman (1979). The study area was located around the gas pipes in high-speed rail project of Qom to Isfahan, from 50. Eto51.0Eand 33.5N to 34.5N. In this study, we used the four-digit seismograph machines of type CMG-6TD. The seismograph machines were arranged around the blast site in such a way that the three-component seismic energy radiation patterns of the blast were recorded in the vertical, radial and transverse directions. Further, a component (the radial component) was always placed in the site of the explosion. The distance between two consecutive samples was determined to be 10 milliseconds. The purpose of this arrangement of seismic machines has been to determine the blast radiation pattern and the impact of maximum velocity of the particles resulting from the explosion on a gas pipe. To study changes of the velocity model in output waveforms, three independent files are prepared, including the velocity model parameters of the region and the physical parameters of the source. We will generate three different outputs and select the output which is more consistent with the real waveform. The process is repeated by increasing and decreasing the source depth by 0.5 m and the amount of explosive materials by 0.5 kg (to impact the energy released). Finally, the best model of waveform, in terms of the closeness to reality, will be determined as a result. After obtaining the waveforms from the modeling and matching with real waveforms, to compare the accuracy of the model waveform with the real waveform the following methods were used. In the first method, with the calculation of the maximum normal correlation amplitude between the two series of models and real data, the difference of the resulting numbers from unity indicates which waveform is better and more accurate in modeling. In this method, the smaller the number is, the higher is the matching of that model with the real data. In the second method, the values of the natural logarithm (Ln) of each waveform and their differences are calculated. The positive difference of the standard deviation of the resulting values from one is considered as a measure of error. According to the second method, a high amount of epicentral distance leads to more consistency between models and real data. This results in blasts recorded from the station No. 4 to be clearly visible. Also, data from the station No. 2 in both methods show the lowest percentage of errors. The effects of the three factors of velocity model, the source depth and the energy of explosion were studied. The results obtained can be summarized as follows. For the velocity model, by decreasing the body wave velocities by 0.5kms−1, the correlation between the waves amplitudes diminishes in such a way that the model wave amplitudes increases with respect to the real wave amplitudes, and vice versa. For the source depth, the source depth alteration is directly related to the amplitudes of the vertical component of the waves, so that a 0.5 m increase in depth cause the amplitude to increase and decrease in depth causes the amplitude to decrease. To study the effect of the amount of the energy released on waveform modeling, having epicentral distance and maximum wave amplitudes, the equivalent magnitude can be calculated in local scale (ml) and by changing that, the released energy can be indirectly incorporated in modeling.