IRANIAN JOURNAL OF GEOPHYSICS
Year:2016 | Volume:10 | Issue:3|Papers Count:11

Time-frequency analysis is an important technique in the seismic data involving nonstationary signal processing and interpretation. Due to the limitations of the Fourier transform in analyzing nonstationary signals, it cannot be used for time-frequency representation. Time-frequency transforms such as short-time Fourier transform, wavelet transform and S-transform are common tools in the study of nonstationary characteristics in the seismic data. Based on regularized least-squares inversion, Liu et al. (2009) have recently proposed a new method of time-varying frequency characterization of nonstationary seismic signals. In this paper, we applied the method of Liu et al. (2009) by designing an invertible nonstationary time-frequency transform called local time-frequency (LTF) transform (Liu and Fomel, 2010). This method generates time-frequency characterization without sliding windows. The LTF transform aims at depicting the nonstationary character of seismic data. This transform uses a Fourier basis to match the target signal under the regularized least-squares norm and provided an invertible time-frequency representation where are the Fourier coefficients and. The use of a non-stationary regression makes it possible for the coefficients to vary with x. In the linear case, they can be obtained by solving the least-squares problem: The solution where denotes the least-squares estimate of m and LT denotes the adjoint operator and S is a smoothing (shaping) operator, was introduced by Fomel (2007) using shaping regularization. The λ scaling in this solution controls the relative scaling of the forward operator L. The key idea is to minimize the error between the input signal and all of its Fourier components simultaneously using a regularized nonstationary regression (Fomel, 2009) with control on time resolution. The transform can provide LTF representations for common seismic data interpretation tasks such as Q factor estimation. Seismic waves lose energy by traveling through the earth due to a variety of phenomena such as attenuation. Attenuation refers to the loss of energy caused by phenomena other than geometrical spreading, and depends on the characteristics of the transmitting medium. Generally, attenuation is expressed in terms of quality factor (Q) which is a dimensionless parameter and is inversely proportional to attenuation coefficient. Experiments show that the quality factor is controlled by the elastic properties of formation and its fluid content. Hence, as one of the most important attributes in seismic exploration, it is used to directly identify hydrocarbon reserves on seismic sections (Hedlin et al., 2001). In this paper, we present a procedure using the differences in seismic reflection time-frequency spectra to estimate relative seismic attenuation in a reservoir of carbonate rock with fractures and voids. It is difficult to determine the seismic reflections at the top and bottom of the reservoir, required by the conventional amplitude ratio and the frequency spectra ratio methods. But in this study, we use the difference of seismic reflection time-frequency coefficients to estimate the relative seismic attenuation in such reservoirs. Gu and Stewart (2006) considered a special case where the incident spectrum has a Gaussian distribution, and discussed the reflection centroid of frequency downshift. Here, we give the analytical expression of the differences of reflection spectra corresponding to two frequencies with an attenuation coefficient. We show that the differences of reflection spectra corresponding to the two frequencies that are symmetrical to the centroid of frequency and separated by twice the incident signal’s standard deviation can be used to calculate the attenuation coefficient. As a time-frequency representation tool, the LTF transform of Fomel et al. (2010) is used to study seismic wave attenuation coefficient in synthetic and real field data examples.

In this study, using similarity laws and regression methods, numerical modeling of wave height estimations has been carried out for the data recorded by a station in the southern part of the Caspian Sea during the year 2008. Numerical simulations require powerful computers and their results may be accurate and reliable if proper physics is used. However, semi-empirical models can easily calculate wave heights but in practice are less accurate. In the latter methods, the surface temperature and instability effects have little impact on the effective wave height calculation. The nature of these methods is such that usually the effects of the instability of the surface layer are not considered in the calculation of wave heights. In experimental methods, the effects of variations of frictional velocity and the stability parameter are important factors which are considered as coefficients and directly applied in the equations. In this paper, the correlations between the wave height measurements and the friction velocity and the stability function are considered. These correlations are evaluated for exponential, logarithmic, polynomial and inverse polynomial functions. In order to calculate the friction velocity and the stability functions, a bulk aerodynamic model similar to the one used in Fairall et al. (2003) is used. It includes an iterative algorithm with an initial guess for the friction velocity and Monin–Obukhov length. The accuracies of the friction velocity and Monin–Obukhov length are then increased by iteration. The wave height measurements are analyzed by multi-variable nonlinear regression on the outputs of the bulk model. This analysis shows that, for functions that have an acceptable correlation coefficient (exponential, logarithmic, polynomial and inverse polynomial functions), the stability function can reduce the wave height error. The following conclusions can then be made: a) In the semi-empirical models, the effect of stability function can be directly entered into the equation and the uncertainty in the coefficient can be reduced. b) The stability function in the atmospheric surface layer has an important role in the air–sea interaction; hence its consideration can increase the accuracy of the wave height prediction. Results show that compared to the warm months (stable cases), in the cold months (unstable cases) the stability function had smaller effects on wave height prediction. Results of the analyses are summarized as: a) In the cold months (October–December) in the port of Amir Abad, the conditions are more unstable, while in relatively warm months (July and August) there are more stable conditions. b) The correlation coefficient between the friction velocity and the wave amplitude in the cold months is less than that for the warm months. c) The correlation between the stability fiction and the wave height in the cold months is smaller than that in the warm months. d) Regression error rates for the friction velocity in the cold months are higher than those for the warm months of the year. Wave height calculation error in December (with a maximum error of linear regression with friction velocity) decreased down to ten percent by calculating the nonlinear regression with a mixture of the third-degree polynomial, inverse polynomial and logarithmic functions for the stability. This reduction in the wave height calculation error is simply due to the effect of stability function. The friction velocity seems to have little impact on the results of the regression method. Although the results are only for the data of the year 2008 recorded in the Amir Abad Port, but the results are not limited to a specific time and place and show the impact of stability function on wave height estimations. The main achievement of this paper is a reduction of error in wave height estimation that is achieved by taking into account the stability function. However, much higher reduction in the wave hight errors cannot be expected, because a large part of the waves characteristics are caused by wind effects at different times and places. Since the effect of the surrounding area of Amir Abad Port has not been considered in the present work, only a reduction of the calculated wave height error was possible.

The East Atlantic-West Russia (EA–WR) teleconnection pattern is one of the low-frequency atmospheric phenomena that affects Europe and Asia, especially in the cold season. In this study, the effects of EA–WR on the climate of Southwest Asia are investigated using the NCEP/NCAR reanalysis dataset from 1950 to 2012 for winter months (Dec. to Feb.) and the monthly indices taken from the Climate Prediction Center (CPC). Because of the large zonal extension of the EA–WR teleconnection pattern, all data north of 20°N are taken into account for the analysis. In this paper, the method of composite maps is employed. Considering the critical positive and negative months, the average state of the troposphere is studied for each of the two phases from the synoptic viewpoint. To this end, a month is considered to be a critical positive (negative) month, if the monthly index of EA–WR is higher (lower) than the long-term mean value of EA–WR index plus (minus) its standard deviation. In this way, among 189 winter months during the 63-year period from 1950 to 2012, 26 positive critical months and 29 negative critical months are identified. For the part of analysis based on outgoing long wave radiation for which the data is available from 1974 onwards, there exist 20 positive critical months and 18 negative critical months. The composite map analyses include 500 hPa geopotential height and its anomaly, mean sea level pressure, 300 hPa wind field, 1000-500 hPa thickness, the outgoing long wave radiation and the Eady's parameter for the growth rate of baroclinic eddies: where is the Coriolis parameter or inertial frequency, is the buoyancy frequency, and is the magnitude of the vertical wind shear. The growth rate is evaluated and compared between the two phases at 800 hPa. In the critical positive months of EA–WR, in 500 hPa geopotential height field, there is a trough from the western part of Russia to the Middle East and a ridge over the eastern part of the North Atlantic. In critical negative months of EA–WR, however, there is a ridge over the western part of Russia, a trough over Europe and a dominantly zonal flow is observed over the Middle East. In the critical positive months, the subtropical jet stream over the southwest of Asia is stronger, and at the same time, the exit region of the polar front is extended to the border between Europe and Asia, which together with the wind field anomaly result in significant cold air advection to the northwest of Iran. Furthermore, the stronger subtropical jet over the Southwest of Asia, the southwest of Iran, and Saudi Arabia in the critical positive months is associated with increased amounts of Eady' s parameter for the growth rate and thus baroclinic instability. Overall, results point to a significant effect of EA–WR on the climate of the southwest Asia. Compared with the negative phase, the presence of a mid-tropospheric trough of geopotential height as well as the upper-tropospheric wind anomaly in the form of a stronger subtropical jet stream provide a better ground for the development and organization of synoptic systems and their impact on the climate of Iran in the positive phase. The dynamical effects mentioned are also helped by a more suitable lower-tropospheric moisture transport in the positive phase.

During the last decades, GPS has been used in many applications of geodesy, geophysics and surveying. The ellipsoidal height which is provided by GPS is a geometric height and lacks a physical meaning. So it cannot be used in most of engineering applications. Therefore, we need another type of height known as orthometric height which is one of the physical heights. The main advantage of leveling is its high accuracy but on the other hand, leveling measurements involve large amounts of time, cost and labor work. So another way to reach orthometric heights from GPS measurements is to determine geodetic height correcting surface of the area. Depending on the data availability and accuracy requirements, there are two basic approaches to transform the ellipsoidal heights to orthometric heights which include the gravimetric approach and geometric approach. The gravimetric approach uses the gravity data to determine the correcting surface. The geometric approach is to use the GPS/leveling data and interpolation methods in order to determine the correcting surface. In the present research, the geometric approach is applied to determine the geodetic height correcting surface from GPS/Leveling data. To do so, some reference points with known ellipsoidal and orthometric heights are used. Since there is a linear relation between the triplet of ellipsoidal, orthometric and geoid heights, one can calculate the geoid undulation by means of that relation and finally calculate the geoid undulation in every point of the area by the use of interpolation methods. First, basic definitions and equations of the radial basis function are explained. The method of radial basis functions is a global interpolation method used for interpolation of scattered data. This method was developed by Richard L. Hardy. Different kernels can be used as the kernel of radial basis function like multiquadric, Gaussian, thin plate spline and linear. In this paper, multiquadric and thin plate spline kernels are used to determine the correcting surface. The multiquadric method contains a parameter called shape parameter which is defined by the user and in actual applications, it affects the accuracy of the method. There are different ways to determine the shape parameter, among which the cross-validation procedure is explained. Then, the artificial neural network is defined as the other method to determine the correcting surface. Neural networks have different types but the multilayer perceptron neural network is the method commonly used in interpolation applications. Therefore, in this research, a three-layer perceptron neural network is used. Finally, as the case study, the GPS/leveling data of Tehran is analyzed. The ellipsoidal heights and orthometric heights of 147 benchmark points distributed all over the area are used in order to calculate the geoid undulations. And then the geodetic height correcting surface is determined using the multiquadric and thin plate spline and the artificial neural network methods. At last, their root-mean-square (RMS) values are compared with each other and the method with the lowest RMS is chosen as the most accurate method. The results show that the thin plate spline method leads to better results in Tehran.

In recent years, substantial amounts of research work have been devoted to using highly accurate numerical methods in the numerical solution of complex flow fields with multi-scale structures. The compact finite-difference methods are simple and powerful ways to attain the purpose of high accuracy and low computational costs. Compact schemes, compared with the traditional explicit finite difference schemes of the same order, have proved to be significantly more accurate along with the benefit of using smaller stencil sizes, which can be essential in treating non-periodic boundary conditions. Applications of some families of the compact schemes to spatial differencing of some idealized models of the atmosphere and oceans show that the compact finite difference schemes are promising methods for the numerical simulation of the atmosphere–ocean dynamics problems. This work is devoted to the application of the combined compact finite-difference method to the numerical solution of the gravity current problem. The two-dimensional incompressible Boussinesq equations constitute the governing equations that are used here for the numerical simulation of such flows. The focus of this work is on the application of the sixth-order combined compact finite difference method to spatial differencing of the vorticity-stream function-temperature formulation of the governing equations. First, we express formulation of the governing equations in dimensionless form. Then, we discretize the governing equations in time and space. For the spatial differencing of the governing equations, the sixth-order combined compact finite difference scheme is used and the classical fourth-order Runge–Kutta is used to advance the Boussinesq equations in time. Details of spatial differencing of the boundary conditions required to generate stable numerical solutions are presented. Furthermore, the details of development and implementation of appropriate no-slip boundary conditions, compatible with the sixth-order combined compact method, are presented. To assess the numerical accuracy, the Stommel ocean circulation model with known exact analytical solution is used as a linear prototype test problem. The performance of the sixth-order combined compact method is then compared with the conventional second-order centered and the fourth-order compact finite difference schemes. The global error estimations indicate the better performance of the sixth-order combined compact method over the conventional second-order centered and the fourth-order compact in term of accuracy. The two-dimensional planar and cylindrical lock-exchange flow configurations are used to conduct the numerical experiments using the governing Boussinesq equations. In this work, we used the no-penetration boundary conditions for temperature and no-slip boundary conditions for vorticity at walls compatible with the sixth-order combined compact scheme. The results are then compared qualitatively with the results presented by other researchers. Quantitative and qualitative comparisons of the results of the present work with the other published results for the planar lock-exchange flow indicate the better performance of the sixth-order combined compact scheme for the numerical solution of the two-dimensional incompressible Boussinesq equations over the second-order centered and the fourth-order compact methods. Hence, such methods can be used in numerical modelling of large-scale flows in the atmosphere and ocean with higher resolutions.

The variability of sea surface temperature (SST) is used as a valuable climate index for the prediction of precipitation in far and near areas from the sea. The prediction of SST in the north western of the Indian Ocean is the main goal of this study. This water region including 81 gridpoints with 2o×2o grid in the geographical location of 10-30N and 45-76oE. The SST was extracted from the National Oceanic and Atmospheric Administration (NOAA) for the period 1951–2007. The principal components analysis technique was used to identify the main patterns of SST and data reduction. The PCA performed was based on the correlation matrix. The number of row and column of the input file for the correlation matrix was, respectively, the number of months and gridpoints. The four principal components that explained 98% of the SST total variance were extracted. The first, second, third and fourth principal components explained 79, 8.9, 5.9 and 4% of the SST total variance, respectively. These four principal components as four regions over the area of interest were studied. The first, second and third regions were geographically located in 16-24oN and 58-72oE, 10-14oN and 48-76oE, and 14-16oN and 50-74oE. Also, the fourth region was the Persian Gulf. The spatial average of SST within each region was considered as a regional index. As the linear stochastic models, the “seasonal auto-regressive integrated moving average” (SARIMA) models were used to predict the monthly time series of the regional indices of SST patterns. The dataset for the monthly time scale for the 1951–2000 period was used to construct SARIMA models for each region. There is a linear trend in SST time series over three regions which indicate that the monthly SST over these regions is non-stationary. Since ARMA models prefer stationary time series data as their input files, a differencing procedure was considered as a smart approach for transforming these non-stationary series into the stationary ones. On the basis of the corrected Akaike information criterion (AIC) and significant coefficients, the best seasonal auto-regressive integrated moving average model was separately selected for each region. The auto-correlation function plots of the residuals for the selected models have indicated that the residuals are uncorrelated. The selected model for each region had a minimum value of AIC and its parameters were significantly different from zero. For example, SARIMA (1, 1, 0) × (1, 1, 0) 12 model was identified for SST time series over the northwestern parts of the study area. As the independent data of training period, the SST time series for the 2001–2007 period was predicted at lead times ranging from one to 12 months and then was evaluated. For example, the value of Pearson correlation between the observed and the predicted SST over the northwestern parts of the study area with SARIMA (1, 1, 0) × (1, 1, 0) 12 model for the test period (84 months) was 0.94. Also, the corresponding root mean square error was 0.46 degrees Celsius. In all of the regions, the correlation coefficient between the observed and the predicted SST for the independent dataset is higher than 0.9. Therefore, the time series models have a valuable ability in forecasting the monthly time series of SST in each region.

Induced polarization (IP) method is a main geophysical method in deposits exploration. As an extension of the IP method, the spectral induced polarization (SIP) has been used extensively in mineral prospecting and increasingly in environmental investigations, hydro-geophysics, archaeo-geophysics, bio-geophysics. The reason for this extensive use is that SIP measurements are sensitive to the low-frequency capacitive properties of rocks and soils. One major limitation of SIP method is electromagnetic (EM) coupling effect. In SIP method, the amplitude and phase components of the earth’s resistivity are measured in a frequency range typically from 0.001 Hz to 10 kHz. At low frequencies, the inductive coupling effects impact the spectrum Ohmic responses and normal polarization effect of the subsurface material. In SIP literature, there are three types of the EM coupling effects: the first is the EM coupling effect removal methods from SIP field data. In the second method, the mutual impedance of the earth is calculated using the Cole-Cole equation as IP dispersion of the earth. SIP data and mutual impedance are compared using an inversion algorithm in order to recover the earth IP parameters. Since the SIP method employs alternating fields using grounded wires, this method should be characterized as an EM method. The third method uses a current cable arrangement in order to reduce the EM coupling effects from SIP data. Many different models have been proposed for the description of the dispersive behavior of the IP. However, the most widely used model is the Cole–Cole model. This model describes the resistivity dispersion observed in the field data from areas with metallic mineral content. It is also used to estimate various subsurface properties of nonmetallic soil and rocks in IP frequency domain investigations (SIP). A multiple Cole–Cole model is typically a more general and proper model than a single Cole–Cole model for describing IP data with various dispersion ranges caused either by multiple-length scales in sediments or by coupling effects in the IP measurements. The Cole–Cole model parameters are widely used to interpret both time- and frequency-domain induced polarization data. Among many studies in which the Cole–Cole parameters are estimated from the SIP measurements on soils and rocks, a majority of them use least squares (deterministic) methods. The previous studies have shown that the geometry of an array such as electrode spacing (e.g., dipole–dipole electrode array) has an important effect on mutual impedance. In this study, by using the dipole–dipole array on a homogeneous polarizable half-space, the electromagnetic coupling effect on mutual impedance is investigated. The aim of this work is an investigation of the Cole–Cole parameters effects on the mutual impedance of a polarizable half-space. Since Sunde’s mutual impedance equation is widely used for an impolarizable earth (real resistivity), the effect of a polarizable earth has been less investigated. We use the Nyquist plot to show the mutual impedance response of theoretical and field data. The results show that if the Cole–Cole parameters including time constant (τ), frequency constant (c) or chargeability (m) of the half-space are small, the IP response is very small compared with the EM coupling response and thus the Cole–Cole parameters recovered from the inversion algorithms are less reliable. In practice, the above-mentioned terms occur when there are small particles of ore, extended grain size distribution of ore or low-grade ore in porphyry deposit. It is worth mentioning that the chargeability of the earth in environmental investigations also has a small value.

Precipitation is a vital component of the global water and energy cycle with large variations in space and time. The observational datasets that are based on meteorological stations data usually serve as the main sources of precipitation data. However, because of uneven distribution in space, such datasets may not be directly applicable to some problems. Furthermore, there are gaps in the data as there may be times for which the precipitation has not been recorded by some metrological stations for various technical reasons. In the last decades, several gridded precipitation databases have been developed by researchers or institutes. The main aim of creating these databases is to serve user requirements and solve the problems mentioned above. The even distribution in space of the gridded precipitation data and their availability are two very important factors. These databases are critical for many studies including climate change and numerical weather prediction (NWP) applications, management of water resources, agriculture, and disaster management. The Global Precipitation Climatology Centre (GPCC) has been established in the year 1989 at the request of the World Meteorological Organization (WMO). It is constructed by the Deutscher Wetterdienst (DWD, National Meteorological Service of Germany) as a German contribution to the World Climate Research Programme (WCRP). The precipitation data of GPCC are freely available via the website http: //gpcc.dwd.de at 2.5º × 2.5º, 1º×1º, and 0.5 º × 0.5 º resolutions. The aim of this research is to evaluate the accuracy of GPCC database over Iran by comparing it with two national databases, the Asfezari and that of the synoptic stations called Stations hereafter. The monthly precipitation data from the three databases including GPCC, Asfezari and Stations have been used from January 1962 to the end of December 2010. To evaluate the accuracy of the estimated GPCC precipitation data, first the spatial resolutions of the three databases have been synchronized by the nearest neighbor algorithm. The high-resolution database is converted to the low-resolution database in order to select spatial pixels and carry out the comparisons. Seven accuracy evaluation indices have been used. The results indicate a high temporal correlation between the estimated precipitation of GPCC and the observed precipitation by the Asefazri and the Stations databases. The results of applying accuracy evaluation indices to the precipitation time series show that in addition to high temporal correlation, quantitatively the estimated precipitations are also very similar to the observed precipitations. Although in some regions the estimated precipitation values are contaminated with bias, but overall the estimated precipitation error is low compared to the total precipitation received. In a spatial viewpoint, the highest accuracy is observed over the western parts of Zagros mountain range and the northeast of the country. Over these regions, the index of agreement and the coefficient of determination are close to unity. The highest relative root mean square (rms) is observed over the dry interior regions and the Lut desert. The relative rms is low in the regions with high precipitation when compared with the rest of the regions. In a temporal viewpoint, the highest correlation between the precipitation time series is observed in rainy months. Based on the results from the Nash–Sutcliffe efficiency index, over most regions of the country, using the estimated precipitations by GPCC is preferable to applying the mean precipitation amounts. The results of this study confirmed the finding of other researches about the accuracy of the estimated precipitation by GPCC database.

The assessment of Guilan’s meteorological radar data (PPI and SRI) during autumn 2013 and winter 2014 shows the presence of a large number of pixels, mostly located in Guilan mountainous areas, involving incorrect reflectivities. For example, the daily average of incorrect reflectivity on the 26th of December 2014 calculated by PPI (plan position indicator) products of Guilan’s meteorological radar resulted in 47.014 dBZ in Lahijan synoptic station which would be equal to 31.6 mm/hr rain rate through Marshall-Palmer Z-R relationship. The importance of these results could be appreciated when we note that Lahijan synoptic station had a daily gauge rainfall height equal to zero on the same day. Thus, according to the position of incorrect reflectivities of Guilan’s meteorological radar, this study aims to evaluate the ground clutter influence on incorrect reflectivity appearance. The method used is to study the overlapping of the beam blocking map and the daily average of incorrect reflectivity map by changing the elevation angle of the radar beam at a clear air day. The results show that the ground clutter is the main factor responsible for incorrect reflectivities in %75.34 pixels of the image. Using comparison of the mean distance of the blocked pixels and the pixels containing the ground clutter to the radar antenna leads us to conclude that the real amount of the elevation angle is equal to -0.20°. The method adopted by assuming 5 km donut-shape regions around the radar antenna and using it to determine the radar elevation angle is one of the innovations of this study. Furthermore, to ensure the working of the actions undertaken for the clutter contamination mitigation, one requires a numerical threshold to determine the effectiveness of the actions. For this purpose, therefore, a method is developed based on the correlation between the radar’s 24-hr rainfall heights and their respective gauge data. The method provides a maximum threshold in dBZ unit for the clutter contamination while Z=200 R1.6 is applied. The numerical thresholds resulted from the application of the method for the meteorological stations of Rasht (the Airport), Rasht (the Agriculture Faculty), Anzali, Talesh, Lahijan, Jirandeh, Masouleh, Deylaman and Manjil are equal to -3.748, -1.638, -6.074, 14.952, -3.611, -1.482, 12.466, 4.872 and 6.852 in dBZ, respectively. Also, a value of 9.536 dBZ is obtained for all the stations considered together. The latter value can be extended to all pixels of the radar image. As the conclusion, in spite of the fact that on the 26th of December 2014 as a non-rainy day, the average clutter was considerable, a close correlation between the radar’s rainfall heights and their respective gauge data is observed during the 15 days studied. The Pearson correlation coefficients (Pcc) in the meteorological stations of Rasht (the Airport), Rasht (the Agriculture Faculty), Anzali, Lahijan, Jirandeh, Deylaman and Manjil are significant at 1% statistical significance level, and also Pcc values in the meteorological stations of Talesh and Masouleh are larger than that can be excluded. Therefore, the clutter contamination of Guilan’s meteorological radar data cannot be considered an impediment to its application in practical purposes. However, for further progress, a higher degree of accuracy would be required.

Delineation of homogeneous precipitation sub-regions featured with different time variabilities is very important for large countries such as Iran, which are characterized by complex topography and different climates. Very rare efforts have been devoted to identify modes of monthly precipitation variability in Iran and delineating sub-regions having different temporal variabilities of precipitation. On the other hand, most studies of precipitation regionalization in Iran have used very limited and unevenly scattered stations across the country; thus making it necessary to identify the most realistic precipitation sub-regions for Iran using almost all available stations. As such, 155 synoptic stations with relatively regular distribution over Iran, mostly having full data records for the 25 years common period of 1990–2014, were used for identifying an updated precipitation regionalization of the country. The cubic root transformed monthly precipitation of the considered stations were used as input for an S-mode principal component analysis (PCA) applied to the inter-stations correlation matrix (300×155) that is composed of 155 stations and 300 cubic root transformed monthly precipitation. The computed Kaiser–Meyer–Olkin measure of sampling adequacy for the considered matrix with a value of 0.98 indicates that the considered matrix is marvelous for a PCA application. The first five leading significant PCs accounting for approximately 80% of total variance of the dataset were considered for further analysis based on the Scree plot and the sampling errors of the PCs (North et al., 1982). To better characterize the underlying spatial structure of the considered data matrix, the retained PCs were then rotated using varimax orthogonal criteria. The five leading varimax rotated loadings were mapped to present spatial modes of monthly precipitation variability across the country and precipitation sub-regions borders were delineated using the maximum loading value approach (Comrie and Glenn, 1998; Miller and Goodrich, 2007; Chen et al., 2009). The maps of varimax rotated loadings well represent areas characterized by different modes of precipitation variability and regimes. The five precipitation sub-regions identified using maximum loading values of the varimax rotated components are the Caspian Sea region, the northwestern, the western, the central-eastern, and the central-northeastern of the country. The Caspian Sea region featured with maximum precipitation in autumn and relatively regular distribution of precipitation throughout the year includes the coastal areas of the Caspian Sea and the northern faces of the Alborz Mountain in northern Iran. The north-western sub-region is distinguished from the rest of the country for its identical precipitation regime characterized by maximum precipitation in spring and relatively uniform precipitation all over the year. The three remained precipitation sub-regions of Iran are characterized with a much shorter rainy season, which maximizes in the winter time. The western sub-region encompasses mountainous areas of western Iran as well as the lowlands of the southwestern country. The central-eastern sub-region differs from the western sub-region due to its shorter rainy season and much lower precipitation values in all of the months, but similarly, its maximum precipitation occurs in January. Finally, the central-northeastern precipitation sub-region receives its maximum precipitation in March as opposed to the two aforementioned sub-regions which peak in January. The independence of the identified precipitation sub-regions was examined by applying the Kolmogorov–Smirnov non-parametric test to the regional anomalies of annual precipitation series; the result proved that all the sub-regions are statistically different at 99% confidence level. The identified precipitation sub-regions can serve as a tool for a better water resources management in the country.