University of Tehran Press Journal of the Earth and Space Physics 2538-371X 40 1 2014 04 21 3D Tomography of Isfahan and Sharekord regions using local earthquakes 3D Tomography of Isfahan and Sharekord regions using local earthquakes 1 14 36692 10.22059/jesphys.2014.36692 FA Kourosh Emami M.Sc. Student, Earth Physics Department, Institute of Geophysics, University of Tehran, Iran Esmaeil Bayramnejad Assistant Professor, Earth Physics Department, Institute of Geophysics, University of Tehran, Iran Mohammad Reza Gheitanchi Professor, Earth Physics Department, Institute of Geophysics, University of Tehran, Iran Journal Article 2013 07 14 ^*نگارنده رابط:          تلفن: 61118230-021         دورنگار: 88630479-021                               E-mail: ebayram@ut.ac.ir   One of the most important purposes in seismology is determination of crustal velocity using earthquakes data that have been recorded by regional and local seismic stations. The more precise the process is carried out, the better shall be the results reached in various studies in the area including earthquake locating, seismicity of area, seismic zone mapping or the determination of the plane of faults causing earthquakes Imaging velocity structure of Earth's interior using travel times inversion commonly called seismic tomography which is usually done two or three-dimensionally. . This method has extensively been used in recent decades by researchers. Since the seismic waves are associated with valuable information about direction of propagation and environmental properties, using earthquakes as natural seismic sources are very useful in seismic tomography as well as the artificial sources such as limited and controlled explosions, air guns and bore-hole sources.. The seismic tomography characterizes the size, geometry and extent of velocity anomalies. In this method subsurface structure is modeled initially by several parameters and improved by the inversion of seismic travel times data. In this study the crustal velocity structure was determined using three-dimensional inversion of local earthquakes travel times recorded by seismic networks of Institute of Geophysics University of Tehran (IGUT) occurred within the period from 2000 and 2012 in the study area. The study area is bounded within 31^oE to 34^oE and 50^oN to 53^oN. We used the VELEST software in one-dimensional modeling section. The procedure of study is minimizing the differences between observed and calculated travel time by applying the initial obtained model. This software simultaneously optimizes the earthquake locations, crustal velocity model and station corrections using the Joint-Hypocenter-Determination (JHD) method.  The initial estimates for P waves velocity and crust thickness of the region are achieved using the travel-time curve of primary phases of all earthquakes occurred in the area. Then the larger relative earthquakes are selected and the best crustal one-dimensional model was derived by simultaneously inverse modeling method using this data set and VELEST algorithm. This method can be considered as one of the useful methods in study of 1D crust structure. The results proposed a 4 layers model of crust in which the P wave velocity is equal to 5.4 km/s for depths less than 5 km,   6.0 km/s for depths from 5 km to 20 km,  6.2 km/s for depths of 20 km to 32 km and finally  6.9 km/s for depths of 32 km to 47 km. The thickness of crust and P_nvelocity are respectively obtained 47 km and 7.9 km/s.  The aim of this work is obtaining an optimal crust model that can aid to improve seismic data and can be used to determine the next earthquake locating. Then the obtained crustal model is used as an initial model to study of three-dimensional inverse modeling of crust in the region by using FAST algorithm. All the earthquakes relocated using new obtained model. In this study a data set recorded by the 8 seismic stations of Isfahan and Shahrekord networks were used. The resolution of final solution of 3D model was investigated using synthetic dataset (checkerboard model) that shows fair resolution for various depths. The lateral variations of the main resolved structures in the model obtained are highly correlated with the faulting systems in the region. ^*نگارنده رابط:          تلفن: 61118230-021         دورنگار: 88630479-021                               E-mail: ebayram@ut.ac.ir   One of the most important purposes in seismology is determination of crustal velocity using earthquakes data that have been recorded by regional and local seismic stations. The more precise the process is carried out, the better shall be the results reached in various studies in the area including earthquake locating, seismicity of area, seismic zone mapping or the determination of the plane of faults causing earthquakes Imaging velocity structure of Earth's interior using travel times inversion commonly called seismic tomography which is usually done two or three-dimensionally. . This method has extensively been used in recent decades by researchers. Since the seismic waves are associated with valuable information about direction of propagation and environmental properties, using earthquakes as natural seismic sources are very useful in seismic tomography as well as the artificial sources such as limited and controlled explosions, air guns and bore-hole sources.. The seismic tomography characterizes the size, geometry and extent of velocity anomalies. In this method subsurface structure is modeled initially by several parameters and improved by the inversion of seismic travel times data. In this study the crustal velocity structure was determined using three-dimensional inversion of local earthquakes travel times recorded by seismic networks of Institute of Geophysics University of Tehran (IGUT) occurred within the period from 2000 and 2012 in the study area. The study area is bounded within 31^oE to 34^oE and 50^oN to 53^oN. We used the VELEST software in one-dimensional modeling section. The procedure of study is minimizing the differences between observed and calculated travel time by applying the initial obtained model. This software simultaneously optimizes the earthquake locations, crustal velocity model and station corrections using the Joint-Hypocenter-Determination (JHD) method.  The initial estimates for P waves velocity and crust thickness of the region are achieved using the travel-time curve of primary phases of all earthquakes occurred in the area. Then the larger relative earthquakes are selected and the best crustal one-dimensional model was derived by simultaneously inverse modeling method using this data set and VELEST algorithm. This method can be considered as one of the useful methods in study of 1D crust structure. The results proposed a 4 layers model of crust in which the P wave velocity is equal to 5.4 km/s for depths less than 5 km,   6.0 km/s for depths from 5 km to 20 km,  6.2 km/s for depths of 20 km to 32 km and finally  6.9 km/s for depths of 32 km to 47 km. The thickness of crust and P_nvelocity are respectively obtained 47 km and 7.9 km/s.  The aim of this work is obtaining an optimal crust model that can aid to improve seismic data and can be used to determine the next earthquake locating. Then the obtained crustal model is used as an initial model to study of three-dimensional inverse modeling of crust in the region by using FAST algorithm. All the earthquakes relocated using new obtained model. In this study a data set recorded by the 8 seismic stations of Isfahan and Shahrekord networks were used. The resolution of final solution of 3D model was investigated using synthetic dataset (checkerboard model) that shows fair resolution for various depths. The lateral variations of the main resolved structures in the model obtained are highly correlated with the faulting systems in the region. https://jesphys.ut.ac.ir/article_36692_d0618c769cba518c8d8a23cd7db7d818.pdf

University of Tehran Press Journal of the Earth and Space Physics 2538-371X 40 1 2014 04 21 Estimation of the kinematic source parameters and frequency independent shear wave Quality factor from acceleration records of the Ahar-Varzagan earthquake 2012 Estimation of the kinematic source parameters and frequency independent shear wave Quality factor from acceleration records of the Ahar-Varzagan earthquake 2012 15 27 36693 10.22059/jesphys.2014.36693 FA Neda Masominia M.Sc. Student of Geophysics, Earth Physics Department, Institute of Geophysics, University of Tehran, Iran Habib Rahimi Assistant Professor, Earth Physics Department, Institute of Geophysics, University of Tehran, Iran Mahdi Rezapour Associate Professor, Earth Physics Department, Institute of Geophysics, University of Tehran, Iran Journal Article 2013 08 18 ^*نگارنده رابط:          تلفن: 61118240-021         دورنگار: 88630479-021                                E-mail: rahimih@ut.ac.ir   Two relatively large earthquakes on the 11^th of August of 2012, struck the region of Varzagan and Ahar County, East Azerbaijan Province, NW Iran. The devastation caused by earthquakes in different regions and prediction of strong ground motions from large earthquakes has attracted the attention of seismologists to investigate the attenuation characteristics of the region and source characteristics of earthquakes to better understand the seismic hazards in different regions. Several recent and historical catastrophic earthquakes have destroyed different parts of East Azerbaijan Province. The seismic hazard map of Iran [published by the Building and Housing Research Center (BHRC)] shows that most of the cities in this Province are located within the high or very high relative risk areas. The quality of manmade constructions, especially in small villages and towns, is usually poor. Residences are generally built without considering seismic design regulations, so they are highly vulnerable and will collapse under shaking caused by even moderate earthquakes. A displacement spectrum contains valuable information regarding the source and medium characteristics. The source spectrum of an earthquake can be approximated by the omega-square model (Brune, 1970), which has ω² decay of high frequencies above the corner frequency. The source displacement spectrum can be estimated from a displacement record after correcting with diminution function, which accounts for the geometrical spreading and anelastic attenuation. The anelastic attenuation of seismic waves is characterized by a dimensionless quantity called quality factor Q (Knopoff, 1964). So far a few studies have been carried out to understand the attenuation characteristics of the Iranian crust. Examples include the work by Nutlii (1980), Michell (1995) and Rahimi & Hamzehloo (2008). An analysis scheme for obtaining source parameters and quality factor Q using the generalized inversion has been presented in this paper. The work presented here is approximately based on the technique of Fletcher (1995) and Joshi (2006) that used inversion methods. In this paper, the Brune’s source model (Brune, 1970) is used together with the propagation filter. This study uses the acceleration data of the Ahar-Varzagan main shock recorded by Building and Housing Research Center (BHRC) strong ground motion network. Our main objectives are: (i) to compute the source parameters of these earthquakes using the acceleration data, and (ii) to compute the frequency-independent shear wave quality factor in the recorded stations. In this study to get shear wave quality factor and source parameter in the near field, we used the strong motion data of the two earthquakes. Source parameters are estimated from that are related to two events with moment magnitudes of 6.5 and 6.4 in the hypocentral distance range from 22 to 206 Km. In this approach the theoretical S-wave displacement spectra, conditioned by frequency-independent Q, was fitted with the observed displacement spectra. Therefore corner frequency, moment magnitude and frequency-independent Q for each record are estimated simultaneously and the error estimate is given as the root-mean-square over all the frequencies. The source terms estimated here are , 3.26E+25(dyn-cm), , 0.17(Hz), , 8.19(km), 86.74, =32.02(bar), ∆ =110.4, ∆ =50.39(cm), =9.13, =5.51(sec)) and estimated moment magnitude ( , = 6.2) agree well with values obtained from telesiesmic wave of Harvard University. Estimate of path-average crustal shear –wave quality factors give a range of Q= 71 to 501 for frequency band of 0.01 to 15 Hz, that for near stations has a low value (high absorption) and for others at the further distance it has high value (low absorption), which shows good agreement with high-frequency absorption in near field. Independent estimates of Q at various stations give its average value of 276. ^*نگارنده رابط:          تلفن: 61118240-021         دورنگار: 88630479-021                                E-mail: rahimih@ut.ac.ir   Two relatively large earthquakes on the 11^th of August of 2012, struck the region of Varzagan and Ahar County, East Azerbaijan Province, NW Iran. The devastation caused by earthquakes in different regions and prediction of strong ground motions from large earthquakes has attracted the attention of seismologists to investigate the attenuation characteristics of the region and source characteristics of earthquakes to better understand the seismic hazards in different regions. Several recent and historical catastrophic earthquakes have destroyed different parts of East Azerbaijan Province. The seismic hazard map of Iran [published by the Building and Housing Research Center (BHRC)] shows that most of the cities in this Province are located within the high or very high relative risk areas. The quality of manmade constructions, especially in small villages and towns, is usually poor. Residences are generally built without considering seismic design regulations, so they are highly vulnerable and will collapse under shaking caused by even moderate earthquakes. A displacement spectrum contains valuable information regarding the source and medium characteristics. The source spectrum of an earthquake can be approximated by the omega-square model (Brune, 1970), which has ω² decay of high frequencies above the corner frequency. The source displacement spectrum can be estimated from a displacement record after correcting with diminution function, which accounts for the geometrical spreading and anelastic attenuation. The anelastic attenuation of seismic waves is characterized by a dimensionless quantity called quality factor Q (Knopoff, 1964). So far a few studies have been carried out to understand the attenuation characteristics of the Iranian crust. Examples include the work by Nutlii (1980), Michell (1995) and Rahimi & Hamzehloo (2008). An analysis scheme for obtaining source parameters and quality factor Q using the generalized inversion has been presented in this paper. The work presented here is approximately based on the technique of Fletcher (1995) and Joshi (2006) that used inversion methods. In this paper, the Brune’s source model (Brune, 1970) is used together with the propagation filter. This study uses the acceleration data of the Ahar-Varzagan main shock recorded by Building and Housing Research Center (BHRC) strong ground motion network. Our main objectives are: (i) to compute the source parameters of these earthquakes using the acceleration data, and (ii) to compute the frequency-independent shear wave quality factor in the recorded stations. In this study to get shear wave quality factor and source parameter in the near field, we used the strong motion data of the two earthquakes. Source parameters are estimated from that are related to two events with moment magnitudes of 6.5 and 6.4 in the hypocentral distance range from 22 to 206 Km. In this approach the theoretical S-wave displacement spectra, conditioned by frequency-independent Q, was fitted with the observed displacement spectra. Therefore corner frequency, moment magnitude and frequency-independent Q for each record are estimated simultaneously and the error estimate is given as the root-mean-square over all the frequencies. The source terms estimated here are , 3.26E+25(dyn-cm), , 0.17(Hz), , 8.19(km), 86.74, =32.02(bar), ∆ =110.4, ∆ =50.39(cm), =9.13, =5.51(sec)) and estimated moment magnitude ( , = 6.2) agree well with values obtained from telesiesmic wave of Harvard University. Estimate of path-average crustal shear –wave quality factors give a range of Q= 71 to 501 for frequency band of 0.01 to 15 Hz, that for near stations has a low value (high absorption) and for others at the further distance it has high value (low absorption), which shows good agreement with high-frequency absorption in near field. Independent estimates of Q at various stations give its average value of 276. https://jesphys.ut.ac.ir/article_36693_36eb49f8bb1c04379dd4ac3ae0f9d9ee.pdf

University of Tehran Press Journal of the Earth and Space Physics 2538-371X 40 1 2014 04 21 Static correction using blind channel identification Static correction using blind channel identification 29 43 36694 10.22059/jesphys.2014.36694 FA Seyed Hosein Seyed Aghamiry M.Sc. in Geophysics, Physics Department, Faculty of Dezfoul, Technical and Vocational Univercity, Iran Hamid Reza Siahkoohi Professor, Earth Physics Department, Institute of Geophysics, University of Tehran, Iran Journal Article 2011 09 28 Near-surface variations can be very complex and may distort amplitudes and arrival times of the reflections events from target reflectors. Near-surface complexities include topographic variations, near-surface irregularities, variations in soil conditions and the weathered layer. These perturbations generally have a significant impact on seismic recordings. Although there is a general agreement that near-surface distortions are very complex and we usually rely on a rather simplified parameterization to compensate for these perturbations. Determinatin of time shifts is generally referred to as static corrections and residual static correction. Underlying concept of static corrections is the assumption that a simple time shift of an entire seismic trace will yield the seismic record that would have been observed if the geophones had been placed on the reference datum. Hence, static time shifts corrections are assumed to be surface consistent. Surface consistencymeans that the effects associated with a particular source or receiver affect all wave types similarly, regardless of the direction of propagation. Conventional methods of static time shift corrections need information on velocities and depths of near-surface layers to determine and compensate the time shifts. These methods rely on simple models for near-surface layers. In this paper, we develop an approach to compensate for complex time shift using blind channel identification, as it does not use near-surface information. The blind channel identification deals with the recovery of either the input signal or the channel response from the observed transmitted signal only. This method differs from conventional methods for seismic deconvolution. The latter resolve the undetermined nature of the problem by making assumptions about the reflectivity sequence (whiteness, sparsity) and/or the seismic wavelet (minimum phase/ zero phase). The blind channel identification method does not rely on these assumptions. It uses multichannel recordings to fully constrain the problem and is therefore purely data driven. Many recent blind channel estimation techniques exploitsubspace structures of observation. The key idea in subspace methods of blind channel identification that the channel vector (or part of the channel vector) is in a one dimensional subspace of a block of noiseless observations. These methods, which are often referred to as subspace algorithms, have the attractive property that the channel estimates can often be obtained in a closed form from optimizing a quadratic cost function. We use blind channel identification to estimate for near-surface source and receiver perturbations. These perturbations are parameterized as finite-impulse response (FIR) filters, and are referred to as the channels. Because the channels describe the near-surface perturbations, we can estimate time shifts from correlation of the channels. We applied the method to synthetic data and to part of a field data set acquired in an area with significant near-surface heterogeneity. The application of new static corrections greatly improves the trace-to-trace consistency in prestack data. The procedure delineates reflection events that are difficult to detect prior to the application of new static corrections. Based on these results, we conclude that the new static corrections can successfully remove complex time shifts from land seismic data. The field data example demonstrates that the new static corrections can greatly enhance the imaging capabilities of land seismic data.     Near-surface variations can be very complex and may distort amplitudes and arrival times of the reflections events from target reflectors. Near-surface complexities include topographic variations, near-surface irregularities, variations in soil conditions and the weathered layer. These perturbations generally have a significant impact on seismic recordings. Although there is a general agreement that near-surface distortions are very complex and we usually rely on a rather simplified parameterization to compensate for these perturbations. Determinatin of time shifts is generally referred to as static corrections and residual static correction. Underlying concept of static corrections is the assumption that a simple time shift of an entire seismic trace will yield the seismic record that would have been observed if the geophones had been placed on the reference datum. Hence, static time shifts corrections are assumed to be surface consistent. Surface consistencymeans that the effects associated with a particular source or receiver affect all wave types similarly, regardless of the direction of propagation. Conventional methods of static time shift corrections need information on velocities and depths of near-surface layers to determine and compensate the time shifts. These methods rely on simple models for near-surface layers. In this paper, we develop an approach to compensate for complex time shift using blind channel identification, as it does not use near-surface information. The blind channel identification deals with the recovery of either the input signal or the channel response from the observed transmitted signal only. This method differs from conventional methods for seismic deconvolution. The latter resolve the undetermined nature of the problem by making assumptions about the reflectivity sequence (whiteness, sparsity) and/or the seismic wavelet (minimum phase/ zero phase). The blind channel identification method does not rely on these assumptions. It uses multichannel recordings to fully constrain the problem and is therefore purely data driven. Many recent blind channel estimation techniques exploitsubspace structures of observation. The key idea in subspace methods of blind channel identification that the channel vector (or part of the channel vector) is in a one dimensional subspace of a block of noiseless observations. These methods, which are often referred to as subspace algorithms, have the attractive property that the channel estimates can often be obtained in a closed form from optimizing a quadratic cost function. We use blind channel identification to estimate for near-surface source and receiver perturbations. These perturbations are parameterized as finite-impulse response (FIR) filters, and are referred to as the channels. Because the channels describe the near-surface perturbations, we can estimate time shifts from correlation of the channels. We applied the method to synthetic data and to part of a field data set acquired in an area with significant near-surface heterogeneity. The application of new static corrections greatly improves the trace-to-trace consistency in prestack data. The procedure delineates reflection events that are difficult to detect prior to the application of new static corrections. Based on these results, we conclude that the new static corrections can successfully remove complex time shifts from land seismic data. The field data example demonstrates that the new static corrections can greatly enhance the imaging capabilities of land seismic data.     https://jesphys.ut.ac.ir/article_36694_2e8f8fdffba0ebb520aaa5d037f3c30e.pdf

University of Tehran Press Journal of the Earth and Space Physics 2538-371X 40 1 2014 04 21 Demonstrating buried channels using principal component analysis Demonstrating buried channels using principal component analysis 45 56 36695 10.22059/jesphys.2014.36695 FA Mahdi Sadeghi M.Sc. Student, School of Mining, Petroleum and Geophysics, Shahrood University of Technology, Iran Amin Roshandel Assistant Professor, School of Mining, Petroleum and Geophysics, Shahrood University of Technology, Iran Hamid Reza Siahkoohi Professor, Earth Physics Department, Institute of Geophysics, University of Tehran, Iran Alireza Heydarian Ph.D. Student, Institute of Geophysics, University of Tehran, Iran Journal Article 2012 08 05 Spectral decomposition of time series has a significant role in seismic data processing and interpretation. Since the earth acts as a low-pass filter, it changes frequency content of passing seismic waves. Conventional representing methods of signals in time domain and frequency domain cannot show time and frequency information simultaneously. Time-frequency transforms upgraded spectral decomposition to a new level and can show time and frequency information simultaneously. Time-frequency transforms generate high volume of spectral components, which contain useful information about the reservoir and can be decomposed into single frequency volumes. These single frequency volumes can overload the limited space of computer hard disk and are not easy for an interpreter to investigate them individually; therefore, it is important to use methods to decrease volume with no information lost, so frequency slices are separated from these volumes and used for interpretation. An expert interpreter can achieve some information about channel content and lateral variation is of thickness. Since different frequencies contain different types of information (low frequencies are sensible to channel content and high frequencies are sensible to channel boundaries), these slices cannot show this information simultaneously. Therefore RGB images can be produced by plotting three different frequency slices against red, green and blue components. An RGB image, sometimes referred to as a true color image, it is an image that defines red, green, and blue color components for each individual pixel and has intensity between 0 and 1. Although this method obviates some drawbacks of single frequency plots, but it uses only three slices and practically ignores a big part of information and the frequency choice is not clear, so different choices will result to different images. Principal component is a statistical method for identifying patterns in data and expressing them in a way to highlight their similarities and differences. In order to find major patterns in data this technique reduces the number of dimensions of data without the loss of information. Principal component analysis introduces new set of orthogonal axes through data set called “eigenvectors” which data variance along them is maximized and have the importance proportional to their corresponding eigenvalues. The projection of single frequency slices onto eigenvectors is called “principal component (PC) bands”. The amount of total variance that each PC band represents is proportional to its eigenvalue, thus after normalizing the total sum of all eigenvalues, each eigenvalue represents the percentage of total spectral variance that its corresponding principal component can represent. So the first PC band (having largest eigenvalue) best represents the spectral variance in data, the second PC band (having the second largest eigenvalue) best represents the spectral variance in data which is not represented by the first PC and so on. Therefore PC bands with the smallest eigenvalues will represent a small portion of variance and can be deduced as random noise. So choosing the PC bands with the largest eigenvalues can be an effective way for data denoising, image processing and in our case determining the major trends in data set. We can represent more than 80 percent of spectral variation by plotting three largest principal components against red, green and blue components in a RGB image. In this paper, we applied spectral decomposition on land seismic data of an oil field in south-west of Iran using short time Fourier transform (STFT) and S transform. Then we constructed single frequency slices and investigated them. We produced RGB images by color stacking method and improved interpretation. Finally we used principal component analysis to use all the frequency bandwidth. Our results showed that PCA based images showed channel and its branches in a more precise manner than the other methods.     Spectral decomposition of time series has a significant role in seismic data processing and interpretation. Since the earth acts as a low-pass filter, it changes frequency content of passing seismic waves. Conventional representing methods of signals in time domain and frequency domain cannot show time and frequency information simultaneously. Time-frequency transforms upgraded spectral decomposition to a new level and can show time and frequency information simultaneously. Time-frequency transforms generate high volume of spectral components, which contain useful information about the reservoir and can be decomposed into single frequency volumes. These single frequency volumes can overload the limited space of computer hard disk and are not easy for an interpreter to investigate them individually; therefore, it is important to use methods to decrease volume with no information lost, so frequency slices are separated from these volumes and used for interpretation. An expert interpreter can achieve some information about channel content and lateral variation is of thickness. Since different frequencies contain different types of information (low frequencies are sensible to channel content and high frequencies are sensible to channel boundaries), these slices cannot show this information simultaneously. Therefore RGB images can be produced by plotting three different frequency slices against red, green and blue components. An RGB image, sometimes referred to as a true color image, it is an image that defines red, green, and blue color components for each individual pixel and has intensity between 0 and 1. Although this method obviates some drawbacks of single frequency plots, but it uses only three slices and practically ignores a big part of information and the frequency choice is not clear, so different choices will result to different images. Principal component is a statistical method for identifying patterns in data and expressing them in a way to highlight their similarities and differences. In order to find major patterns in data this technique reduces the number of dimensions of data without the loss of information. Principal component analysis introduces new set of orthogonal axes through data set called “eigenvectors” which data variance along them is maximized and have the importance proportional to their corresponding eigenvalues. The projection of single frequency slices onto eigenvectors is called “principal component (PC) bands”. The amount of total variance that each PC band represents is proportional to its eigenvalue, thus after normalizing the total sum of all eigenvalues, each eigenvalue represents the percentage of total spectral variance that its corresponding principal component can represent. So the first PC band (having largest eigenvalue) best represents the spectral variance in data, the second PC band (having the second largest eigenvalue) best represents the spectral variance in data which is not represented by the first PC and so on. Therefore PC bands with the smallest eigenvalues will represent a small portion of variance and can be deduced as random noise. So choosing the PC bands with the largest eigenvalues can be an effective way for data denoising, image processing and in our case determining the major trends in data set. We can represent more than 80 percent of spectral variation by plotting three largest principal components against red, green and blue components in a RGB image. In this paper, we applied spectral decomposition on land seismic data of an oil field in south-west of Iran using short time Fourier transform (STFT) and S transform. Then we constructed single frequency slices and investigated them. We produced RGB images by color stacking method and improved interpretation. Finally we used principal component analysis to use all the frequency bandwidth. Our results showed that PCA based images showed channel and its branches in a more precise manner than the other methods.     https://jesphys.ut.ac.ir/article_36695_43f2090942a06f395ee54a7985d74124.pdf

University of Tehran Press Journal of the Earth and Space Physics 2538-371X 40 1 2014 04 21 Qualitative characterization of core samples porosity using (MRI) image (a carbonate reservoir in southern Iran) Qualitative characterization of core samples porosity using (MRI) image (a carbonate reservoir in southern Iran) 57 67 36696 10.22059/jesphys.2014.36696 FA Ezatollah Kazemzadeh Assistant Professor, Research Institute of Petroleum Industry (RIPI), Tehran, Iran Amir Salehi Masters Student, Islamic Azad University - Science and Research Branch, Tehran, Iran Seyed Jamal Sheikhzakariai Assistant Professor, Islamic Azad University - Science and Research Branch, Tehran, Iran Armin Afrough Research Assistant, Research Institute of Petroleum Industry (RIPI), Tehran, Iran Journal Article 2012 10 15 Qualitative analysis of petroleum reservoir rocks is still one of the most important topics of core laboratories and directly affects hydrocarbon-in-place, fluid flow, and prediction of field performance. Magnetic Resonance Imaging (MRI) as a non-invasive, millimeter resolution technique which only images fluids in porous media fits this purpose perfectly. Magnetic resonance is a radio frequency spectrometry, based on excitation of nucleus energy levels which can exploit a wide range of information on the saturation fluid, geometry of pores, and diffusion. Using magnetic gradients and signal encoding, this can be used as a tomography technique. Determining porosity regardless of rock lithology, reservoir rock quality, bound and free water (which presents production potential) and potential permeable tight beds are applications of this method, from among many others. Recent advances in imaging techniques along with new software and processing methods has resulted in exploiting the images containing valuable physical information. The main aims of this study are qualitative investigation of core sample, determining the number and distribution of vugs, presenting three-dimensional models and comparing this method with other conventional methods. We acquired images of cores adequate for revealing characteristics of matrix, vugs, and other different porosity types in carbonate rocks as well as their interaction. The images in the form of matrices of MR signal were analyzed using both image analysis and physics of MR signal. Studied samples are a selection from one of southern Iranian carbonate reservoirs, cleaned using Soxhlet extraction, dried in oven and saturated with synthetic reservoir brine. Rock type, pore-filling fluid, the MRI imaging hardware and software, pulse sequence, and image processing affect such study results. So, impacts of several items were considered and best available pulse sequence parameters were set. Presence of ferromagnetic minerals adversely affects image quality, so samples bearing these minerals are very difficult to image in high fields and have to be imaged with special pulse sequences and systems. Samples without ferromagnetic minerals, as in most carbonate rocks, can be imaged in high fields, so because of the superior quality of high field imaging, this method is used in our carbonate samples. Validity of MRI images was verified using histogram analysis of water, rock and air segments of the images subsequent to acquisition. Reference fluids, brine (the pore-filling fluid) and air, helped us in matching and comparing histograms and check the validity of signal from porous sample. In image analysis we utilized histogram, field of view, and segmentation techniques. Results of this analysis after using physical models of MRI signal in porous media led to numerical and visual models of rock samples. In addition to visual models of porosity, we prepared visual models of mean-T2 of invaluable to quantitative and qualitative study of porosity, ultimately resulting in determining fraction of vuggy, moldic, and inter-particle porosity. This model was constructed by superposition of image slices. Inter-particle porosity cannot be determined from sub-millimeter MRI image analysis, so it is calculated from the physics of MR signal. We determined the accuracy of the method in comparison with other conventional experiments such as helium porosimetry and petrographic image analysis which revealed reasonable accuracy of the method in determining porosity types and visualization. Effect of several items in the accuracy of this method is proposed. First, gravimetry porosity is always lower compared to the helium porosimetry, and porosity calculated from MRI imaging only is affected by water saturation of the sample, so MRI porosity should always be compared with gravimetry porosity. Second, MRI images are not only sensitive to water in porous samples, but also sensitive to pore size of the sample. Water in very fine pores is not shown in images unless echo time is set to 1-3ms, not possible in our study because of hardware limits. Third, magnetic field inhomogeneity can account for up to nine percent of signal intensity. This study obtained a new perspective in using MRI imaging for qualitative study of porosity and determining share of porosity types.       Qualitative analysis of petroleum reservoir rocks is still one of the most important topics of core laboratories and directly affects hydrocarbon-in-place, fluid flow, and prediction of field performance. Magnetic Resonance Imaging (MRI) as a non-invasive, millimeter resolution technique which only images fluids in porous media fits this purpose perfectly. Magnetic resonance is a radio frequency spectrometry, based on excitation of nucleus energy levels which can exploit a wide range of information on the saturation fluid, geometry of pores, and diffusion. Using magnetic gradients and signal encoding, this can be used as a tomography technique. Determining porosity regardless of rock lithology, reservoir rock quality, bound and free water (which presents production potential) and potential permeable tight beds are applications of this method, from among many others. Recent advances in imaging techniques along with new software and processing methods has resulted in exploiting the images containing valuable physical information. The main aims of this study are qualitative investigation of core sample, determining the number and distribution of vugs, presenting three-dimensional models and comparing this method with other conventional methods. We acquired images of cores adequate for revealing characteristics of matrix, vugs, and other different porosity types in carbonate rocks as well as their interaction. The images in the form of matrices of MR signal were analyzed using both image analysis and physics of MR signal. Studied samples are a selection from one of southern Iranian carbonate reservoirs, cleaned using Soxhlet extraction, dried in oven and saturated with synthetic reservoir brine. Rock type, pore-filling fluid, the MRI imaging hardware and software, pulse sequence, and image processing affect such study results. So, impacts of several items were considered and best available pulse sequence parameters were set. Presence of ferromagnetic minerals adversely affects image quality, so samples bearing these minerals are very difficult to image in high fields and have to be imaged with special pulse sequences and systems. Samples without ferromagnetic minerals, as in most carbonate rocks, can be imaged in high fields, so because of the superior quality of high field imaging, this method is used in our carbonate samples. Validity of MRI images was verified using histogram analysis of water, rock and air segments of the images subsequent to acquisition. Reference fluids, brine (the pore-filling fluid) and air, helped us in matching and comparing histograms and check the validity of signal from porous sample. In image analysis we utilized histogram, field of view, and segmentation techniques. Results of this analysis after using physical models of MRI signal in porous media led to numerical and visual models of rock samples. In addition to visual models of porosity, we prepared visual models of mean-T2 of invaluable to quantitative and qualitative study of porosity, ultimately resulting in determining fraction of vuggy, moldic, and inter-particle porosity. This model was constructed by superposition of image slices. Inter-particle porosity cannot be determined from sub-millimeter MRI image analysis, so it is calculated from the physics of MR signal. We determined the accuracy of the method in comparison with other conventional experiments such as helium porosimetry and petrographic image analysis which revealed reasonable accuracy of the method in determining porosity types and visualization. Effect of several items in the accuracy of this method is proposed. First, gravimetry porosity is always lower compared to the helium porosimetry, and porosity calculated from MRI imaging only is affected by water saturation of the sample, so MRI porosity should always be compared with gravimetry porosity. Second, MRI images are not only sensitive to water in porous samples, but also sensitive to pore size of the sample. Water in very fine pores is not shown in images unless echo time is set to 1-3ms, not possible in our study because of hardware limits. Third, magnetic field inhomogeneity can account for up to nine percent of signal intensity. This study obtained a new perspective in using MRI imaging for qualitative study of porosity and determining share of porosity types.       https://jesphys.ut.ac.ir/article_36696_7ddd00ed8d1809f07b426aa9346c5677.pdf

University of Tehran Press Journal of the Earth and Space Physics 2538-371X 40 1 2014 04 21 An approach to height datum unification based on local gravity field modeling using radial base function, case study: height datum unification of leveling network of class 1 in Iran An approach to height datum unification based on local gravity field modeling using radial base function, case study: height datum unification of leveling network of class 1 in Iran 69 81 36697 10.22059/jesphys.2014.36697 FA Abdolreza Safari Associate Professor, Department of Surveying and Geomatics Engineering, University College of Engineering, University of Tehran, Iran 0000-0001-5938-5468 Mohammad Ali Sharifi Assistant Professor, Department of Surveying and Geomatics Engineering, University College of Engineering, University of Tehran, Iran 0000-0003-0745-4147 Ismaeil Foroughi M.Sc. student of Geodesy, Department of Surveying and Geomatics Engineering, University College of Engineering, University of Tehran, Iran Hadi Amin M.Sc. student of Geodesy, Department of Surveying and Geomatics Engineering, University College of Engineering, University of Tehran, Iran Journal Article 2013 08 24 One of the most important problems in geodesy is the unification of height datum. Generally in geodesy; there are two types of height systems, the geometrical height based on ellipsoid and the physical height based on gravity-defined surface (Zhang et al, 2009).Local height datum is determined according to Mean Sea Level (MSL). In regarding to mismatch of mean sea level and geoid, on the one hand, and height datum unification requirement on the other hand, this paper defines an approach to height datum determination of Iran based on local gravity field modeling. Although there are so many algorithms to local gravity field modeling, the radial base functions (RBF) is one of the well known methods to precise local gravity field modeling using variety of data sets. Two groups of data sets, gravity acceleration observation and satellite altimetry data, are used in this paper, to determine the quasi geoid height at reference tide gauge; then the height reference of leveling network can be calculated. According to Runge-Kutla algorithm, potential anomaly on the earth can be defined as following:   Where the expansion coefficients are (scale coefficients) and Bjerhammar sphere is a sphere with radius R, which is entirely inside the topographic masses of the earth,  are the set of radial basis functions with following representation:   Where  are points inside and outside of the Bjerhammar respectively, is the Legendre polynomial function of degree  and  are the Legendre coefficients, the point <em>y </em>is called the centre of the RBF. As we can see Radial Base Functions (RBF) have some unknown parameters that should be determined precisely: the location of RBF center, shape (depth parameter) and scale coefficient. If these parameters are selected correctly we could have good representation of potential anomaly field. Examples of linear functionals used in local gravity field modeling are gravity anomalies and gravity disturbances . After linearization and spherical approximation, these functionals are related to the potential anomaly as                                                         We used Poisson-wavelet as RBF kernel, which is defined as follows:                                               Where n is the order of Poisson wavelet kernel,  is the operator norm. Significant points in this paper are: precise positioning of the reference bench mark point by GPS system, calculation of geoid height of reference bench mark using both geodetic and orthometric heights by leveling method, quasi geoid height determination using potential anomaly value at bench mark and convert it to geoid height and finally height datum unification using comparison of the leveling height and the height from local gravity field modeling.     One of the most important problems in geodesy is the unification of height datum. Generally in geodesy; there are two types of height systems, the geometrical height based on ellipsoid and the physical height based on gravity-defined surface (Zhang et al, 2009).Local height datum is determined according to Mean Sea Level (MSL). In regarding to mismatch of mean sea level and geoid, on the one hand, and height datum unification requirement on the other hand, this paper defines an approach to height datum determination of Iran based on local gravity field modeling. Although there are so many algorithms to local gravity field modeling, the radial base functions (RBF) is one of the well known methods to precise local gravity field modeling using variety of data sets. Two groups of data sets, gravity acceleration observation and satellite altimetry data, are used in this paper, to determine the quasi geoid height at reference tide gauge; then the height reference of leveling network can be calculated. According to Runge-Kutla algorithm, potential anomaly on the earth can be defined as following:   Where the expansion coefficients are (scale coefficients) and Bjerhammar sphere is a sphere with radius R, which is entirely inside the topographic masses of the earth,  are the set of radial basis functions with following representation:   Where  are points inside and outside of the Bjerhammar respectively, is the Legendre polynomial function of degree  and  are the Legendre coefficients, the point <em>y </em>is called the centre of the RBF. As we can see Radial Base Functions (RBF) have some unknown parameters that should be determined precisely: the location of RBF center, shape (depth parameter) and scale coefficient. If these parameters are selected correctly we could have good representation of potential anomaly field. Examples of linear functionals used in local gravity field modeling are gravity anomalies and gravity disturbances . After linearization and spherical approximation, these functionals are related to the potential anomaly as                                                         We used Poisson-wavelet as RBF kernel, which is defined as follows:                                               Where n is the order of Poisson wavelet kernel,  is the operator norm. Significant points in this paper are: precise positioning of the reference bench mark point by GPS system, calculation of geoid height of reference bench mark using both geodetic and orthometric heights by leveling method, quasi geoid height determination using potential anomaly value at bench mark and convert it to geoid height and finally height datum unification using comparison of the leveling height and the height from local gravity field modeling.     https://jesphys.ut.ac.ir/article_36697_0b656371d2990c2722855528807c16c2.pdf

University of Tehran Press Journal of the Earth and Space Physics 2538-371X 40 1 2014 04 21 Paleomagnetism of Damavand volcano over the last half million years Paleomagnetism of Damavand volcano over the last half million years 83 93 36698 10.22059/jesphys.2014.36698 FA Mohammad Hamedpour Darabi Assistant Professor, Physics Department, Hormozgan University, Iran Ali Akbar Bahmanizadeh M.Sc. Student, Earth Physics Department, Institute of Geophysics, University of Tehran, Iran Behrooz Oskooi Assistant Professor, Earth Physics Department, Institute of Geophysics, University of Tehran, Iran 0000-0003-3065-194X Journal Article 2011 09 21 Damavand volcano with 5671 m height from sea level which covers 400 square km of area is the highest peak in the Middle East. It is located at 50 km to the east of Capital city of Tehran, Iran. It is in the fumarolic stage and exemplifies ongoing Quaternary volcanic activity. Damavan volcano was influenced by regional tectonics during its history. Allenbach (1966) believed that the volcano originates from a fault that already existed in the sedimentary basin which allowed the magma to rise. Darvishzadeh (1985) believed that the compression movement of Iran plate that influenced the Alborz Mountains which commenced the two curved faults (Vara-rud and Ask) that joined under the cone of the volcano and let the magma rise to the earth's crust. Frequent eruptions during ~ 1.8 Ma of Damavand volcanos make the opportunity to study the Earth’s magnetic field records in its stacked lava to clear out tectonic movement of the cone. Davidson et al., (2004) reported three major phases of volcanic activity during the past 1.8 Ma. They present radiometric age dating (U-TH)/ He, for several samples from different lava deposits around the cone. The youngest deposits of lava are distributed on the western flank of the cone and most of them are Trachytes. On the southern flank of the cone in which lava deposits are of the Trachy-andesites, show the oldest age among the collected samples. We followed Davidson et al., 2004, sampling sites to use their dated ages for our study. We used a water supplied petrol drill to collect 200 oriented samples from 10 sites around the cone, from north-west to the north-east. All samples were thermally demagnetized and results were plotted on orthogonal diagrams and principal component analysis (PCA) was carried out to extract primary magnetic remnant components. Orthogonal diagrams of the sites D1- D10 mostly show two magnetic components of 100-400 ̊ ^C and 400-600 ̊ ^C. sites D10, D9, D8 and D6 also show a third hematite bearing component with curie temperature around 675 ̊ ^C. VRM components exist in very low temperatures with a maximum boundary of 300 ̊ ^C, they have omitted from the data list during component analysis. Magnetic susceptibility variations with temperature during thermal treatments show a stable mineral composition till 600 ̊ ^C therefore we have not included the components over 600 ̊ ^C. Magnetic directions and their palaeomagnetic pole positions were calculated. The calculated α95 for most of the mean directions were low therefore mean directions and consequently pole positions are reliable. Magnetic inclination and declination variations are plotted according to the Davidson <em>et al</em>., 2004, from old to the youngest ages. Our results show two Earth’s magnetic polarity anomalies in 7000 and 27000 years ago, which have not been shown in GPTS for Late Pleistocene and Holocene normal polarity chron (Cande & Kent, 1992). Magnetic inclinations for these two sites show reverse polarity however, magnetic inclination only for 27000 y.a. event have a 180̊ rotation and it has a rotation of 90̊ for the 7000 y.a. event. Therefore a reversal for the 27000 y.a. event and an excursion for the 7000 y.a. event is probable. However we don’t see such record in the previous works therefore these might be the effect of a non-dipole in this area in that time. We have used window method (Besse, & Courtillot, 2002) for averaging paleomagnetic poles in three ranges of age of 7.0-7.2 Ky , 25-65 Ky and 250- 400 Ky. Poles arrangement shows an anti-clockwise rotation of Damavand cone with an average of  7/0 ̊  per 1 Ky. Many faults which are distributed around the cone can be the evidences for such a rotation for the cone.     Damavand volcano with 5671 m height from sea level which covers 400 square km of area is the highest peak in the Middle East. It is located at 50 km to the east of Capital city of Tehran, Iran. It is in the fumarolic stage and exemplifies ongoing Quaternary volcanic activity. Damavan volcano was influenced by regional tectonics during its history. Allenbach (1966) believed that the volcano originates from a fault that already existed in the sedimentary basin which allowed the magma to rise. Darvishzadeh (1985) believed that the compression movement of Iran plate that influenced the Alborz Mountains which commenced the two curved faults (Vara-rud and Ask) that joined under the cone of the volcano and let the magma rise to the earth's crust. Frequent eruptions during ~ 1.8 Ma of Damavand volcanos make the opportunity to study the Earth’s magnetic field records in its stacked lava to clear out tectonic movement of the cone. Davidson et al., (2004) reported three major phases of volcanic activity during the past 1.8 Ma. They present radiometric age dating (U-TH)/ He, for several samples from different lava deposits around the cone. The youngest deposits of lava are distributed on the western flank of the cone and most of them are Trachytes. On the southern flank of the cone in which lava deposits are of the Trachy-andesites, show the oldest age among the collected samples. We followed Davidson et al., 2004, sampling sites to use their dated ages for our study. We used a water supplied petrol drill to collect 200 oriented samples from 10 sites around the cone, from north-west to the north-east. All samples were thermally demagnetized and results were plotted on orthogonal diagrams and principal component analysis (PCA) was carried out to extract primary magnetic remnant components. Orthogonal diagrams of the sites D1- D10 mostly show two magnetic components of 100-400 ̊ ^C and 400-600 ̊ ^C. sites D10, D9, D8 and D6 also show a third hematite bearing component with curie temperature around 675 ̊ ^C. VRM components exist in very low temperatures with a maximum boundary of 300 ̊ ^C, they have omitted from the data list during component analysis. Magnetic susceptibility variations with temperature during thermal treatments show a stable mineral composition till 600 ̊ ^C therefore we have not included the components over 600 ̊ ^C. Magnetic directions and their palaeomagnetic pole positions were calculated. The calculated α95 for most of the mean directions were low therefore mean directions and consequently pole positions are reliable. Magnetic inclination and declination variations are plotted according to the Davidson <em>et al</em>., 2004, from old to the youngest ages. Our results show two Earth’s magnetic polarity anomalies in 7000 and 27000 years ago, which have not been shown in GPTS for Late Pleistocene and Holocene normal polarity chron (Cande & Kent, 1992). Magnetic inclinations for these two sites show reverse polarity however, magnetic inclination only for 27000 y.a. event have a 180̊ rotation and it has a rotation of 90̊ for the 7000 y.a. event. Therefore a reversal for the 27000 y.a. event and an excursion for the 7000 y.a. event is probable. However we don’t see such record in the previous works therefore these might be the effect of a non-dipole in this area in that time. We have used window method (Besse, & Courtillot, 2002) for averaging paleomagnetic poles in three ranges of age of 7.0-7.2 Ky , 25-65 Ky and 250- 400 Ky. Poles arrangement shows an anti-clockwise rotation of Damavand cone with an average of  7/0 ̊  per 1 Ky. Many faults which are distributed around the cone can be the evidences for such a rotation for the cone.     https://jesphys.ut.ac.ir/article_36698_3039d41a24ba453f9e6ec2478497eee4.pdf

University of Tehran Press Journal of the Earth and Space Physics 2538-371X 40 1 2014 04 21 Integration of Geoelecrtical information layers by fuzzy method to choose the best point for drilling: a case study Hamyj, Birjand Integration of Geoelecrtical information layers by fuzzy method to choose the best point for drilling: a case study Hamyj, Birjand 95 105 36699 10.22059/jesphys.2014.36699 FA Mohammad Shahi Ferdows Ph.D. Student of Mining Exploration, Amir Kabir University, Tehran, Iran Mohammad Bicharanlo M.Sc. Student of Mining Exploration, Sahand University of Technology, Tabriz, Iran Rashed Poormirzaee Ph.D. Student of Mining Exploration, Sahand University of technology, Tabriz, Iran Journal Article 2013 02 25 Resistivity and induced polarization methods are used in exploration of porphyry metals for years. The resistivity method is used in the study of horizontal and vertical discontinuities in the electrical properties of the ground. The induced polarization method makes use of the capacitive action of the subsurface to locate zones where conductive minerals are disseminated within their host rocks. The simplicity of the equipment, the lower cost of the survey compared to the other methods and the abundance of interpretation methods make it a popular method. There are many methods to interpret resistivity and induced polarization data. Inversion method is one of the most popular methods. This method was reported as early as the 1930s. Visual and analytical methods are used for the interpretations which are used for simple structure such as faults. These methods require a certain degree of symmetry and are suitable only for simple geological situations. Complex resistivity distributions cannot be solved by analytical methods and must numerical techniques be used. But these data alone cannot determine the location of anomalies precisely. One of the main concerns of interpreters is the selection of optimum drilling points by the results of surveys. The selection of optimum drilling point have an important impact on the reduction the cost and risk, The selection of these points can be done by integration of geoelectrical and structural data. By the progress of the computer science, many methods for processing of geophysical data were also developed. These methods include fuzzy, neural networks and genetic algorithms methods. Fuzzy method is based on following steps: 1-layer fuzzy information, 2- fuzzy inference, 3- deterministic output. Fuzzy layers of information are done by membership function, and there are different kinds of membership functions. In this study resistivity and induced polarization data have been collected using dipole-dipole array in mining index of Hamyj mine located in Birjand. Hamyj mine index that is the result of remote sensing and economic geology surveys show promising mineralization in this mine. This area comprises of volcanic rocks as andesit, dasite, volcanic breccia, metamorphic tuff. The dasit extensions in the study area are more frequent. For in detail study and modeling of the mineralization the geoelctrical study with Electrod spacing of 20m has been used. Then inversion of geoelecrtical data for a profile based on Gaussian- Newton and Newton’s methods were carried out. In this study, the models have been integrated based on knowledge-driven fuzzy logic method in a GIS environment. The best profile was chosen based on the Spriman correlation. Because the resistivity and induced polarization have negative correlation, so the profile 3 was selected as the best profile by highest negative correlation. Results of inverse modeling show that in this profile the resistivity of the left and right side of this profile are different. This difference in resistivity represents the chalcopyrite in the study area. The boundary of the difference is an attribute to a fault that is located 230m from the left side. Also induced polarization data have shown a high value of chalcopyright. In fact the fault existence in the mentioned location was approved by geological surveys. The effect of layers information of resistivity and induced polarization a decision, fuzzy functions were selected for them. Also a large membership function was used for the induced polarization. Also, small membership function was used for resistivity data. For faults in this study a 5 meters distance as buffer zone was applied. The fuzzy gamma operator was also used. This operator is a combination of multiplication and addition operator. The value of fuzzy gamma operator was selected 0.85 for the integration of layers information. The result of the fuzzy map, was converted to binary map and on this map, decision maker could easily select the optimum area. Finally, based on the results of the fuzzy modeling, the best point for drilling was proposed. Using this type of integration in order to determine the best drilling points in for porphyry copper exploration, appears to be a trusty method. Based on the fuzzy modeling for the study area the best drilling point was chosen near the fault in the study area. The results show that the optimum exploration borehole location is the point that is located at a distance 240m from the left side of the profile. Also in this study the maximum depth of drilling was calculated 50m. The applied fuzzy method in this surveyed profile with 400m length and proposed a 60m length as a promising mineralization area with high confidence. Also, in this method expert viewpoint could be applied in the procedure of integration. This method has high flexibility in integration of also information layers and could simultaneously use several information layers. This technique has simple mathematical relationship. By the results of this study this method can be used for other geophysical data.     Resistivity and induced polarization methods are used in exploration of porphyry metals for years. The resistivity method is used in the study of horizontal and vertical discontinuities in the electrical properties of the ground. The induced polarization method makes use of the capacitive action of the subsurface to locate zones where conductive minerals are disseminated within their host rocks. The simplicity of the equipment, the lower cost of the survey compared to the other methods and the abundance of interpretation methods make it a popular method. There are many methods to interpret resistivity and induced polarization data. Inversion method is one of the most popular methods. This method was reported as early as the 1930s. Visual and analytical methods are used for the interpretations which are used for simple structure such as faults. These methods require a certain degree of symmetry and are suitable only for simple geological situations. Complex resistivity distributions cannot be solved by analytical methods and must numerical techniques be used. But these data alone cannot determine the location of anomalies precisely. One of the main concerns of interpreters is the selection of optimum drilling points by the results of surveys. The selection of optimum drilling point have an important impact on the reduction the cost and risk, The selection of these points can be done by integration of geoelectrical and structural data. By the progress of the computer science, many methods for processing of geophysical data were also developed. These methods include fuzzy, neural networks and genetic algorithms methods. Fuzzy method is based on following steps: 1-layer fuzzy information, 2- fuzzy inference, 3- deterministic output. Fuzzy layers of information are done by membership function, and there are different kinds of membership functions. In this study resistivity and induced polarization data have been collected using dipole-dipole array in mining index of Hamyj mine located in Birjand. Hamyj mine index that is the result of remote sensing and economic geology surveys show promising mineralization in this mine. This area comprises of volcanic rocks as andesit, dasite, volcanic breccia, metamorphic tuff. The dasit extensions in the study area are more frequent. For in detail study and modeling of the mineralization the geoelctrical study with Electrod spacing of 20m has been used. Then inversion of geoelecrtical data for a profile based on Gaussian- Newton and Newton’s methods were carried out. In this study, the models have been integrated based on knowledge-driven fuzzy logic method in a GIS environment. The best profile was chosen based on the Spriman correlation. Because the resistivity and induced polarization have negative correlation, so the profile 3 was selected as the best profile by highest negative correlation. Results of inverse modeling show that in this profile the resistivity of the left and right side of this profile are different. This difference in resistivity represents the chalcopyrite in the study area. The boundary of the difference is an attribute to a fault that is located 230m from the left side. Also induced polarization data have shown a high value of chalcopyright. In fact the fault existence in the mentioned location was approved by geological surveys. The effect of layers information of resistivity and induced polarization a decision, fuzzy functions were selected for them. Also a large membership function was used for the induced polarization. Also, small membership function was used for resistivity data. For faults in this study a 5 meters distance as buffer zone was applied. The fuzzy gamma operator was also used. This operator is a combination of multiplication and addition operator. The value of fuzzy gamma operator was selected 0.85 for the integration of layers information. The result of the fuzzy map, was converted to binary map and on this map, decision maker could easily select the optimum area. Finally, based on the results of the fuzzy modeling, the best point for drilling was proposed. Using this type of integration in order to determine the best drilling points in for porphyry copper exploration, appears to be a trusty method. Based on the fuzzy modeling for the study area the best drilling point was chosen near the fault in the study area. The results show that the optimum exploration borehole location is the point that is located at a distance 240m from the left side of the profile. Also in this study the maximum depth of drilling was calculated 50m. The applied fuzzy method in this surveyed profile with 400m length and proposed a 60m length as a promising mineralization area with high confidence. Also, in this method expert viewpoint could be applied in the procedure of integration. This method has high flexibility in integration of also information layers and could simultaneously use several information layers. This technique has simple mathematical relationship. By the results of this study this method can be used for other geophysical data.     https://jesphys.ut.ac.ir/article_36699_3dd62fbc6670019466249a809185a876.pdf

University of Tehran Press Journal of the Earth and Space Physics 2538-371X 40 1 2014 04 21 A second order closure model for calculating the turbulence fluxes in atmospheric boundary layer A second order closure model for calculating the turbulence fluxes in atmospheric boundary layer 107 122 36700 10.22059/jesphys.2014.36700 FA Majid Farahani Assistant Professor, Space Physics Department, Institute of Geophysics, University of Tehran, Iran Ali Jalali Ph.D. Student, Western Ontario University, Waterloo, Canada Mohammad Ali Saghafi Research Assistant, Institute of Geophysics, University of Tehran, Iran Journal Article 2011 09 13 Turbulent fluxes in the atmospheric boundary layer are calculated using different models. One of the most popular and operational models in this groundwork is the second order turbulent model of Mellor and Yamada that has been widely used to simulate the planetary boundary layer (PBL). Despite its popularity, it has been shown that this model has several weeknesses that affect its accuracy. The main deficiency of this model is its inability in distinguishing between vertical and horizontal components of turbulent kinetic energy. Also, this model can not accurately predict the height of the boundary layer. This deficiency in turn can lead to errors in the parameterization of physical processes. Canuto and his colleagues attempt to overcome these deficiencies by introducing a modified model that is an improved version of Mellor and Yamada model. This later model in theory has proven to have more reliable results in computing turbulent fluxes. To verify the results of this theoretical model, the model was coded and verified using some of the available data sets. The data sets used in this study include, the Sodar data of Institute of Geophysics of University of Tehran and Mehrabad airport sounding data. Also the standard Large Edy Simulation data set (LES) introduced in GABLS (Cuxart, 2006) project was examined. For better certainty of the new model performance, the Nakanishi's LES data that used in his study was used and applied to the model. The data used here present a boundary layer of the atmosphere in neutral or Ekman layer condition. Potential temperature does not exist in the Sodar data set and therefore it was computed using measured laps rate in the same data set. The output of model includes the main flux terms, ,, . Results show that the model produces acceptable values for fluxes and have almost the same vertical profile and trend. There are some differences in the predicted values of parameters. It seems that this model could recognize eddies of boundary layer with more accuracy. The differences between the vertical and horizontal TKE distribution is clear. Some of the differences are explained. The results of model as long as the boundary layer conditions match to those that theory proposed or assumed, are reasonable. The analytical data feed to the model results in better output. Compare to a well-known intercomparison program (GABLS project) the results are quite good. The PBL height is defined as the height at which the turbulent kinetic energy or the magnitude of the momentum flux decreases to a small fraction of the corresponding surface value. Results of experiments show that boundary layer height is estimated well in this model. The critical Richardson number predicted by the model is around unity rather than 0.2 as given by previous models. The larger critical Richardson number is in agreement with LES data. Experiments revealed that the presented model is not suitable for boundary layers in which the Richardson numbers within them approach 1 or are larger than 1.     Turbulent fluxes in the atmospheric boundary layer are calculated using different models. One of the most popular and operational models in this groundwork is the second order turbulent model of Mellor and Yamada that has been widely used to simulate the planetary boundary layer (PBL). Despite its popularity, it has been shown that this model has several weeknesses that affect its accuracy. The main deficiency of this model is its inability in distinguishing between vertical and horizontal components of turbulent kinetic energy. Also, this model can not accurately predict the height of the boundary layer. This deficiency in turn can lead to errors in the parameterization of physical processes. Canuto and his colleagues attempt to overcome these deficiencies by introducing a modified model that is an improved version of Mellor and Yamada model. This later model in theory has proven to have more reliable results in computing turbulent fluxes. To verify the results of this theoretical model, the model was coded and verified using some of the available data sets. The data sets used in this study include, the Sodar data of Institute of Geophysics of University of Tehran and Mehrabad airport sounding data. Also the standard Large Edy Simulation data set (LES) introduced in GABLS (Cuxart, 2006) project was examined. For better certainty of the new model performance, the Nakanishi's LES data that used in his study was used and applied to the model. The data used here present a boundary layer of the atmosphere in neutral or Ekman layer condition. Potential temperature does not exist in the Sodar data set and therefore it was computed using measured laps rate in the same data set. The output of model includes the main flux terms, ,, . Results show that the model produces acceptable values for fluxes and have almost the same vertical profile and trend. There are some differences in the predicted values of parameters. It seems that this model could recognize eddies of boundary layer with more accuracy. The differences between the vertical and horizontal TKE distribution is clear. Some of the differences are explained. The results of model as long as the boundary layer conditions match to those that theory proposed or assumed, are reasonable. The analytical data feed to the model results in better output. Compare to a well-known intercomparison program (GABLS project) the results are quite good. The PBL height is defined as the height at which the turbulent kinetic energy or the magnitude of the momentum flux decreases to a small fraction of the corresponding surface value. Results of experiments show that boundary layer height is estimated well in this model. The critical Richardson number predicted by the model is around unity rather than 0.2 as given by previous models. The larger critical Richardson number is in agreement with LES data. Experiments revealed that the presented model is not suitable for boundary layers in which the Richardson numbers within them approach 1 or are larger than 1.     https://jesphys.ut.ac.ir/article_36700_b44724074d7242865399e596ce97728b.pdf

University of Tehran Press Journal of the Earth and Space Physics 2538-371X 40 1 2014 04 21 Wind shear and stratification effects on baroclinic adjustment Wind shear and stratification effects on baroclinic adjustment 123 136 36701 10.22059/jesphys.2014.36701 FA Ali Mohammadi Instructor of Meteorology, Imam-Khomeini Marine University, Noshahr, Iran Ali Reza Mohebalhojeh Accosiate Professor, Institute of Geophysics, University of Tehran, Iran Journal Article 2013 01 21 The baroclinic adjustment is defined as the neutralization of the mean flow by eddies. With the amplitude growth of the baroclinic waves generated by an unstable mean flow, the eddy fluxes increase and act to reduce the meridional temperature gradient or vertical wind shear of the mean flow to subcritical value. This is believed to constitute the basic mechanism of baroclinic adjustment. The concomitant reduction of the potential vorticity (PV) along isentropic surfaces during the baroclinic adjustment process can be accommodated within the quasi-geostrophic theory by removing meridional gradient of quasigeostrophic PV in the interior of the domain and satisfying the requirements of the Charny–Stern–Pedlosky theorem. The result is a group of models that extend the work of Eady and thus called Eady-like models. In addition to uniform interior PV, these models share a rigid lid as the upper boundary mimicking the tropopause.  In fact, the upper rigid helps to resolve a major issue regarding an unphysical aspect of the state obtained through  baroclinic adjustment. For the problem with a constant buoyancy frequencyfor the mean flow, as results of previous studies suggest, the major issue is that baroclinic neutrality requires that the surface value of buoyancy frequency become infinite or exceedingly large, even if frictional damping is taking into account. The mathematical tractability, physical simplicity, and yet the ability to explain some of the important features of baroclinic eddies in the real atmosphere and oceans, have made the Eady-like models of baroclinic instability the ideal setting to study baroclinic adjustment. In particular, the presence of a short-wave cut-off for instability, an immediate effect of using uniform interior PV, can provide a route to reach baroclinic neutrality. This can be achieved through limiting the interaction between the lower and upper edge waves by, for example, moving the lid to upper levels and thus displacing the short-wave cut-off to longer wavelengths. Further, confining the baroclinic eddies in meridional direction leads to a lower limit for the horizontal wavenumber.  Making the latter lower limit equal to the short-wave cut-off leads to a state of neutrality for all of the zonal modes. Within the extent of Eady-like models, this paper is devoted to study the effect of vertical wind shear and stratification on the baroclinic adjustment. Disregarding the density variation with height, analytical solutions are obtained in the form of the Bessel and Neumann functions for the instability problem with a quadratic profile for the zonal mean flow and a general form for the height variation of the squared buoyancy frequency.  The instability problem is investigated by varying  and the level of maximum zonal mean flow, that is, the location of the jet. Varying the latter parameter help us determine the sensitivity to the changes of vertical shear in the upper boundary relative to that in the lower boundary.  The relation between the eddy flux of PV and the stratification parameters including is also explored at a representative mid-tropospheric level  it is shown that the eddy flux of PV is minimum (maximum) for positive (negative) values of   which correspond to decreasing (increasing) values of with height.  Therefore, it can be concluded that the process of baroclinic adjustment is associated with negative values of Using a properly determined vertical profile for , the meridional gradient of PV is removed in the interior of the domain. It is shown that the adjusted buoyancy frequency decreases with increasing altitude. With stronger winds, the reduction with height in buoyancy frequency becomes more significant. With the meridional gradient of PV eliminated, the necessary condition for instability based on the Charny–Stern–Pedlosky criterion requires that the wind shear on the lower and upper boundaries be of the same sign. However, even in the case of equal-sign wind shear on the two boundaries, neutrality may be achieved depending on the structure of the wind shear throughout the domain and not just on the boundaries. In the Eady problem, wind shear has no effect in the condition for instability. Considering the effects of wind shear on the adjusted buoyancy frequency, it is shown that instability can occur if wind shear increases even when the surface value of is. Without considering the variation with height of, as the case with no rigid lid, neutrality can only happen for exceedingly large values of surface buoyancy frequency. It is shown that for baroclinic adjustment to occur, the two following conditions are necessary: 1- a decreasing with height of the buoyancy frequency; 2- a weak wind shear. Nevertheless, baroclinic adjustment can also occur in the presence of strong wind shear, if there is sufficiently large surface value of . Further, it is shown that the effect of changing the wind shear difference between the lower and upper boundaries is compensated by a corresponding change in the reduction of buoyancy frequency with height.     The baroclinic adjustment is defined as the neutralization of the mean flow by eddies. With the amplitude growth of the baroclinic waves generated by an unstable mean flow, the eddy fluxes increase and act to reduce the meridional temperature gradient or vertical wind shear of the mean flow to subcritical value. This is believed to constitute the basic mechanism of baroclinic adjustment. The concomitant reduction of the potential vorticity (PV) along isentropic surfaces during the baroclinic adjustment process can be accommodated within the quasi-geostrophic theory by removing meridional gradient of quasigeostrophic PV in the interior of the domain and satisfying the requirements of the Charny–Stern–Pedlosky theorem. The result is a group of models that extend the work of Eady and thus called Eady-like models. In addition to uniform interior PV, these models share a rigid lid as the upper boundary mimicking the tropopause.  In fact, the upper rigid helps to resolve a major issue regarding an unphysical aspect of the state obtained through  baroclinic adjustment. For the problem with a constant buoyancy frequencyfor the mean flow, as results of previous studies suggest, the major issue is that baroclinic neutrality requires that the surface value of buoyancy frequency become infinite or exceedingly large, even if frictional damping is taking into account. The mathematical tractability, physical simplicity, and yet the ability to explain some of the important features of baroclinic eddies in the real atmosphere and oceans, have made the Eady-like models of baroclinic instability the ideal setting to study baroclinic adjustment. In particular, the presence of a short-wave cut-off for instability, an immediate effect of using uniform interior PV, can provide a route to reach baroclinic neutrality. This can be achieved through limiting the interaction between the lower and upper edge waves by, for example, moving the lid to upper levels and thus displacing the short-wave cut-off to longer wavelengths. Further, confining the baroclinic eddies in meridional direction leads to a lower limit for the horizontal wavenumber.  Making the latter lower limit equal to the short-wave cut-off leads to a state of neutrality for all of the zonal modes. Within the extent of Eady-like models, this paper is devoted to study the effect of vertical wind shear and stratification on the baroclinic adjustment. Disregarding the density variation with height, analytical solutions are obtained in the form of the Bessel and Neumann functions for the instability problem with a quadratic profile for the zonal mean flow and a general form for the height variation of the squared buoyancy frequency.  The instability problem is investigated by varying  and the level of maximum zonal mean flow, that is, the location of the jet. Varying the latter parameter help us determine the sensitivity to the changes of vertical shear in the upper boundary relative to that in the lower boundary.  The relation between the eddy flux of PV and the stratification parameters including is also explored at a representative mid-tropospheric level  it is shown that the eddy flux of PV is minimum (maximum) for positive (negative) values of   which correspond to decreasing (increasing) values of with height.  Therefore, it can be concluded that the process of baroclinic adjustment is associated with negative values of Using a properly determined vertical profile for , the meridional gradient of PV is removed in the interior of the domain. It is shown that the adjusted buoyancy frequency decreases with increasing altitude. With stronger winds, the reduction with height in buoyancy frequency becomes more significant. With the meridional gradient of PV eliminated, the necessary condition for instability based on the Charny–Stern–Pedlosky criterion requires that the wind shear on the lower and upper boundaries be of the same sign. However, even in the case of equal-sign wind shear on the two boundaries, neutrality may be achieved depending on the structure of the wind shear throughout the domain and not just on the boundaries. In the Eady problem, wind shear has no effect in the condition for instability. Considering the effects of wind shear on the adjusted buoyancy frequency, it is shown that instability can occur if wind shear increases even when the surface value of is. Without considering the variation with height of, as the case with no rigid lid, neutrality can only happen for exceedingly large values of surface buoyancy frequency. It is shown that for baroclinic adjustment to occur, the two following conditions are necessary: 1- a decreasing with height of the buoyancy frequency; 2- a weak wind shear. Nevertheless, baroclinic adjustment can also occur in the presence of strong wind shear, if there is sufficiently large surface value of . Further, it is shown that the effect of changing the wind shear difference between the lower and upper boundaries is compensated by a corresponding change in the reduction of buoyancy frequency with height.     https://jesphys.ut.ac.ir/article_36701_9664a35ab26bcb8351d981e4decc2162.pdf

University of Tehran Press Journal of the Earth and Space Physics 2538-371X 40 1 2014 04 21 Electrical charge transfer modeling (lightning) in cloud and its implementation in a one-dimensional prognostic cloud model Electrical charge transfer modeling (lightning) in cloud and its implementation in a one-dimensional prognostic cloud model 137 148 36702 10.22059/jesphys.2014.36702 FA Maryam Gharaylou Assistant Professor, Space Physics Department, Institute of Geophysics, University of Tehran, Iran Nafiseh Pegahfar Assistant Professor, Iranian National Institute for Oceanography and Atmospheric Science, Tehran, Iran Abbas Ali Bidokhti, Professor, Space Physics Department, Institute of Geophysics, University of Tehran, Iran Journal Article 2013 05 04 One of the most dazzling events in the atmosphere is lightning. During updrafts in the life cycle of cumulus clouds, collision of graupels and ice crystals in the presence of liquid water results in vertical separation of electrical charges and lightning. There are four types of lightening depending on the location of discharge. The first type is cloud-to-ground lightning or fork lightning that happens due to electrical discharge between the cloud base and the negatively charged earth. Where the regions with opposite electrical charge within a cloud are connected, intra-cloud lightning occurs. The third one is lightning between clouds with opposite electric charge namely cloud-to-cloud lightening or sheet or heat lightening and the last one is known as cloud-to-air lightening. In spite of the fact that lightning is considered as a part of severe weather systems, but it is hard to be predicted in short-term prediction. A thunderstorm could contain several tens of Coulombs of charge. The negative charge region has a temperature of -5 to -10⁰C while the positive charge region is located 2-3km upper than the negative charge location. In the gravity separation theory (principle of this research) some microphysical processes lead to charge separation. Negative charges are carried by heavier particles (cloud droplets, ice crystals and ions). Therefore during precipitation, negative charges accumulate in the lower levels while positive charges moved upwards by updrafts within the cloud. In this paper, it is shown that the charge transfer due to interaction of charged particles (collision of graupels and ice crystals) is the most important ionization process. Two processes of non-inductive and inductive ionization are almost considered. In non-inductive ionization process, collision of hydrometeors results in charge separation.  Whilst in the inductive one the existence of an external field induces polarization, and then charge separation occurs. However, non-inductive process that happens because of collision between the graupels and ice crystals in the presence of liquid water has the most significant role in the charge transfer (Sanders et al. 1991, Miller et al. 2001). In this work, the non-inductive ionization mechanism is applied.   In this research outputs of a one-dimensional cloud model were used (including vertical velocity, mixing ratios of graupels and ice, liquid water content, terminal velocity of graupels and temperature) to simulate the charge transfer intracloud (at microphysics scale). The vertical one-dimensional cloud model (Explicit Time-dependent Model (ETM)) is based on Chen and Sun (2002) equations (Gharaylou et al., 2009). The cloud in this model is considered as a cylindrical column of air with a constant radius. Non-hydrostatic pressure is assumed within the cloud column while the environment is in hydrostatic equilibrium. In the cloud model, the microphysical processes such as evaporation/sublimation, deposition/condensation, melting, freezing, aggregation, accretion and Bergeron process, entrainment and detrainment and lateral and vertical eddy mixing effects have been considered. In the ETM model convection is initiated using a potential temperature perturbation based on Chen and Sun (2004) relation. To proceed, an idealized sounding was used as the input data. This profile consists of temperature, relative humidity and ambient pressure. Surface temperature and relative humidity are equal to 298 K and 94.5%. Its temperature profile was determined according to dry adiabatic lapse rate below 1 km, saturated adiabatic lapse rate from 1 to 10 km and isothermal for upper levels. The relative humidity increased linearly below 1 km and afterwards decreased with a rate of 5%. The vertical velocity initialized based on Ogura and Takahashi (1971) relations. In this research, the mean charge transferred per collision of graupels and ice crystals was simulated using parametric equations suggested by Sanders et al (1991). The simulation has been done for 70 minutes with 1 second time step. The vertical resolution was typically 250 meter up to 15 km above ground level. The studied parameters consist of relative vertical velocity of graupels, mixing ratios of graupels and ice, electric filed and mean charge transfer per collision. The results showed that lightning happened between 35-50 minutes inward simulation. It is worth to note that simultaneous presence of graupels and ice crystals guarantees the initiation of electric field in the cloud. Once the electric field intensity (positive or negative) exceeded the threshold electric field defined by Marshall et al. (1995), lightning occurred. Decrease of mixing ratios of graupels and ice crystals leads to weakened electric field in the upper levels. Precipitation also results in electric fields form at the lower altitudes. The negative sign of electric field can be inferred from the negative charge transfer in collisions.     One of the most dazzling events in the atmosphere is lightning. During updrafts in the life cycle of cumulus clouds, collision of graupels and ice crystals in the presence of liquid water results in vertical separation of electrical charges and lightning. There are four types of lightening depending on the location of discharge. The first type is cloud-to-ground lightning or fork lightning that happens due to electrical discharge between the cloud base and the negatively charged earth. Where the regions with opposite electrical charge within a cloud are connected, intra-cloud lightning occurs. The third one is lightning between clouds with opposite electric charge namely cloud-to-cloud lightening or sheet or heat lightening and the last one is known as cloud-to-air lightening. In spite of the fact that lightning is considered as a part of severe weather systems, but it is hard to be predicted in short-term prediction. A thunderstorm could contain several tens of Coulombs of charge. The negative charge region has a temperature of -5 to -10⁰C while the positive charge region is located 2-3km upper than the negative charge location. In the gravity separation theory (principle of this research) some microphysical processes lead to charge separation. Negative charges are carried by heavier particles (cloud droplets, ice crystals and ions). Therefore during precipitation, negative charges accumulate in the lower levels while positive charges moved upwards by updrafts within the cloud. In this paper, it is shown that the charge transfer due to interaction of charged particles (collision of graupels and ice crystals) is the most important ionization process. Two processes of non-inductive and inductive ionization are almost considered. In non-inductive ionization process, collision of hydrometeors results in charge separation.  Whilst in the inductive one the existence of an external field induces polarization, and then charge separation occurs. However, non-inductive process that happens because of collision between the graupels and ice crystals in the presence of liquid water has the most significant role in the charge transfer (Sanders et al. 1991, Miller et al. 2001). In this work, the non-inductive ionization mechanism is applied.   In this research outputs of a one-dimensional cloud model were used (including vertical velocity, mixing ratios of graupels and ice, liquid water content, terminal velocity of graupels and temperature) to simulate the charge transfer intracloud (at microphysics scale). The vertical one-dimensional cloud model (Explicit Time-dependent Model (ETM)) is based on Chen and Sun (2002) equations (Gharaylou et al., 2009). The cloud in this model is considered as a cylindrical column of air with a constant radius. Non-hydrostatic pressure is assumed within the cloud column while the environment is in hydrostatic equilibrium. In the cloud model, the microphysical processes such as evaporation/sublimation, deposition/condensation, melting, freezing, aggregation, accretion and Bergeron process, entrainment and detrainment and lateral and vertical eddy mixing effects have been considered. In the ETM model convection is initiated using a potential temperature perturbation based on Chen and Sun (2004) relation. To proceed, an idealized sounding was used as the input data. This profile consists of temperature, relative humidity and ambient pressure. Surface temperature and relative humidity are equal to 298 K and 94.5%. Its temperature profile was determined according to dry adiabatic lapse rate below 1 km, saturated adiabatic lapse rate from 1 to 10 km and isothermal for upper levels. The relative humidity increased linearly below 1 km and afterwards decreased with a rate of 5%. The vertical velocity initialized based on Ogura and Takahashi (1971) relations. In this research, the mean charge transferred per collision of graupels and ice crystals was simulated using parametric equations suggested by Sanders et al (1991). The simulation has been done for 70 minutes with 1 second time step. The vertical resolution was typically 250 meter up to 15 km above ground level. The studied parameters consist of relative vertical velocity of graupels, mixing ratios of graupels and ice, electric filed and mean charge transfer per collision. The results showed that lightning happened between 35-50 minutes inward simulation. It is worth to note that simultaneous presence of graupels and ice crystals guarantees the initiation of electric field in the cloud. Once the electric field intensity (positive or negative) exceeded the threshold electric field defined by Marshall et al. (1995), lightning occurred. Decrease of mixing ratios of graupels and ice crystals leads to weakened electric field in the upper levels. Precipitation also results in electric fields form at the lower altitudes. The negative sign of electric field can be inferred from the negative charge transfer in collisions.     https://jesphys.ut.ac.ir/article_36702_f063d53aa3a5c845e1b246dcdd0772fa.pdf

University of Tehran Press Journal of the Earth and Space Physics 2538-371X 40 1 2014 04 21 Simulation of inertia-gravity waves using WRF model over Iran: a case study Simulation of inertia-gravity waves using WRF model over Iran: a case study 149 166 36703 10.22059/jesphys.2014.36703 FA Mojgan AmirAmjadi M. Sc. Graduated Student of Meteorology, Space Physics Department, Institute of Geophysics, University of Tehran, Iran Alireza Mohebalhojeh Associate Professor, Space Physics Department, Institute of Geophysics, University of Tehran, Iran Mohammad Mirzaei Assistant Professor, Space Physics Department, Institute of Geophysics, University of Tehran, Iran Journal Article 2013 05 18 A fluid, which is stable under the action of buoyancy, can oscillate under the influence of buoyancy and Coriolis forces. The resulting ageostrophic oscillations with a frequency between Coriolis and buoyancy frequencies are called inertia-gravity waves, hereafter IGWs. In this study we investigate the generation and propagation of IGWs over Iran. To this end, the Weather Research and Forecasting (WRF) mesoscale model is used to simulate an event with noticeable IGWs within a synoptic system accompanied by rain and snow, from 6 to 9 February, 2012. The NCEP FNL data (final analyses) are used for this simulation in a domain of 8000 km × 6375 km extent and a medium resolution (grid interval Δs=25 km). The model has 35 levels in vertical direction with its top at 10 hPa (~30 km). A time step of 150 seconds is used and the model is run for 72 hours initialized at 12UTC 6 February and continued until 12UTC 9 February 2012. An implicit Rayleigh damping layer is used to prevent unphysical wave reflection from the upper boundary of the computational domain (Klemp et al, 2008). In order to investigate generation and propagation of IGWs and to determine the effects of each of energy sources in the region, four different model runs are designed and performed with/without orography/moisture. In order to facilitate the conduct of the study, the main area of interest is classified into three regions where the IGWs are most active. The first area is located in the northwest of the country, in the vicinity of the maximum wind speed of the midlatitude jet stream and the left side of the jet streak. The second area, including the central and southern Zagros Mountain, is located in the vicinity of unbalanced flow and the right side of the jet streak. The third area is located geographically at the conjunction of Zagros and Alborz mountain ranges, is at the exit region of the jet stream and along a rain belt with significant cloud coverage as in the second area. The fundamental quantities of wave like the wave frequency and periods, the intrinsic phase speed, the group velocity and horizontal and vertical wavelengths are obtained based on the horizontal divergence field as the main quantity. This is possible, because the procedure avoids explicit treatment of the background field, which has zero divergence, and is applicable to waves of arbitrary wavelength. Previous studies have shown that the most important energy sources for IGWs are: jet streams, fronts, convection and orography. Furthermore, the IGWs propagate in the atmosphere with a phase speed of 15 to 35 ms^-1, vertical wavelength of 500 m to 15 km and horizontal wavelength of 50 to 1000 km. In the reference run (with the topography and moisture included), the distribution of the horizontal divergence clearly shows the waves closely follow the major topography and propagate nearly perpendicular to it. Estimation of wave properties shows that a high- frequency wave with  is emitted in this case. The quantity  is the estimate for wave frequency scaled by inertial frequency. In the first and second areas, the typical values of the horizontal wavelengths are from 150 to 175 km, the vertical wavelength from 5 to 6 km and intrinsic phase speed from 15 to 19 ms^-1. In the third area, the IGWs travel with a phase speed of about 15 to 21 ms^-1 and horizontal and vertical wavelengths of 100 to 120 km and 5 to 6 km, respectively. Based on the characteristics of the IGWs and their propagation, it can be inferred that the waves are generated in the troposphere by jet-front mechanisms due to topography and then undergo deformation. Some of the waves generated in the upper troposphere cross the tropopause and propagate well into the stratosphere. Disregard of the way wave are generated, this transfer of activity from troposphere to the stratosphere is a common phenomena.  The results of this study are in good agreement with those obtained by Zhang and Koch in 2000 who studied simulation over Rockies Mountain in Montana and Dakota. They estimated a single wave packet with 3 or 4 waves, a wavelength of 150 km and phase speed of 15.2 ms^-1.  Zulick and Peters in 2006 identified two wave packets in their study. The first wave packet included large waves with wavelengths of about 500 km and period of approximately 10 hour which were classified as sub-synoptic (or meso-α). The second wave packet, however, consisted of small waves with much higher frequency than the first wave packet, and wavelengths of almost 200 km and period of nearly 5 hour, which were classified as meso-β (mesoscale-β). Having investigated the possible energy sources in the region, the following conclusions can be made. In the run without either orography or moisture, large-amplitude waves are observed that carry energy upward to the stratosphere and downward to the troposphere. Entering a positive feedback, the downward propagating waves emitted by the jet can facilitate convective activity at lower levels and play a role in enhancing precipitation. Although adding mountains to the physics of the model affects wave characteristics like intensity, amplitude, and direction of propagation, the presence of all of the factors is necessary for describing the actual wave generation and propagation processes.       A fluid, which is stable under the action of buoyancy, can oscillate under the influence of buoyancy and Coriolis forces. The resulting ageostrophic oscillations with a frequency between Coriolis and buoyancy frequencies are called inertia-gravity waves, hereafter IGWs. In this study we investigate the generation and propagation of IGWs over Iran. To this end, the Weather Research and Forecasting (WRF) mesoscale model is used to simulate an event with noticeable IGWs within a synoptic system accompanied by rain and snow, from 6 to 9 February, 2012. The NCEP FNL data (final analyses) are used for this simulation in a domain of 8000 km × 6375 km extent and a medium resolution (grid interval Δs=25 km). The model has 35 levels in vertical direction with its top at 10 hPa (~30 km). A time step of 150 seconds is used and the model is run for 72 hours initialized at 12UTC 6 February and continued until 12UTC 9 February 2012. An implicit Rayleigh damping layer is used to prevent unphysical wave reflection from the upper boundary of the computational domain (Klemp et al, 2008). In order to investigate generation and propagation of IGWs and to determine the effects of each of energy sources in the region, four different model runs are designed and performed with/without orography/moisture. In order to facilitate the conduct of the study, the main area of interest is classified into three regions where the IGWs are most active. The first area is located in the northwest of the country, in the vicinity of the maximum wind speed of the midlatitude jet stream and the left side of the jet streak. The second area, including the central and southern Zagros Mountain, is located in the vicinity of unbalanced flow and the right side of the jet streak. The third area is located geographically at the conjunction of Zagros and Alborz mountain ranges, is at the exit region of the jet stream and along a rain belt with significant cloud coverage as in the second area. The fundamental quantities of wave like the wave frequency and periods, the intrinsic phase speed, the group velocity and horizontal and vertical wavelengths are obtained based on the horizontal divergence field as the main quantity. This is possible, because the procedure avoids explicit treatment of the background field, which has zero divergence, and is applicable to waves of arbitrary wavelength. Previous studies have shown that the most important energy sources for IGWs are: jet streams, fronts, convection and orography. Furthermore, the IGWs propagate in the atmosphere with a phase speed of 15 to 35 ms^-1, vertical wavelength of 500 m to 15 km and horizontal wavelength of 50 to 1000 km. In the reference run (with the topography and moisture included), the distribution of the horizontal divergence clearly shows the waves closely follow the major topography and propagate nearly perpendicular to it. Estimation of wave properties shows that a high- frequency wave with  is emitted in this case. The quantity  is the estimate for wave frequency scaled by inertial frequency. In the first and second areas, the typical values of the horizontal wavelengths are from 150 to 175 km, the vertical wavelength from 5 to 6 km and intrinsic phase speed from 15 to 19 ms^-1. In the third area, the IGWs travel with a phase speed of about 15 to 21 ms^-1 and horizontal and vertical wavelengths of 100 to 120 km and 5 to 6 km, respectively. Based on the characteristics of the IGWs and their propagation, it can be inferred that the waves are generated in the troposphere by jet-front mechanisms due to topography and then undergo deformation. Some of the waves generated in the upper troposphere cross the tropopause and propagate well into the stratosphere. Disregard of the way wave are generated, this transfer of activity from troposphere to the stratosphere is a common phenomena.  The results of this study are in good agreement with those obtained by Zhang and Koch in 2000 who studied simulation over Rockies Mountain in Montana and Dakota. They estimated a single wave packet with 3 or 4 waves, a wavelength of 150 km and phase speed of 15.2 ms^-1.  Zulick and Peters in 2006 identified two wave packets in their study. The first wave packet included large waves with wavelengths of about 500 km and period of approximately 10 hour which were classified as sub-synoptic (or meso-α). The second wave packet, however, consisted of small waves with much higher frequency than the first wave packet, and wavelengths of almost 200 km and period of nearly 5 hour, which were classified as meso-β (mesoscale-β). Having investigated the possible energy sources in the region, the following conclusions can be made. In the run without either orography or moisture, large-amplitude waves are observed that carry energy upward to the stratosphere and downward to the troposphere. Entering a positive feedback, the downward propagating waves emitted by the jet can facilitate convective activity at lower levels and play a role in enhancing precipitation. Although adding mountains to the physics of the model affects wave characteristics like intensity, amplitude, and direction of propagation, the presence of all of the factors is necessary for describing the actual wave generation and propagation processes.       https://jesphys.ut.ac.ir/article_36703_5f65779385b869a20f71cc0bc77311cd.pdf

University of Tehran Press Journal of the Earth and Space Physics 2538-371X 40 1 2014 04 21 Detection of climate change effect on meteorological droughts in northwest of Iran Detection of climate change effect on meteorological droughts in northwest of Iran 167 184 36704 10.22059/jesphys.2014.36704 FA Mahdi Ghamghami Ph.D. Student of Agro-Meteorology, University of Tehran, Iran Nozar Ghahreman Associate Professor, Department of Irrigation and Reclamation Engineering, University of Tehran, Iran Somayeh Hejabi Ph.D. Student of Agro-Meteorology, University of Tehran, Iran Journal Article 2013 07 16 Climate change that the human faces is a somewhat unavoidable phenomenon. Successful management of water resources needs recognition and perception of climate change in order to cope with water scarcity. The water scarcity is created by natural forcings such as drought, which is affected by regional climate. In other words, variation of climate variables as a result of climate change leads to variations in drought severity and frequency. Since climate change scenarios are based on assumption of increasing, decreasing or non-significant trend in climatic means, It is expected that the effects of these assumptions would be reflected in the prediction of meteorological phenomena like  drought. In Markov analysis, these variations are determined as change in transfer probability function values or shift in drought severity class, which are both important in management decisions. For instance, by increasing the temperature or decreasing the rainfall it is expected that occurrence of a drought event under certain conditions would be more probable. In this study, the outputs of three General Circulation Models (GCMs) namely; ECHO-G, CGCM3T63 HADCM3 under three climate change scenarios were downscaled using a non-parametric approach for simulation of rainfall and temperature series during 2011-2040 in northwest of Iran. This downscaling approach is combination of two techniques i.e. Kernel probability density function estimator (KDE) and Strategic Re-sampling method by which predicted variations of GCM outputs are extended and transformed to generated time series of a given future period. In KDE method, A probability density function is defined with center value of i^th observation from series (<em>x_i , i=1,...,n</em>). Contribution of each observation in estimation of probability density function of <em>i</em>_th observation is estimated by this Kernel function. The main parameter of this function is the bandwidth which is, by mathematical definition, a distance on x-axis in which the function variation is insignificant. Firstly a random normal kernel is selected and its average is considered as the base vector. Selection probability of each vector is 1/n. Then by calculation of cumulative probability and comparison with a random number between 0 and 1, one of the normal kernels is selected for rest of the simulations. The strategic re-sampling method uses a rule for generating series with specific feature such as increasing frequency of warmer or more rainy days. The criteria for such features are selected by the user based on the outputs of GCMs. Considering its semi random nature, this approach cannot be used alone for regional climate change simulations and should be combined with a weather generator such that the applied rule should be run on observed or historical series. Then, the outputs are feed in weather generator for generating a completely random series coincide with climatic scenario. After simulation of climate, Reconnaissance Drought Index (RDI) was used for monitoring drought during two periods 1971-2000 and 2011 to 2040 in northwest of Iran. This index uses the ratio of precipitation and evapotranspiration (calculated by Thorntwait method), hence as the index becomes smaller, more severe would be the drought. Thus, the necessary variables for RDI estimation are monthly mean temperature and total rainfall. For RDI calculation, firstly, the precipitation (prec) and potential evapotranspiration (PET) are calculated cumulatively with determination of the moving window value, and then, RDI values are obtained as logarithm of cumulative prec to PET ratio. Four classes are considered for RDI including: normal class (larger than -1), moderately drought class (-1 to -1.5), severe drought class (-1.5 to -2) and very severe drought (less than -2). Taking into account the length of the dry periods in the arid regions of the country, the Reconnaissance Drought Index in 6-month timescale was used for drought monitoring.  Markov analysis was applied for calculation of transfer probability and corresponding drought severity classes with three steps forward to assess the sequences of climatic assumptions on management early warnings. Behavior of a Markov model is determined by a series of probabilities in transition from one state to another namely transition probabilities. These probabilities may vary by climate change. The first-order Markov chain model was employed to predict drought condition up to 3-step ahead. This model was fitted on the RDI series at all stations of interest, and it was identified that it can represent the probabilistic behavior of drought over northwest of Iran. Research findings are presented in three parts of downscaling method implementation, RDI monitoring and Markov analysis. The weather generator model was successful for simulation of monthly normals including means and standard deviation. Also, strategic re-sampling technique as aligning method was successful for simulation of deviations from normal. Drought monitoring with RDI showed a water tension resonance in second 15 years of 1971-2000 periods. Likewise, in part of Markov analysis, findings of this study revealed that under conditions of increasing temperature and decreasing rainfall, as the worst case, the effect of climate change on meteorological drought would appear as the class shift, and in most of the study stations under this scenario, increased duration of extremely drought (class 4) was forecasted, even 2 steps ahead, which is important in water resource management.     Climate change that the human faces is a somewhat unavoidable phenomenon. Successful management of water resources needs recognition and perception of climate change in order to cope with water scarcity. The water scarcity is created by natural forcings such as drought, which is affected by regional climate. In other words, variation of climate variables as a result of climate change leads to variations in drought severity and frequency. Since climate change scenarios are based on assumption of increasing, decreasing or non-significant trend in climatic means, It is expected that the effects of these assumptions would be reflected in the prediction of meteorological phenomena like  drought. In Markov analysis, these variations are determined as change in transfer probability function values or shift in drought severity class, which are both important in management decisions. For instance, by increasing the temperature or decreasing the rainfall it is expected that occurrence of a drought event under certain conditions would be more probable. In this study, the outputs of three General Circulation Models (GCMs) namely; ECHO-G, CGCM3T63 HADCM3 under three climate change scenarios were downscaled using a non-parametric approach for simulation of rainfall and temperature series during 2011-2040 in northwest of Iran. This downscaling approach is combination of two techniques i.e. Kernel probability density function estimator (KDE) and Strategic Re-sampling method by which predicted variations of GCM outputs are extended and transformed to generated time series of a given future period. In KDE method, A probability density function is defined with center value of i^th observation from series (<em>x_i , i=1,...,n</em>). Contribution of each observation in estimation of probability density function of <em>i</em>_th observation is estimated by this Kernel function. The main parameter of this function is the bandwidth which is, by mathematical definition, a distance on x-axis in which the function variation is insignificant. Firstly a random normal kernel is selected and its average is considered as the base vector. Selection probability of each vector is 1/n. Then by calculation of cumulative probability and comparison with a random number between 0 and 1, one of the normal kernels is selected for rest of the simulations. The strategic re-sampling method uses a rule for generating series with specific feature such as increasing frequency of warmer or more rainy days. The criteria for such features are selected by the user based on the outputs of GCMs. Considering its semi random nature, this approach cannot be used alone for regional climate change simulations and should be combined with a weather generator such that the applied rule should be run on observed or historical series. Then, the outputs are feed in weather generator for generating a completely random series coincide with climatic scenario. After simulation of climate, Reconnaissance Drought Index (RDI) was used for monitoring drought during two periods 1971-2000 and 2011 to 2040 in northwest of Iran. This index uses the ratio of precipitation and evapotranspiration (calculated by Thorntwait method), hence as the index becomes smaller, more severe would be the drought. Thus, the necessary variables for RDI estimation are monthly mean temperature and total rainfall. For RDI calculation, firstly, the precipitation (prec) and potential evapotranspiration (PET) are calculated cumulatively with determination of the moving window value, and then, RDI values are obtained as logarithm of cumulative prec to PET ratio. Four classes are considered for RDI including: normal class (larger than -1), moderately drought class (-1 to -1.5), severe drought class (-1.5 to -2) and very severe drought (less than -2). Taking into account the length of the dry periods in the arid regions of the country, the Reconnaissance Drought Index in 6-month timescale was used for drought monitoring.  Markov analysis was applied for calculation of transfer probability and corresponding drought severity classes with three steps forward to assess the sequences of climatic assumptions on management early warnings. Behavior of a Markov model is determined by a series of probabilities in transition from one state to another namely transition probabilities. These probabilities may vary by climate change. The first-order Markov chain model was employed to predict drought condition up to 3-step ahead. This model was fitted on the RDI series at all stations of interest, and it was identified that it can represent the probabilistic behavior of drought over northwest of Iran. Research findings are presented in three parts of downscaling method implementation, RDI monitoring and Markov analysis. The weather generator model was successful for simulation of monthly normals including means and standard deviation. Also, strategic re-sampling technique as aligning method was successful for simulation of deviations from normal. Drought monitoring with RDI showed a water tension resonance in second 15 years of 1971-2000 periods. Likewise, in part of Markov analysis, findings of this study revealed that under conditions of increasing temperature and decreasing rainfall, as the worst case, the effect of climate change on meteorological drought would appear as the class shift, and in most of the study stations under this scenario, increased duration of extremely drought (class 4) was forecasted, even 2 steps ahead, which is important in water resource management.     https://jesphys.ut.ac.ir/article_36704_9479fd8a0bd1e70eab6e30105a0597d9.pdf

University of Tehran Press Journal of the Earth and Space Physics 2538-371X 40 1 2014 04 21 The tests of IN seeding materials used in Iran for cloud seeding (PV-26) and introduction of a cloud chamber The tests of IN seeding materials used in Iran for cloud seeding (PV-26) and introduction of a cloud chamber 185 200 36705 10.22059/jesphys.2014.36705 FA Kianoush Keshavarzian Ph.D. in Physics, Research Staff of National Research Center for Cloud Seeding, Yazd, Iran Journal Article 2013 06 08 The aim of this research is to report the results of testing of the PV-26 pyro-techniques materials and to report the operating conditions of a cloud and an aerosol chambers and also to introduce and describe these chambers of National Research Center of Cloud Seeding in Iran for further researches. The effect of ice nuclei produced by PV-26 pyro-techniques on artificial clouds has been surveyed under laboratory conditions. These tests show that this materials cause of increase the terminal velocities of ice crystals by  using an exponential fit, and by  using a linear fit. Also having captured using a microscope, the vivid photographs of ice crystals which, in the temperature  and  <br />the relative humidity less than 95%, had grown to the size of about , prove the reliability of the cloud chamber for scientific research. More than 54% of crystals have been greater than 20 microns and fit well in Magono-Lee classification. <br />Numerous tests have been performed to examine the toleration of chambers (especially, the more complicated one, cloud chamber) in extreme conditions and also to examine the capability of cloud chamber to produce reliable data. The results of these tests show that the gradient of the temperature in cloud chamber which is equipped with four temperature, one humidity and one pressure sensor, has a reasonable span. These testes are still being carried out. <br />In case of cloud chamber, the temperature could be varied from temperature of the environment to, the pressure, from pressure of environment to  below that, and the humidity, from humidity of environment to the maximum value that environmental conditions allow. Both of the chambers have centrifugal pump and filters to clean their internal air from background aerosols. The cloud chamber has an extra pump to keep the inside pressure at the desired value (Langsdorf, 1939; Donnan & wright, 1969). This chamber also has four thermometers to control the chiller compressors and monitor the gradient of the temperature. <br />At the present, the laboratory is equipped with two separate measurement tools. One of them consist a (Helium-Neon) laser and a digital detector for visible light which measures the intensity of the laser beam (Wagner et al., 2009). This system works like a turbidimeter and measures the transparency inside the cloud chamber. Having had a greater terminal velocity, bigger ice crystals precipitate faster and leave the internal space of cloud chamber for laser beam to pass (Heymsfield 1972, Heymsfield & Kajikawa 1987). When all other quantities are constant, the time interval in which the clarity of the cloud returns to the value corresponding to empty chamber, could be a criterion of the size of ice crystals (Wagner et al., 2006). <br />The other measurement system consists of an optical microscope and a special chemical gel to inscribe the shape of the crystals. The number of crystals and their shapes could be varied with respect to the physical circumstances like temperature, pressure, the number and type of the aerosols, <em>etc</em> (Mitchell et al., 1990). Although we know the ingredients of this gel, its suitable usage seems to be hard to achieve and there are still difficulties to make it practical. I will give more details in section 5. According to the Magono-Lee classification of natural snow crystals (Magono and Chung 1966), in these experiments, <em>P1b</em>, <em>P1c</em>, <em>C1h</em> and <em>CP1c</em> made the major portion amongst other crystals. <br />  <br />  <br />  The aim of this research is to report the results of testing of the PV-26 pyro-techniques materials and to report the operating conditions of a cloud and an aerosol chambers and also to introduce and describe these chambers of National Research Center of Cloud Seeding in Iran for further researches. The effect of ice nuclei produced by PV-26 pyro-techniques on artificial clouds has been surveyed under laboratory conditions. These tests show that this materials cause of increase the terminal velocities of ice crystals by  using an exponential fit, and by  using a linear fit. Also having captured using a microscope, the vivid photographs of ice crystals which, in the temperature  and  <br />the relative humidity less than 95%, had grown to the size of about , prove the reliability of the cloud chamber for scientific research. More than 54% of crystals have been greater than 20 microns and fit well in Magono-Lee classification. <br />Numerous tests have been performed to examine the toleration of chambers (especially, the more complicated one, cloud chamber) in extreme conditions and also to examine the capability of cloud chamber to produce reliable data. The results of these tests show that the gradient of the temperature in cloud chamber which is equipped with four temperature, one humidity and one pressure sensor, has a reasonable span. These testes are still being carried out. <br />In case of cloud chamber, the temperature could be varied from temperature of the environment to, the pressure, from pressure of environment to  below that, and the humidity, from humidity of environment to the maximum value that environmental conditions allow. Both of the chambers have centrifugal pump and filters to clean their internal air from background aerosols. The cloud chamber has an extra pump to keep the inside pressure at the desired value (Langsdorf, 1939; Donnan & wright, 1969). This chamber also has four thermometers to control the chiller compressors and monitor the gradient of the temperature. <br />At the present, the laboratory is equipped with two separate measurement tools. One of them consist a (Helium-Neon) laser and a digital detector for visible light which measures the intensity of the laser beam (Wagner et al., 2009). This system works like a turbidimeter and measures the transparency inside the cloud chamber. Having had a greater terminal velocity, bigger ice crystals precipitate faster and leave the internal space of cloud chamber for laser beam to pass (Heymsfield 1972, Heymsfield & Kajikawa 1987). When all other quantities are constant, the time interval in which the clarity of the cloud returns to the value corresponding to empty chamber, could be a criterion of the size of ice crystals (Wagner et al., 2006). <br />The other measurement system consists of an optical microscope and a special chemical gel to inscribe the shape of the crystals. The number of crystals and their shapes could be varied with respect to the physical circumstances like temperature, pressure, the number and type of the aerosols, <em>etc</em> (Mitchell et al., 1990). Although we know the ingredients of this gel, its suitable usage seems to be hard to achieve and there are still difficulties to make it practical. I will give more details in section 5. According to the Magono-Lee classification of natural snow crystals (Magono and Chung 1966), in these experiments, <em>P1b</em>, <em>P1c</em>, <em>C1h</em> and <em>CP1c</em> made the major portion amongst other crystals. <br />  <br />  <br />  https://jesphys.ut.ac.ir/article_36705_daf9cfcd0f0bb2f78d944a219243f228.pdf

University of Tehran Press Journal of the Earth and Space Physics 2538-371X 40 1 2014 04 21 A comparison of dust emission schemes in estimation of vertical dust flux in dust source regions of Iran A comparison of dust emission schemes in estimation of vertical dust flux in dust source regions of Iran 201 216 36706 10.22059/jesphys.2014.36706 FA Lida Khosrosereshki M.Sc.in Meteorology, Space Physics Department, Institute of Geophysics, University of Tehran, Iran Parviz Irannejad Associate Professor, Space Physics Department, Institute of Geophysics, University of Tehran, Iran Abbas Ali Aliakbari-Bidokhti Professor, Space Physics Department, Institute of Geophysics, University of Tehran, Iran Journal Article 2013 09 10 The first step in studying dust events is the estimation of dust emission from sources. In this study, the simulations of the vertical dust flux (<em>VDF</em>) over Iran by the MBA (Marticorena, Bergametti, and Alfaro) and Shao schemes are compared. The MBA scheme is energy-based and derived on the basis of energy balance between the kinetic energy of saltating particles and binding energy of the surface during the bombardment. The Shao scheme, on the other hand, is a volume removal- based approach, in which dust emission rate is estimated on the basis of the volume removed by saltating particles when they impact the surface. Because of the difficulty in measuring interparticle binding forces (or energy), Shao scheme uses the plastic pressure. This quantity is directly related to interparticle binding forces (or energy). The plastic pressure, that shows the surface resistance against the penetration of impacting particles, varies between Pa for light spray fine soil and Pa for deep wetted soil. <br />Another difference between the Shao and MBA schemes is in their parameterizations of the threshold friction velocity, at which the wind erosion is initiated. The surface and soil- related factors, which are introduced as roughness and moisture correction functions, are different in Shao and MBA schemes. Final vertical dust flux estimation is different among various versions of the Shao scheme. In this study, we use the version introduced by Shao (2004). The required data to estimate vertical dust flux by the dust emission schemes, including friction velocity, soil moisture, soil texture, surface cover fraction and roughness length were obtained from WRF model. Soil particle size distribution and plastic pressure for different grid squares are estimated based the soil texture data by reviewing the literature.        <br />In this study, dust source regions and emission rates estimated by different parameterization schemes are compared for the 24^th May 2012, when dust event was reported in many stations of Iran. According to results, compared to the Shao scheme, the MBA predicted smaller number of source regions in Iran. The vertical dust flux simulated by MBA is higher than that by the Shao scheme in some common source regions, but lower in some others. After applying roughness and moisture correction functions to the threshold friction velocity, the <em>VDF</em> simulated by the MBA will be zero wherever the Shao schemes’ <em>VDF</em> is zero, but this is not necessarily true the other way. Both of the schemes simulated the maximum <em>VDF </em>at 11am. The highest value of the estimated <em>VDF</em> was  (for the grid centered at , ) by the Shao scheme, while it was   (or the grid centered at , ) for the MBA scheme. <br />The most influential parameters in the estimation of <em>VDF</em> are the roughness length and plastic pressure in MBA and Shao, respectively. Impacts of the friction velocity, and surface and soil- related factors on the estimated vertical dust flux, though complicated, should be considered simultaneously. In every source region, <em>VDF </em>increases when the soil moisture decreases and the friction velocity increases. Furthermore, the <em>VDF</em> variation during a day is mainly affected by daily changes of the friction velocity, because of negligible changes in the soil moisture during the day. Thorough evaluation of these schemes requires experimental and remote sensing data. <br />  <br />  The first step in studying dust events is the estimation of dust emission from sources. In this study, the simulations of the vertical dust flux (<em>VDF</em>) over Iran by the MBA (Marticorena, Bergametti, and Alfaro) and Shao schemes are compared. The MBA scheme is energy-based and derived on the basis of energy balance between the kinetic energy of saltating particles and binding energy of the surface during the bombardment. The Shao scheme, on the other hand, is a volume removal- based approach, in which dust emission rate is estimated on the basis of the volume removed by saltating particles when they impact the surface. Because of the difficulty in measuring interparticle binding forces (or energy), Shao scheme uses the plastic pressure. This quantity is directly related to interparticle binding forces (or energy). The plastic pressure, that shows the surface resistance against the penetration of impacting particles, varies between Pa for light spray fine soil and Pa for deep wetted soil. <br />Another difference between the Shao and MBA schemes is in their parameterizations of the threshold friction velocity, at which the wind erosion is initiated. The surface and soil- related factors, which are introduced as roughness and moisture correction functions, are different in Shao and MBA schemes. Final vertical dust flux estimation is different among various versions of the Shao scheme. In this study, we use the version introduced by Shao (2004). The required data to estimate vertical dust flux by the dust emission schemes, including friction velocity, soil moisture, soil texture, surface cover fraction and roughness length were obtained from WRF model. Soil particle size distribution and plastic pressure for different grid squares are estimated based the soil texture data by reviewing the literature.        <br />In this study, dust source regions and emission rates estimated by different parameterization schemes are compared for the 24^th May 2012, when dust event was reported in many stations of Iran. According to results, compared to the Shao scheme, the MBA predicted smaller number of source regions in Iran. The vertical dust flux simulated by MBA is higher than that by the Shao scheme in some common source regions, but lower in some others. After applying roughness and moisture correction functions to the threshold friction velocity, the <em>VDF</em> simulated by the MBA will be zero wherever the Shao schemes’ <em>VDF</em> is zero, but this is not necessarily true the other way. Both of the schemes simulated the maximum <em>VDF </em>at 11am. The highest value of the estimated <em>VDF</em> was  (for the grid centered at , ) by the Shao scheme, while it was   (or the grid centered at , ) for the MBA scheme. <br />The most influential parameters in the estimation of <em>VDF</em> are the roughness length and plastic pressure in MBA and Shao, respectively. Impacts of the friction velocity, and surface and soil- related factors on the estimated vertical dust flux, though complicated, should be considered simultaneously. In every source region, <em>VDF </em>increases when the soil moisture decreases and the friction velocity increases. Furthermore, the <em>VDF</em> variation during a day is mainly affected by daily changes of the friction velocity, because of negligible changes in the soil moisture during the day. Thorough evaluation of these schemes requires experimental and remote sensing data. <br />  <br />  https://jesphys.ut.ac.ir/article_36706_a0f6825d01e9cb9eecbaca4589bff5ae.pdf