مجله ژئوفیزیک ایران
سال پنجم شماره 2 (پیاپی 9، تابستان 1390)

The normalized full gradient (NFG) method defined by Berezkin (1967, 1973 and 1998)is used for downward continuation maps. Analytical downward continuation is a methodof estimating the field closer to the source, and, consequently, it results in a betterresolution of underground rock distribution. However, the usefulness of this process islimited by the fact that the operation is extremely sensitive to noise. With noise free data,downward continuation is well defined; it is unnecessary to continue below the sourcelevel. In the presence of noise, the amplification of high frequencies is so strong that itquickly masks the information in the original profile. Low-pass Fourier filtering, whilesuppressing such noise, also blurs the signal, overcoming the purpose of sharpening bydownward continuation.Despite the above-mentioned problems, most geophysical experts have long beeninterested in this technique because of its importance to mineral exploration. Furthermore,this method is a fast and cheap way to determine the initial depth of the subsurfacefeatures, especially where there is no other geophysical or well-logging data. A goodanalytical downward continuation process could provide subsurface general images,allowing an enhanced interpretation. Also, analytical downward continuation has theability to determine accurately both the horizontal and vertical extents of geologicalsources.The NFG method nullifies perturbations due to the passage of mass depth duringdownward continuation. The method depends on the downward analytical continuation ofnormalized full gradient values of magnetic data. Analytical continuation discriminatescertain structural anomalies which cannot be distinguished in the observed magnetic field.It can be used to estimate location, depth to the top, and center of the deposit that isapplied also for detecting oil reserviors and tectonic studies. One of the importantparameters to estimate the accurate shape of the deposit is true selection of the harmonicnumber. In this paper, the correct harmonic number is determined and then this methodwill be tested for noise-free and noise-corruption synthetic data. Finally, 2D and 3D ofthis method are applied to real data from Zanjan’s Morvarid mine.The Morvarid mine is located at a distance of 23 km southeast of Zanjan, near Aliabadvillage. Hematite and magnetite are the major minerals in the mine in which magneticexploration has been performed to find anomalies around the mine, near outcrops of theregion. From the IGRF model, the geomagnetic field of this region is calculated at 47680nT (inclination = 54.7°, declination= 4.5°).The results of the method were satisfactory for a dike. After applying the method, ageneral view of deposit structure was obtained that helped us to make a fast interpretation.2D and 3D results of NFG method satisfied the results of trenches and drilling in theregion are under study, as well. The final result of this method shows that the depositbegins from the low depth to approximately less than 100 meters. This modeling yeildedsatisfactory results for drilling in the region. The results of the drillings show that thelowest depth of the deposit is about 6 meters.Finally, 2.5D modeling of th edeposit was done using Modelvision software that hadthe same results. This method can easily be applied for gravity and microgravity data.

The outer part of the upper atmosphere, the ionosphere, is where free electrons formed by solar X-rays and extreme ultraviolet (EUV) radiation play an important role in radio and satellite communication. Knowledge of the critical frequency and maximum electron density of the F2 ionospheric layer distribution is crucial for ionospheric storm studies and the estimation and correction of propagation delays in telecommunications. It has proven that the F2 ionospheric layer can significantly affect the propagation of radiowaves. The ionosphere can be used to reflect radio signals over long distances. Indeed,the ionosphere is an efficient reflector with frequencies below approximately 30 MHz.Also, it can be useful to study the variability of the F2 ionospheric layer to show that thislayer is affected primarily by space weather phenomena, mostly of solar origin such as thesolar zenith angle, solar ionizing radiation, and the solar cycle. These extraterrestrialphenomena control the ionosphere from above. Some other waves enter the ionospherefrom below and can cause some significant changes in the F2 ionospheric layer.The International Reference Ionosphere (IRI) model is a standard model of the ionosphere supported by the Committee on Space Research (COSPAR) and InternationalUnion of Radio Science (URSI). The IRI model has many practical applications in HighFrequency (HF) predictions. The IRI model offers a description of the average ionosphere. It is a mathematical description of the ionosphere as a function of location,time, altitude, solar activity and geomagnetic activity. Periodic updates to this model areessential to maintaining its prediction ability. A large number of independent studies haveevaluated the IRI model in comparisons with direct and indirect ionospheric measurements, those not used in the model development. A favorable comparison withIRI model is often one of the major goals of ionospheric teams all over the world.In this work, to evaluate the latest available ionospheric model, IRI-2007, we haveobtained hourly monthly values of foF2 and NmF2 over the Tehran area (35.4N,51.2E,52.7dip) during low solar activity period, in which the Rz12 (12-month running averagesunspot number) varies between 7.7 and 15.3 from July 2006 to June 2007. Datameasured using the IPS-71 at the ionospheric station at the Institute of Geophysics at theUniversity of Tehran (35.4N, 51.2E, 52.7dip) were used to perform the calculations.Subsequently, the observed critical frequency and maximum electron density of the F2ionospheric layer is compared with IRI-2007 model predictions. To run the IRI2007model, the URSI and CCIR coefficients were used.Our study shows that values of the foF2 and NmF2 parameters have the highest valuesduring the daytime hours and the lowest values occur at pre-sunrise hours. Our studyshows that the IRI-2007’s ability to predict semi-annual anomalies with a maximumelectron density of the F2 ionospheric layer are most accurate in the summer and less soduring the autumn and spring.In general, the predictions obtained with CCIR and URSI are similar. Our results showthat IRI-2007 can successfully predict the critical frequency and maximum electrondensity of the F2 ionospheric layer. Additionally, our study shows that differencesbetween the foF2(OBS) and the foF2(IRI-2007) remains below 12% during all seasons.The best agreement occurs during the summer and winter, and the largest differences are observed in the spring and autumn. The average percentage deviation of a full yearregisters at approximately 6.5% for CCIR and nearly 8% for URSI coeffecients.Moreover, our results show that the percentage deviation between the NmF2(OBS) andthe NmF2(IRI-2007) remains lower than 21% during all seasons. The total averagepercentage deviation of a full year is approximately 13% for the CCIR coefficient andnearly 16.5% for the URSI.

To understand the physics of the F2 ionospheric layer and for ionospheric radio-wavepropagation studies, accurate values of the peak height of the F2 ionospheric layer(hmF2) are needed. The peak height of the F2 ionospheric layer is the most importantparameter after peak density in high frequency (HF) propagation prediction. Diurnalvariation of hmF2 is controlled by solar radiation and thermospheric winds.Thermospheric winds cause pressures that change the ionosphere at different altitudes.Therefore, hmF2 is a key parameter in studying ionospheric dynamics. In addition, theinteraction between radio waves and charged particles of ionospheric plasma can causesignificant changes in the phase, amplitude and polarization of navigation signals. All ofthese changes depend on the value of the total electron content (TEC) in the ionosphere. Data measured by IPS-71 at the ionospheric station at the Institute of Geophysicsof the University of Tehran from July 2006 to June 2007 were used to perform the calculations. As the real height is difficult to acquire directly from ionograms, the Shimazaki and Bilitza formulations are used to calculate the peak height of the F2 ionospheric layer for the period from July 2006 to June 2007. The most commonly employed mathematical function that depicts the electron density profile in the ionosphere is the Chapman function. In this work, we have used the method presented by Ezquer et.al (1992), in which two Chapman profile expressions are utilizedfor topside ionosphere and bottomside ionosphere, to obtain TEC from ionograms.Furthermore, in the present work, the international reference ionosphere (IRI-2007)model with CCIR coefficients was used to predict the peak height of the F2 ionosphericlayer and the TEC for the period from July 2006 to June 2007. The IRI-2007 model is anempirical models that is widely used to show ionospheric changes. The IRI model is astandard model of the ionosphere supported by the Committee on Space Research (COSPAR) and International Union of Radio Science (URSI). The IRI model has many practical applications in High Frequency (HF) predictions. The IRI model is based on a mathematical description of the ionosphere as a function of location, time, altitude, solaractivity and geomagnetic activity. It has two options for the prediction of the F2, F1 and E layer critical frequencies (CCIR and URSI coefficients). The integration of two Chapman profiles shows that almost 2/3 of the TEC of theionosphere is in the topside ionosphere and almost 1/3 is in the bottomside of ionosphere.Additionally, the daytime values of the TEC show intense fluctuation, while theybecomes smoother in the nighttime. Furthermore, during the daytime values of hmF2computed by the Bilitza method are lower than those calculated by the Shimazakimethod. However, during the nighttime the values of hmF2 computed by Bilitza methodare the as same as those computed with Shimazaki method. In addition, our study showsthat around midnight, the values of hmF2 are greater than those of the daytime.Conversely, there is a downward shift in the F2 peak near sunrise. From morning toafternoon, the hmF2 increases again. In the evening an upward drift of the F2-peak isclear.

There is a set of variables affecting the results of an inversion process. Some effects arenegligible while others are significant. It is important to know these essential variables toset them appropriately. In this study, the effect of change in the number and the thickness of the layers, the trade-off parameter, the number of frequencies and frequency content of the signals used for the inversion are studied. Occam’s inversion is employed toreconstruct the resistivity and/or conductivity of the assumed layers.The thickness and number of layers can affect the recovered model. The best way tochoose the number of layers is to increase it until the fitness error flattens. The errordecreases, significantly with the inclusion of the first additional layer then in increasinglylesser degrees with subsequent additions, until it reaches a relatively constant level. Thenumber of layers at which it first reaches this constant level is assumed to be the bestnumber. As the number of layers is increased, the original model is modified slightly. Onthe other hand, the EM signal strength drops exponentially with distance. Hence, it isbetter to increase the thickness of the lower layers to overcome this decrease. Experienceshows that logarithmically increasing thicknesses result in better models. However, if thesufficient number of layers has been chosen, the change in thickness cannot affect therecovered model significantly because a thick layer can be estimated as a sum of thinnerlayers. The trade-off parameter has a critical effect on model parameters. It controls whether more weight is given to minimize the norm of the data misfit or model norm. When this parameter is large, the inversion process tends to produce a smoother model and the data misfit becomes less important. On the other hand, when it is small, the data misfit takes over regardless of a priori information. Generally speaking, the best trade-off parameter is considered to be that at which the data fit is satisfying. The results of this paper show that the effect of the changes on the trade-off parameter is significant and if one chooses it inappropriately, the estimated model will fail to represent the true model.It is clear that the inversion response depends on the number of data used. The moredata used, the better the model produced. The values of frequencies used are alsoimportant because high frequency data can have information only about surface layerswhile low frequency data have more information from deeper layers. The effect of thenumber and values of the frequencies was studied by changing or removing high/lowfrequency. The results show that the change in frequency content causes a considerableeffect on the resolution of inversion. We also studied the effect of the changes of these parameters on the results of inverting a real data set and the results confirm the effects described in the models recovered. In brief, we can say that the effects due to changes in the number and the thickness of layers are negligible, but the effects of changes in the trade-off parameter and the number of frequencies are significant.

Ground vibration is an inevitable consequence of blasting operations in open pit miningand civil engineering projects. A large amount of energy is released in every blastingexercise in the form of wave propagation, which can cause serious damage to thesurrounding environment. Safe design of blasting for underground structures andtunneling are nowadays carried out mainly based on the permissible peak particle velocity(PPV), which has become a major criterion in setting out the appropriate safety standardsenvisaged to prevent structural damage near such operations. Many published studieshave produced correlations for PPV prediction using experimental data The GotvandOlya Dam is one of the major civil engineering projects in the last thirty years in Iran, constructed as the last dam on the Karoon River with the aim of increasingenvironmentally friendly electrical energy of the country, seasonal flood control in theprovince, and a supply of agricultural water for the strategic farming grounds ofKhoozestan. The geological formations of Gotvand include the two main formations ofBakhtiyari and Aghajari. The existence of discontinuities in these formations as well asjoint sets and their low inherent strength caused a number of challenges in blasting designoperations. This dam is 30 km north west of Shoshtar in the Khoozestan province, and 12km away from Gotvand. In this study, 16 records of 3 components produced from 4 blastswere obtained using 4 seismographs of the type PG-2002 at Intac and surge tanks of theGotvand Olya Dam in order to investigate their effects on underground structures. Thethree-component seismometers were of the GS-11D type and the records were analyzedusing DADISP software. The explosion materials used were ANFO and dynamite, andthe maximum weight of the explosions per delay in the Intac and surge tank were in therange of 32.0 to 343.3 kg. The positions of the seismometers were in the range of 56.2 to 145.5 meters away from the centre of blast blocks. The permissible peak particle velocity was taken from the Korean Institute of Geology, Mining and Materials. Using the scale distance, an exponential equation based on the cubic root of the charge weight per delay and coefficient correlation of 85% was proposed for predicting the PPV. In thisrelationship, the geological coefficients of the Gotvand area were estimated to be 20.339and -3.08. Using a genetic algorithm, an improved coefficient correlation between theexperimental and predicted PPV of 89% was obtained. The results from the geneticalgorithm applied for the coefficients of a, b, and n in the equation of PPV=b(d/wn)a were -9.286849, 9.769703 and 0.745959, respectively. Considering this proposed relationship,the minimum distance of 56.2 m to the centre of the blast block and the permissible peak particle velocity of 254 mm/s for a 10 day old concrete, the maximum charge weight per delay was found to be 254.4 kg.

The quality of seismic data varies tremendously, from areas where excellent reflections(or refractions) are obtained to areas in which the most modern equipment, complex fieldtechniques, and sophisticated data processing do not yield usable data. Between theseextremes, lie most areas in which useful results can be obtained. Seismic records aregenerally affected by various types of noise, such as ground rolls, multiples, randomnoise, and reflection and reflected refraction from near surface structures. Random noiseresultiing from random oscillation during data acquisition is one of the most importantand harmful noises that exists in seismic data over all times and frequencies. Many effortshave been made to remove this type of noise from seismic data. The predictive filter is anordinary method commonly used for random noise attenuation from seismic data. Thisfilter can be used in various domains, such as the f-x domain (Haris and White, 1997) andthe discrete Cosine domain (Lu and Liu, 2007). Jones and Levy (1987) removed events which were not coherent trace-to-trace events by means of the Karhunen-Loevetransform. The empirical mode decomposition (EMD) method is an algorithm for the analysis of multicomponent signals that breaks them down into a number of amplitude and frequency modulated zero-mean signals, termed intrinsic mode functions (IMFs). An IMF must fulfill two requirements: (1) the number of extrema and the number of zero crossings are either equal or differ at most by one; (2) at any point, the mean value of the envelope defined by the local maxima and the envelope defined by the local minima is zero. In contrast to conventional decomposition methods such as wavelets, which performthe analysis by projecting the signal under consideration onto a number of predefinedbasis vectors, EMD expresses the signal as an expansion of basic functions that aresignal-dependent and estimated via an iterative procedure called sifting. Apart from thespecific applications of EMD, a more generalized task in which EMD can prove useful issignal denoising (Kopsinis and McLaughlin, 2009). When EMD is used for denoising, theproblem is to identify properly which IMFs contain noise characteristics. Certain modeswill consist mainly of noise, whereas other modes will contain both signal and noisecharacteristics. In the case of white Gaussian noise, the noise-only energy of the modesdecreases logarithmically. The first mode, carrying the highest amount of noise energy,will consist mainly of noise, and the effect of noise should gradually weaken with highermodes. In this paper, a new signal denoising method based on the empirical modedecomposition framework is used to suppress random noises in seismic data. A Noisysignal is decomposed into oscillatory components (IMFs). The empirical modedecomposition denoising method involves filtering each intrinsic mode function andreconstructing s the estimated signal using the processed intrinsic mode functions. Thedirect application of wavelet-like thresholding to the decomposition modes is, inprinciple, wrong and can have catastrophic consequences regarding the continuity of thereconstructed signal. This arises as a result of the special attributes of IMFs; namely, theyresemble an AM/FM modulated sinusoid with zero mean. Consequently, we used theinterval thresholding method instead of direct theresholding method to denoise seismicsignal. the efficiency of the proposed method was tested on both synthetic and real seismic data. In every case, results show that the denoising algorithm can suppress random noise significantly.

Since, the Monin-Obukhov similarity theory (1954) (MOST) has been investigated innumerous studies, especially for the atmospheric surface layer. This theory provides alow-cost method to estimate non-dimensional turbulence intensity and dimensionlesswind gradients that are the bases for the dispersion models. Therefore, this theory hasbeen examined for many surface types, such as water, forest, rural, and urban areas tovalidate the predicted similarity relations.Since the similarity relations were often found for neutral stratification, somemodifications need to be made in order to make them applicable for the other stabilityconditions and, hence, some stability correction functions are applied.The hypothesis of constant surface fluxes with height as a basic assumption of MOST,makes it unpractical for terrains such as forests and urban areas, over which surface fluxesvary horizontally and with regard to height. For these terrains, the local scaling approachwas applied to solve this problem. Similarity functions derived using local scaling wereexamined in numerous experiments and showed good agreement with observationalresults.In this paper, the applicability of MOST was studied over an urban area with complextopography (Tehran) and with a rough sublayer as the lowest part of the surface layer. Forthis, we used data measured by a Sodar (model PA1) located at the Institute ofGeophysics of the University of Tehran for the first 8 months of 2007. Sodar measuredthe speed and direction of wind, standard deviations of three components of the wind, andmomentum fluxes. Also, data from a 100 m tower 150 m away from the Sodar positionwere used. This tower is equipped with four 2D sonic anemometers at 2, 10, 45 and 100m, which measure the speed and direction of wind and the temperature. To validate Sodardata, the u and v components of the wind at 45 m from Sodar were compared with thosefrom the tower at 45 m. The results illustrate good agreement and show similar trends.The first aim of this research was to determine whether the similarity function can befitted to the non-dimensional wind gradient and standard deviations of three componentsof wind. Toward this end, non-dimensional wind gradients and dimensionless turbulentintensities were studied against the stability parameters which were calculated using thelocal similarity approach. Results show that MOST with the local scaling approach can beapplied over this terrain. Therefore, similarity functions for these parameters werecalculated.The second aim of this study was to compare calculated emperical constants for thisterrain with that of previous results that were made over various surface types and fordifferent stability conditions. To achieve an accurate comparison, neutral stratificationswere also studied.The results indicate that emperical constants for the non-dimensional wind gradientand turbulent intensities over Tehran and in the rough sub-layer were larger than those ofother surface types, including those over flat urban areas. This outcome reveals that theknowledge of the effect of Tehran topography, which may have a great influence on airpollution for this city, may be important.Different universal similarity functions and large emperical constants imply thatapplying results from other terrains to predict wind and air pollution over a city as largeand heterogeneous as Tehran may lead to significant errors. The estimated constants forthe empirical universal similarity functions found in this work appear to be necessary forapplications in dispersion studies for Tehran.

Self-potential (SP) prospecting is one of the oldest geoelectrical methods, and it is stillused in many fields of applied geophysics. The SP measurements refer to that part of thenatural electric field which is stationary in time, or nearly so, and whose current sourcesystem is generated and sustained by phenomena occurring underground withingeological structures. From the physical point of view, the common aspect of the manysource models is an electrical charge polarization is setup, which is responsible forelectrical current circulation in conductive rocks. Hence, the detected SP anomalies aresimply the surface evidence of a more or less steady state of electrical polarization (Abdelrahman et al., 1999). The two-dimensional inclined sheet model can be useful inquantitative interpretation of SP data measured in a small area over a buried structure.Estimating the depth upper edges and the dip of a two-dimensional inclined sheet fromself-potential anomalies measured perpendicular to the strike of the structures has drawnconsiderable attention. The advantages of the fixed geometrical methods over the 2-D and3-D continuous modeling are that they require no current density, and the resistivity andrough estimation on depth and the interpretation of isolated SP anomalies is fast andaccurate. Many of the geological structures in mineral exploration can be classified intofour categories: spheres, vertical cylinders, horizontal cylinders and inclined sheets. Thesefour simple geometric forms are convenient approximations of common geologicalstructures often encountered in the interpretation of SP data. The models may not beentirely geologically realistic, but the approximate equivalence is usually sufficient todetermine whether the form and magnitude of calculated SP effects are close enough tothose observed to make the geologic interpretation reasonable. Evaluation of the depthand the width of a 2D plate from self-potential (SP) anomalies measured perpendicular tothe strike of the structure has attracted considerable attention. Several methods have been proposed and discussed by many researchers to interpret the self-potential anomalies caused by 2D inclined plates, including, for example, logarithmic curve matching(Murthy and Haricharan, 1985) and the analysis in the frequency domain using theFourier transform (Rao et al., 1982). In this paper, a new method has been developed toestimate the depth and the half-width of a 2D inclined plate with a self-potential anomaly.This method is the second- moving average for successive window lengths to determinethe depth and half-width of two-dimensional inclined plate. For a fixed window length,the depth is determined iteratively using the solution of a nonlinear equation with theiterative Newton Raphson method for each half-width value (Stanley, 1977). Thecomputed depths are plotted against the half-width values representing a continuouswindow curve. This method can be used for two model and synthetic data with andwithout random noise. After adding 10% random error in the synthetic data, themaximum error in model parameters is less than 3%. The problem of determining theplate parameters from self-potential data of short profile length can be solved using thesecond-moving average method. The method of window curves has been extended usingthe second-moving average self-potential anomaly to determine the depth and half-widthof an inclined plate.

In the study of atmospheric phenomena, different types of parameters are generallyconsidered. The type of parameter applied depends mostly on the scale and severity of the phenomena under consideration. Weather systems with fast rotating wind, such astropical storms or tornadoes, are among the most destructive weather-related phenomenon. Study of these kinds of phenomena and recognition of the causes of theirintensification or weakening are very important. Quantities that are measured in this typeof system vary and the calculations yield different results and properties depending on themethod used. There is a direct relationship between type of system and severity to thethermodynamic instability of the atmosphere in the region of the occurrence. It is typicalto use certain parameters for studying atmospheric instabilities, such as vorticity,potential vorticity. However, these parameters are gauges for measuring the power of asystem related to the value of the instability of a system suitable for a straight flow. Tomeasure the energy of activity of system with a severely rotating flow it is suggested touse another quantity, that of helicity. Thus, surveying sources and sinks of helicity ishighly important in dealing with vortices and hurricanes and rotating systems in general.Besides helicity, helicity flux is also a very valuable tool for studying rotating flows. Themathematical relationship between helicity and helicity-flux are presented in this study.Also, a new method of computing helicity flux is discussed and the values of helicity intwo methods are compared.In theory, the dissipation of helicity occurs in the boundary layer. In fact, thedissipation aspect of helicity involved in the frictional process in the boundary layer ismuch higher than the production of it in buoyancy process. Consequently, there is aseparating surface in the boundary layer above the surface of which the production ishigher than the loss by friction effect, and below it the friction is much higher than thatproduced by buoyancy. The helicity produced in the upper layer moves downward intothe area wherein the turbulent viscose force of the surface is dissipated due to friction.Since the stationary condition of the vortex is used in calculating helicity,, the downwardflux of helicity is a highly accurate parameter of the stationary turbulent vortex.This study calculates the downward helicity flux, sources and sinks of helicity of thehurricane Gonu. Data of hurricane Bonnie were utilized as both a comparison andreference for the study of hurricane Gonu. The results show that dynamical buoyancy isthe main factor in producing helicity and surface friction and the main dissipating factor.Therefore, an increase in dynamical buoyancy or reduction in dissipative friction resultsin the regeneration or enhancement of a rotating system. Furthermore, the increase (ordecrease) of the magnitude of helicity results in an increase (or decrease) in the height ofthe layer at which the maximum helicity flux occurs. It is suggested that the maximumhelicity flux occurs at the top of viscous turbulence boundary layer.

One of the existing problems in geodetic control networks for computing displacement vectors is finding a fixed station which has no deformation within the observation interval. This problem comes from the datum-dependency characteristic of the displacement vector. Accordingly, extraction of invariant parameters of the deformation tensor in geodetic control networks which are independent from datum has an effective role in the accurate interpretation of results in geodynamic studies. In height networks, change of the curvature tensor and its associated invariants have been introduced for assessing vertical deformation. The summation and difference of elongation of change of the curvature tensor are two common invariants of this tensor. Additionally, changes of the Gaussian and mean curvature parameters are two other key invariants with physical interpretations that are used to describe deformation behavior. Although many methods have been proposed to calculate deformation tensor fields on the earth’s surface, few refer to the actual surface of the Earth. Most of these methods formulate the problem on reference surfaces such as projection planes or spheres and, consequently, their results suffer from possible effects of inaccuracy and incompleteness of the mathematical models of projections. In the present study, we used a method of differential geometry that allows deformation analysis of the actual surface of the Earth for a more reliable and accurate estimate of the surface deformation measures. The method takes advantage of the simplicity of 2-dimensional spaces versus 3-dimensional spaces without losing or neglecting information and effect of the third dimension in the final results. Khorasan is a large province in the northeast of Iran where over-extraction of water resources for industrial and agricultural purposes has caused an extensive subsidence in some cities such as Mashhad, Neyshabour and Kashmar. Since the employment of first-order precise leveling network of Iran, which has been utilized twice for over two centuries in these cities, precise leveling has been one of the observations methods for measuring the subsidence. On the other hand, GPS and InSAR observations in these areas revealed the extent and magnitude of the subsidence. In this paper, we estimated the invariant parameters of the tensor of curvature using precise leveling observation to study the height deformation behavior of the province of Khorasan in these subsidence areas. We computed the invariant parameters of change of curvature tensor relative to two distinct datums to evaluate their datum-independency characteristic as well as studying subsidence behavior. Obtained results relative to 2 datums show the unique value and pattern for each parameter and reveal the subsidence area. The computed maximum mean curvature change was approximately, and the maximum Gaussian curvature change was approximately, which confirms the high rate of subsidence in this area. This application reveals the capabilities and strengths of the proposed method and suggests the superseding of these invariant parameters with other geodetic network adjustment results, especially in areas in which a fixed datum is undefined.