Institute of Geophysics, University of TehranJournal of the Earth and Space Physics2538-371X48220220823Improving diffractivity attribute to image faults using tapered local semblance in post-stack domainImproving diffractivity attribute to image faults using tapered local semblance in post-stack domain2432608690910.22059/jesphys.2022.326374.1007334FAMohammadHosseiniPh.D. Student, Islamic Azad University, Science and Research Branch, Tehran, Iran0000-0003-4145-663XHamid RezaSiahkoohiProfessor, Department of Earth Physics, Institute of Geophysics, University of Tehran, Tehran, Iran0000-0002-9227-2972Journal Article20211009Diffractions carry useful and important information about subsurface features such as unconformities, faults, pinch-out, and so on. On the other hand, most of the information is encoded in diffractions. Polarity reversal across diffraction move out curves that are generated from fault’s edges, is a great challenge in seismic diffraction imaging. For the last few decades, several conventional methods in the pre-and post-stack domains, have been carried out for the diffractions characteristics and their locations. But most of these methods were not able to deal with the polarity reversal for diffraction imaging, some of them were time consuming, and needed to have some correction to deal with polarity changes, especially in diffraction caused by fault’s edges. Despite a large amount of research that has been carried out on diffraction imaging, very few studies have been devoted to the challenge of the polarity reversal across move out surfaces. We used the semblance function along with the hyperbolic move out curves for the diffractions that their travel times have been calculated using the double-square-root equation. As we know, using both semblance and Kirchhoff migration for diffraction imaging from fault’s edges, without which taking the polarity reversal into account would fail. This is caused by presence of same number of positive and negative wavelets in the diffraction move out curves. For solving this problem, we divided the global scanning window along hyperbolic move out surfaces into several subdivided window and the local semblance measurements over the sub-windows were performed separately. Every point in image domain is considered as a diffraction point that we call this points as image points. The final semblance measure at each image point is calculated by averaging the semblance measurements from sub-divided smaller windows. We also contaminated the synthetic data with white Gaussian noise, having different signal to noise ratios. Results showed no significant differences due to the fact that random arrivals in seismic data do not influence the semblance measurement. In next step to improve the diffraction imaging, we used tapered local semblance due to interference of diffractions with dominated reflection waves, other data and even other diffractions, especially at far offsets from diffraction’s apex. We called the proposed method as tapered local semblance method. The method weights the data from top to the bottom along the time axis, we also use less number of traces at shallow parts and more traces at deeper parts to reduce the harming effect of the interference. To coup with this task, we introduced a triangle taper to take few number of traces at the early arrival parts and more traces at the late arrival parts, instead of using a box with constant number of traces in the apertures from top to the bottom of the window. We tested several tapers with different angles of apex to determine the optimum one. We evaluated both methods on synthetic data as well as field recorded dataset. Both methods required no polarity reversal corrections to be applied. The obtained results showed the ability of our workflow to having higher resolution and good localization for diffractions from fault’s edges in synthetic data. The results obtained from using the tapered local semblance method on field recorded dataset showed more diffractivity than local semblance method.Diffractions carry useful and important information about subsurface features such as unconformities, faults, pinch-out, and so on. On the other hand, most of the information is encoded in diffractions. Polarity reversal across diffraction move out curves that are generated from fault’s edges, is a great challenge in seismic diffraction imaging. For the last few decades, several conventional methods in the pre-and post-stack domains, have been carried out for the diffractions characteristics and their locations. But most of these methods were not able to deal with the polarity reversal for diffraction imaging, some of them were time consuming, and needed to have some correction to deal with polarity changes, especially in diffraction caused by fault’s edges. Despite a large amount of research that has been carried out on diffraction imaging, very few studies have been devoted to the challenge of the polarity reversal across move out surfaces. We used the semblance function along with the hyperbolic move out curves for the diffractions that their travel times have been calculated using the double-square-root equation. As we know, using both semblance and Kirchhoff migration for diffraction imaging from fault’s edges, without which taking the polarity reversal into account would fail. This is caused by presence of same number of positive and negative wavelets in the diffraction move out curves. For solving this problem, we divided the global scanning window along hyperbolic move out surfaces into several subdivided window and the local semblance measurements over the sub-windows were performed separately. Every point in image domain is considered as a diffraction point that we call this points as image points. The final semblance measure at each image point is calculated by averaging the semblance measurements from sub-divided smaller windows. We also contaminated the synthetic data with white Gaussian noise, having different signal to noise ratios. Results showed no significant differences due to the fact that random arrivals in seismic data do not influence the semblance measurement. In next step to improve the diffraction imaging, we used tapered local semblance due to interference of diffractions with dominated reflection waves, other data and even other diffractions, especially at far offsets from diffraction’s apex. We called the proposed method as tapered local semblance method. The method weights the data from top to the bottom along the time axis, we also use less number of traces at shallow parts and more traces at deeper parts to reduce the harming effect of the interference. To coup with this task, we introduced a triangle taper to take few number of traces at the early arrival parts and more traces at the late arrival parts, instead of using a box with constant number of traces in the apertures from top to the bottom of the window. We tested several tapers with different angles of apex to determine the optimum one. We evaluated both methods on synthetic data as well as field recorded dataset. Both methods required no polarity reversal corrections to be applied. The obtained results showed the ability of our workflow to having higher resolution and good localization for diffractions from fault’s edges in synthetic data. The results obtained from using the tapered local semblance method on field recorded dataset showed more diffractivity than local semblance method.https://jesphys.ut.ac.ir/article_86909_8e3bb7cc1b45f051061287cdad9a67aa.pdfInstitute of Geophysics, University of TehranJournal of the Earth and Space Physics2538-371X48220220823Magnetic and IP/Res data inversion for investigation of the spatial relation between the geophysical models and mineralization in the southern Dalli Cu-Au porphyry depositMagnetic and IP/Res data inversion for investigation of the spatial relation between the geophysical models and mineralization in the southern Dalli Cu-Au porphyry deposit2612758550810.22059/jesphys.2022.329889.1007356FAMohammadHajheidariM.Sc. Student, Department of Mining Engineering, Isfahan University of Technology, Isfahan, Iran0000-0003-0920-7545Sayyed MohammadAbtahi ForooshaniAssistant Professor, Department of Mining Engineering, Isfahan University of Technology, Isfahan, Iran0000-0001-9133-357XHooshangAsadi HaroniAssistant Professor, Department of Mining Engineering, Isfahan University of Technology, Isfahan, Iran0000-0002-5627-3089KeytashMoshtaghianM.Sc. Student, Department of Mining Engineering, Isfahan University of Technology, Isfahan, Iran0000-0002-7499-762XGhazalJanghorbanM.Sc. Student, Department of Mining Engineering, Isfahan University of Technology, Isfahan, Iran0000-0003-2889-5080Journal Article20210905Because of declining high-grade ore deposits and increasing demands for metal resources, exploration of low-grade metal deposits, such as porphyries, have become feasible. Besides, humankind has spent most of the shallow metal ore deposits, and new prospecting projects focus on deeper deposits. Therefore, geophysical methods have gained more attention due to their ability to determine buried ore bodies' physical properties. Hence, most countries, including Iran, make significant investments in the geophysical exploration of deep porphyry deposits. According to widely accepted Lowell and Guilbert's model for porphyry copper deposits, the ore-bearing zones mainly concentrate at the edge of the potassic alteration zone. Pyrite, a highly conductive and chargeable metallic mineral, is a significant attribute in the potassic alteration. The model also states that the high susceptible magnetite-bearing rocks mainly occur at the bottom of the pyrite shell and the ore body. Due to the occurrence and presence of susceptible and conductive metallic minerals such as magnetite and pyrite in the potassic zone near to the ore body in the copper and gold porphyry deposits, the use of magnetometry, resistivity, and inducing polarization methods give reliable information about the location, depth, and shape of the deposits. For instance, in this research, we focus on the magnetic and IP/Res data in the southern Dalli porphyry deposit, with promising Cu-Au indices, which is located at Euromieh-Dokhtar ore-bearing zone Markazi Province. First, we applied standard processing techniques to remove the aliasing and regional effect in the magnetic data. Then, using the analytic signal technique, we showed the concentration of the magnetic sources over the study area. We also applied the power spectrum and Euler deconvolution techniques to the magnetic data and estimated the magnetic sources' depths. The estimated depth from the power spectrum is between the estimated depth from Euler deconvolution for possible sources with step and pillar shapes. Next, we used the average estimated depth from each of the depth estimation techniques in a three-dimensional magnetic data inversion as the depth of the sources in depth weighting. Also, we studied the inversion results via combining the cross-section of the magnetic susceptibility model along the boreholes and the lithology and geochemical information from core samples analysis. The results indicate that the higher grades for gold and copper occur at the edge of the magnetic sources and possible magnetite mineralization zones. The inversion results using the depth weighting with the depth extracted from the power spectrum show the best correlation and spatial relation with the geochemical data. Besides the magnetic data inversion, applying Oldenburg and Li algorithms for two-dimensional inverse modeling, we extracted the underground bodies' resistivity and chargeability model along with a IP/Res profile in the study area. The resulting chargeability models show a significant relationship with the presence of gold and copper mineralization. We also compared the resulting two-dimensional resistivity and changeability models with their corresponding magnetic susceptibility at the cross-sections along with the IP/Res. The comparison shows that the possible mineralization zones coincide with larger magnetic susceptibility values, high chargeability and low resistivity. The results show good accordance with Lowell and Guilbert's model. Also, highly susceptible rock in the shallower depth indicates that the erosion process has destroyed most possible orebody.Because of declining high-grade ore deposits and increasing demands for metal resources, exploration of low-grade metal deposits, such as porphyries, have become feasible. Besides, humankind has spent most of the shallow metal ore deposits, and new prospecting projects focus on deeper deposits. Therefore, geophysical methods have gained more attention due to their ability to determine buried ore bodies' physical properties. Hence, most countries, including Iran, make significant investments in the geophysical exploration of deep porphyry deposits. According to widely accepted Lowell and Guilbert's model for porphyry copper deposits, the ore-bearing zones mainly concentrate at the edge of the potassic alteration zone. Pyrite, a highly conductive and chargeable metallic mineral, is a significant attribute in the potassic alteration. The model also states that the high susceptible magnetite-bearing rocks mainly occur at the bottom of the pyrite shell and the ore body. Due to the occurrence and presence of susceptible and conductive metallic minerals such as magnetite and pyrite in the potassic zone near to the ore body in the copper and gold porphyry deposits, the use of magnetometry, resistivity, and inducing polarization methods give reliable information about the location, depth, and shape of the deposits. For instance, in this research, we focus on the magnetic and IP/Res data in the southern Dalli porphyry deposit, with promising Cu-Au indices, which is located at Euromieh-Dokhtar ore-bearing zone Markazi Province. First, we applied standard processing techniques to remove the aliasing and regional effect in the magnetic data. Then, using the analytic signal technique, we showed the concentration of the magnetic sources over the study area. We also applied the power spectrum and Euler deconvolution techniques to the magnetic data and estimated the magnetic sources' depths. The estimated depth from the power spectrum is between the estimated depth from Euler deconvolution for possible sources with step and pillar shapes. Next, we used the average estimated depth from each of the depth estimation techniques in a three-dimensional magnetic data inversion as the depth of the sources in depth weighting. Also, we studied the inversion results via combining the cross-section of the magnetic susceptibility model along the boreholes and the lithology and geochemical information from core samples analysis. The results indicate that the higher grades for gold and copper occur at the edge of the magnetic sources and possible magnetite mineralization zones. The inversion results using the depth weighting with the depth extracted from the power spectrum show the best correlation and spatial relation with the geochemical data. Besides the magnetic data inversion, applying Oldenburg and Li algorithms for two-dimensional inverse modeling, we extracted the underground bodies' resistivity and chargeability model along with a IP/Res profile in the study area. The resulting chargeability models show a significant relationship with the presence of gold and copper mineralization. We also compared the resulting two-dimensional resistivity and changeability models with their corresponding magnetic susceptibility at the cross-sections along with the IP/Res. The comparison shows that the possible mineralization zones coincide with larger magnetic susceptibility values, high chargeability and low resistivity. The results show good accordance with Lowell and Guilbert's model. Also, highly susceptible rock in the shallower depth indicates that the erosion process has destroyed most possible orebody.https://jesphys.ut.ac.ir/article_85508_69f706200da40de1981dccc8ecebd77c.pdfInstitute of Geophysics, University of TehranJournal of the Earth and Space Physics2538-371X48220220823Determination of 3D seismic wave velocity in Zagros collision zoneDetermination of 3D seismic wave velocity in Zagros collision zone2772928700410.22059/jesphys.2022.330058.1007358FAAmirTalebiPh.D. Graduated, Department of Seismology, Institute of Geophysics, University of Tehran, Tehran, Iran0000-0002-0588-4525HabibRahimiAssociate Professor, Department of Seismology, Institute of Geophysics, University of Tehran, Tehran, Iran0000-0002-2085-1043AliMoradiAssociate Professor, Department of Seismology, Institute of Geophysics, University of Tehran, Tehran, Iran0000-0002-0836-9027Journal Article20210926The Zagros orogenic belt was formed approximately 12 million years ago due to the convergence between the Arabian and Eurasian plates upon the closing of the Neo-Tethys Ocean. The Zagros is categorized as one of the youngest such settings on Earth, at an early stage of this collision. Many geophysical multiscale studies have been performed in the Zagros region based on different seismic and non-seismic data. Based on these studies, it can be concluded that the Zagros thrust belt has a crustal thickness of 45 ± 3 km, whereas beneath the Sanandaj-Sirjan zone, the Moho depth significantly increases up to 65 3 ± km. Among the many geophysical studies of Zagros and surrounding areas, local earthquake tomography (LET), which uses travel time data of both stations and earthquakes located in the study area, has never been performed for the entire Zagros. In this research, a 3D velocity model of body waves has been extracted using the information of the arrival time of 7783 earthquakes in the period of 2006 to 2018, recorded in the National Seismological Center and the broadband seismic network of Iran. The dataset used for tomography consists of 123,575 P- and 11,520 S-picks from 7783 events with magnitude greater than 2.5. We used the LOTOS code (Koulakov, 2009a) developed for simultaneous inversion for the 3D distributions of the P and S wave velocity anomalies and source locations. In the first step, LOTOS determines initial source locations using tabulated values of travel times previously calculated in a 1-D velocity model. The iterative algorithm of tomographic inversion includes the following steps: (1) Source relocations in the updated 3-D velocity structure based on the ray tracing bending method, (2) calculation of the first derivative matrix and (3) simultaneous inversion for P and S wave velocity anomalies, earthquake source parameters (4 parameters for each source), and station corrections. The inversion uses the LSQR method39. The distribution of estimated 3D velocity models correlates well with tectonic and geological conditions. The Vp and Vs anomalies, which are obtained independently, appear to be almost identical in the crust (depths smaller than 45 km). According to the results, the low velocity anomaly observed in the obtained models in the upper crust can be interpreted due to the presence of Cambrian-Miocene sediments with a thickness of at least 10 km that are spread throughout the Zagros. According to the obtained velocity models in the vertical sections, the Moho depth in the Sanandaj-Sirjan area increases significantly compared to the Zagros region. This increase in Moho depth is related to the subduction of the Arabic plate below the micro-continent of Central Iran, which increases the thickness of the crust (double crust) in the Sanandaj-Sirjan region. Using LOTOS code, the optimal one-dimensional velocity model for the whole Zagros collision zone is also presented. In this model, we can distinguish a ∼10 km thick sedimentary (V<sub>p </sub>∼4.90 km s<sup>-1</sup>), the upper crust down to ∼30 km (V<sub>p </sub>∼ 5.54 km s<sup>-1</sup>) and the lower crust down to ∼45 km (V<sub>p </sub>∼6.30 km s<sup>-1</sup>).The Zagros orogenic belt was formed approximately 12 million years ago due to the convergence between the Arabian and Eurasian plates upon the closing of the Neo-Tethys Ocean. The Zagros is categorized as one of the youngest such settings on Earth, at an early stage of this collision. Many geophysical multiscale studies have been performed in the Zagros region based on different seismic and non-seismic data. Based on these studies, it can be concluded that the Zagros thrust belt has a crustal thickness of 45 ± 3 km, whereas beneath the Sanandaj-Sirjan zone, the Moho depth significantly increases up to 65 3 ± km. Among the many geophysical studies of Zagros and surrounding areas, local earthquake tomography (LET), which uses travel time data of both stations and earthquakes located in the study area, has never been performed for the entire Zagros. In this research, a 3D velocity model of body waves has been extracted using the information of the arrival time of 7783 earthquakes in the period of 2006 to 2018, recorded in the National Seismological Center and the broadband seismic network of Iran. The dataset used for tomography consists of 123,575 P- and 11,520 S-picks from 7783 events with magnitude greater than 2.5. We used the LOTOS code (Koulakov, 2009a) developed for simultaneous inversion for the 3D distributions of the P and S wave velocity anomalies and source locations. In the first step, LOTOS determines initial source locations using tabulated values of travel times previously calculated in a 1-D velocity model. The iterative algorithm of tomographic inversion includes the following steps: (1) Source relocations in the updated 3-D velocity structure based on the ray tracing bending method, (2) calculation of the first derivative matrix and (3) simultaneous inversion for P and S wave velocity anomalies, earthquake source parameters (4 parameters for each source), and station corrections. The inversion uses the LSQR method39. The distribution of estimated 3D velocity models correlates well with tectonic and geological conditions. The Vp and Vs anomalies, which are obtained independently, appear to be almost identical in the crust (depths smaller than 45 km). According to the results, the low velocity anomaly observed in the obtained models in the upper crust can be interpreted due to the presence of Cambrian-Miocene sediments with a thickness of at least 10 km that are spread throughout the Zagros. According to the obtained velocity models in the vertical sections, the Moho depth in the Sanandaj-Sirjan area increases significantly compared to the Zagros region. This increase in Moho depth is related to the subduction of the Arabic plate below the micro-continent of Central Iran, which increases the thickness of the crust (double crust) in the Sanandaj-Sirjan region. Using LOTOS code, the optimal one-dimensional velocity model for the whole Zagros collision zone is also presented. In this model, we can distinguish a ∼10 km thick sedimentary (V<sub>p </sub>∼4.90 km s<sup>-1</sup>), the upper crust down to ∼30 km (V<sub>p </sub>∼ 5.54 km s<sup>-1</sup>) and the lower crust down to ∼45 km (V<sub>p </sub>∼6.30 km s<sup>-1</sup>).https://jesphys.ut.ac.ir/article_87004_b12667bb305a06f4092fb9b7375fcdac.pdfInstitute of Geophysics, University of TehranJournal of the Earth and Space Physics2538-371X48220220823Paleostrees analysis and Evaluation of Movement potential of Dochah Fault, Central IranPaleostrees analysis and Evaluation of Movement potential of Dochah Fault, Central Iran2933088544110.22059/jesphys.2022.327965.1007359FAMohadeseAjamiM.Sc. Graduated, School of Geology, College of Science, University of Tehran, Tehran, Iran0000-0003-2508-5728RezaNozaemAssistant Professor, School of Geology, College of Science, University of Tehran, Tehran, Iran0000-0002-2191-6988MohsenEliassiAssociate Professor, School of Geology, College of Science, University of Tehran, Tehran, Iran0000-0002-1913-5226SaeedMadanipourAssistant Professor, Department of Geology, Faculty of Science, Tarbiat Modares University, Tehran, Iran0000-0003-2175-9780SaeedHaj AminiInstructor, School of Geology, College of Science, University of Tehran, Tehran, Iran0000-0001-7274-4760VahidTavakoliAssociate Professor, School of Geology, College of Science, University of Tehran, Tehran, Iran0000-0002-4182-9259KosarShadramM.Sc. Student, Department of Geology, Faculty of Science, Tarbiat Modares University, Tehran, Iran0000-0002-1463-5571Journal Article20210929Qom region is one the significant area insight of geological features in Central Iran. Several researches have studied about the Cenozoic strata in terms of sedimentology, Stratigraphy and paleontology but, few structural detail data are available from this area. The most important exposure of the rock unites at the west of the Qom city is related to the Eocene volcanics, Lower Red, Qom and Upper Red Formations. Major structures at this area are Kamar Kuh and Mil anticlines, Yazdan syncline, Dochah and Sefid Kuh faults. Dochah Fault with E-W trending and ~70° dipping to the northward placed at the northwest termination of Qom-Zefreh Fault as a recent sinistral strike slip fault. This fault with ~15 km length separate Mil anticline from Yazdan syncline and eliminates the southern limb of Dochah overturned anticline. In this study, we focused on the Dochah Fault damaged zone in order to paleostress analysis using geometric and kinematic characteristics of fault slip data, which is achieved from the deformed Qom and Upper Red Formations. For this purpose, 100 fault slip data with precise and accurate geometric and kinematic characteristics have been measured in the field and analyzed with software Dasiy and Rotax methods. In order to determine the sense of shearing of the faults, the criteria of Petit (1987) and Doblas (1998) have been used. While the trend of the major structures is east-west but, most of slip data is related the transverse oblique slip faults, because the Dochah Fault passes through the soft materials of Lower Red Formation and consequently it is not possible or too hard to find the slicken line. Meanwhile, our results indicate the magnitudes of the axes of the maximum and minimum principal stress (σ1, σ3) as 030/05 and 285/05, Geometric and kinetic structural analysis related to the dochah fault and according to the spatial arrangement of the main stress axes indicate the readiness of the left-hand section on the right-hand section, especially in the western parts of the region (Caspian) attributed. oblate shape of field stress ellipsoid shape (R~0.7). Based on the field stress ellipsoid shape and the rotation of the fault data regarding the Anderson's theory for the compressive stress regime, the stress transition trajectory map has been prepared. The arrangement of maximum stress trajectories is consistent with the general stress regime in the Iranian crust and is consistent with the activity of the Dochah Fault. Different criteria have been proposed to evaluate the activity of a fault in terms of seismicity. In experimental studies, there are various estimates for selection of the part of the fault that the movement rediscovers for each tectonic seismic zone. Here, the possibility of moving Dochah Fault has been estimated by the method of Lee et al. (1997). In this method, the angular relationship between maximum principal stress axis (σ1) and the pole of the fault plane considered in order to evaluate the Fault Movement Potential (FMP) based on equation “FMP=f (G, σ)”. The angle between maximum principal stress axis (σ1) and the pole of Dochah Fault (θ) is equal to ~40° and so FMP=0.33 based on equation FMP= (θ-30°) ⁄ (30°) if θ∈[30°,60°]. This value of FMP indicates the low seismic potential of Dochah fault for movement and creating earthquakes.Qom region is one the significant area insight of geological features in Central Iran. Several researches have studied about the Cenozoic strata in terms of sedimentology, Stratigraphy and paleontology but, few structural detail data are available from this area. The most important exposure of the rock unites at the west of the Qom city is related to the Eocene volcanics, Lower Red, Qom and Upper Red Formations. Major structures at this area are Kamar Kuh and Mil anticlines, Yazdan syncline, Dochah and Sefid Kuh faults. Dochah Fault with E-W trending and ~70° dipping to the northward placed at the northwest termination of Qom-Zefreh Fault as a recent sinistral strike slip fault. This fault with ~15 km length separate Mil anticline from Yazdan syncline and eliminates the southern limb of Dochah overturned anticline. In this study, we focused on the Dochah Fault damaged zone in order to paleostress analysis using geometric and kinematic characteristics of fault slip data, which is achieved from the deformed Qom and Upper Red Formations. For this purpose, 100 fault slip data with precise and accurate geometric and kinematic characteristics have been measured in the field and analyzed with software Dasiy and Rotax methods. In order to determine the sense of shearing of the faults, the criteria of Petit (1987) and Doblas (1998) have been used. While the trend of the major structures is east-west but, most of slip data is related the transverse oblique slip faults, because the Dochah Fault passes through the soft materials of Lower Red Formation and consequently it is not possible or too hard to find the slicken line. Meanwhile, our results indicate the magnitudes of the axes of the maximum and minimum principal stress (σ1, σ3) as 030/05 and 285/05, Geometric and kinetic structural analysis related to the dochah fault and according to the spatial arrangement of the main stress axes indicate the readiness of the left-hand section on the right-hand section, especially in the western parts of the region (Caspian) attributed. oblate shape of field stress ellipsoid shape (R~0.7). Based on the field stress ellipsoid shape and the rotation of the fault data regarding the Anderson's theory for the compressive stress regime, the stress transition trajectory map has been prepared. The arrangement of maximum stress trajectories is consistent with the general stress regime in the Iranian crust and is consistent with the activity of the Dochah Fault. Different criteria have been proposed to evaluate the activity of a fault in terms of seismicity. In experimental studies, there are various estimates for selection of the part of the fault that the movement rediscovers for each tectonic seismic zone. Here, the possibility of moving Dochah Fault has been estimated by the method of Lee et al. (1997). In this method, the angular relationship between maximum principal stress axis (σ1) and the pole of the fault plane considered in order to evaluate the Fault Movement Potential (FMP) based on equation “FMP=f (G, σ)”. The angle between maximum principal stress axis (σ1) and the pole of Dochah Fault (θ) is equal to ~40° and so FMP=0.33 based on equation FMP= (θ-30°) ⁄ (30°) if θ∈[30°,60°]. This value of FMP indicates the low seismic potential of Dochah fault for movement and creating earthquakes.https://jesphys.ut.ac.ir/article_85441_32937f18e71b21bec524cb08b13a43fd.pdfInstitute of Geophysics, University of TehranJournal of the Earth and Space Physics2538-371X48220220823Analysis and prediction of EOP time series using LSHE+ARMA methodAnalysis and prediction of EOP time series using LSHE+ARMA method3093238545110.22059/jesphys.2022.333545.1007380FAMohammad AliSharifiAssociate Professor, Department of Surveying and Geomatics Engineering, Faculty of Engineering, University of Tehran, Tehran, Iran0000-0003-0745-4147ShayanShirafkan.Sc. Student, Department of Surveying and Geomatics Engineering, Faculty of Engineering, University of Tehran, Tehran, Iran0000-0002-2098-3900Seyyed MohsenKhazraeiPh.D. Graduated, Department of Geomatics Engineering, Faculty of Civil Engineering and Transportation, University of Isfahan, Iran0000-0003-3543-3722Ali RezaAmiri SimkoeiProfessor, Department of Geomatics Engineering, Faculty of Civil Engineering and Transportation, University of Isfahan, Iran0000-0002-2952-0160Journal Article20211115The rotation of the solid Earth with respect to inertial space is not constant due to the changes of external gravitational forces and internal dynamics. Earth orientation parameters (EOP), including, the Earth’s polar motion (PM), Anomalies in the Earth’s angular velocity and celestial pole offsets (CPO), describe these irregularities in the Earth’s rotation. Anomalies in the axis defined by the celestial intermediate pole (CIP) with respect to the Z axis of the terrestrial reference system are named as PM. The CPO are expressed as the deviations, dX and dY, between the observed CIP and the conventional CIP position. The difference between the smoothed principal form of universal time UT1 and the coordinated universal time UTC denotes the Earth’s rotation angle, which together with the xp, yp terrestrial pole coordinates, forms a set of Earth orientation parameters (EOP). In addition to the other EOP, the length of day (LOD) is used to model the Anomalies in the Earth’s rotation rate. LOD is the difference between the duration of the day measured by space geodesy and nominal day of 86,400 s duration.<br />Generally, EOP are the parameters that provide the rotation of the International Terrestrial Reference System (ITRS) to the International Celestial Reference System (ICRS) as a function of time. However, the EOP are computed using the modern space geodetic techniques such as Doppler Orbitography and Radiopositioning Integrated by Satellite (DORIS), Satellite Laser Ranging (SLR), Very Long Baseline Interferometry (VLBI) and the Global Navigation Satellite System (GNSS), they are unavailable to the real-time applications due to the data processing complexities. Accurate and rapid EOP predictions are required for different fields like precise orbit determinations of artificial Earth satellites, positional astronomy, space navigation and geophysical phenomena.<br />There are many different methods for analysis and prediction of EOP time series including deep learning methods, least square (LS) with autoregressive (AR) and also Singular Spectrum Analysis as a non-parametric method.<br />In this research Least Square Harmonic Estimation analysis is used to investigate the frequencies of EOP. First, the solid and ocean tide terms are modeled based on IERS technical notes. These effects are removed from LOD time series. The remained time series are named as LODR time series. The univariate time series analysis is then applied to the LODR time series and multivariate analysis is used for detecting the PM periodic patterns. Applying these methods to the 40 years of observations of EOP (since 1 January 1980 to 31 December 2020) revealed the Chandler, annual, semi Chandler, semi-annual and annual signals as the main periodic signals in the EOP time series. The functional model is then formed using all detected signals in order to model the deterministic variations of EOP time series.<br />In order to model the remained non-deterministic variations, an ARMA (Autoregressive Moving Average) model is fitted to the least square residuals. The Akaike's Information Criterion (AIC) is used to investigate the optimized order of ARMA model.<br />The EOP is then predicted for the first 20 days of 2021, using the pre-identified functional model for the deterministic part and the ARMA model for the non-deterministic part of the time series variations. For the prediction of LOD time series, after creating the functional model of LODR time series, the solid and ocean tide terms are added to the functional model of LODR.<br />Finally, in order to validate the accuracy of the proposed method, a comparison is made with an EOP prediction study that used the ANN (Artificial Neural Network) and ANFIS (Adaptive Network Based Fuzzy Inference System) methods for short term prediction of EOP.<br />The result shows that the accuracy of the proposed method is better than the previous study and the method can be used for accurate prediction of EOP time series.The rotation of the solid Earth with respect to inertial space is not constant due to the changes of external gravitational forces and internal dynamics. Earth orientation parameters (EOP), including, the Earth’s polar motion (PM), Anomalies in the Earth’s angular velocity and celestial pole offsets (CPO), describe these irregularities in the Earth’s rotation. Anomalies in the axis defined by the celestial intermediate pole (CIP) with respect to the Z axis of the terrestrial reference system are named as PM. The CPO are expressed as the deviations, dX and dY, between the observed CIP and the conventional CIP position. The difference between the smoothed principal form of universal time UT1 and the coordinated universal time UTC denotes the Earth’s rotation angle, which together with the xp, yp terrestrial pole coordinates, forms a set of Earth orientation parameters (EOP). In addition to the other EOP, the length of day (LOD) is used to model the Anomalies in the Earth’s rotation rate. LOD is the difference between the duration of the day measured by space geodesy and nominal day of 86,400 s duration.<br />Generally, EOP are the parameters that provide the rotation of the International Terrestrial Reference System (ITRS) to the International Celestial Reference System (ICRS) as a function of time. However, the EOP are computed using the modern space geodetic techniques such as Doppler Orbitography and Radiopositioning Integrated by Satellite (DORIS), Satellite Laser Ranging (SLR), Very Long Baseline Interferometry (VLBI) and the Global Navigation Satellite System (GNSS), they are unavailable to the real-time applications due to the data processing complexities. Accurate and rapid EOP predictions are required for different fields like precise orbit determinations of artificial Earth satellites, positional astronomy, space navigation and geophysical phenomena.<br />There are many different methods for analysis and prediction of EOP time series including deep learning methods, least square (LS) with autoregressive (AR) and also Singular Spectrum Analysis as a non-parametric method.<br />In this research Least Square Harmonic Estimation analysis is used to investigate the frequencies of EOP. First, the solid and ocean tide terms are modeled based on IERS technical notes. These effects are removed from LOD time series. The remained time series are named as LODR time series. The univariate time series analysis is then applied to the LODR time series and multivariate analysis is used for detecting the PM periodic patterns. Applying these methods to the 40 years of observations of EOP (since 1 January 1980 to 31 December 2020) revealed the Chandler, annual, semi Chandler, semi-annual and annual signals as the main periodic signals in the EOP time series. The functional model is then formed using all detected signals in order to model the deterministic variations of EOP time series.<br />In order to model the remained non-deterministic variations, an ARMA (Autoregressive Moving Average) model is fitted to the least square residuals. The Akaike's Information Criterion (AIC) is used to investigate the optimized order of ARMA model.<br />The EOP is then predicted for the first 20 days of 2021, using the pre-identified functional model for the deterministic part and the ARMA model for the non-deterministic part of the time series variations. For the prediction of LOD time series, after creating the functional model of LODR time series, the solid and ocean tide terms are added to the functional model of LODR.<br />Finally, in order to validate the accuracy of the proposed method, a comparison is made with an EOP prediction study that used the ANN (Artificial Neural Network) and ANFIS (Adaptive Network Based Fuzzy Inference System) methods for short term prediction of EOP.<br />The result shows that the accuracy of the proposed method is better than the previous study and the method can be used for accurate prediction of EOP time series.https://jesphys.ut.ac.ir/article_85451_73a8064254f6503bd5a6f6f226107733.pdfInstitute of Geophysics, University of TehranJournal of the Earth and Space Physics2538-371X48220220823Analysis of behavioral pattern the basic parameters in foreshocks with target for the prediction of big earthquakes in IranAnalysis of behavioral pattern the basic parameters in foreshocks with target for the prediction of big earthquakes in Iran3253458692110.22059/jesphys.2022.333795.1007382FAAliSaketPh.D. Graduated, Faculty of Geoscience, Kharazmi University, Tehran, Iran0000-0003-3129-4761Seyed MahmoudFatemi AghdaProfessor, Faculty of Geoscience, Kharazmi University, Tehran, Iran0000-0002-4368-273xHosseinSadeghiAssociate Professor, Department of Geology, Faculty of Science, Ferdowsi University Mashhad, Mashhad, Iran0000-0002-7283-3082AhmadFahimifarProfessor, Faculty of civil engineering and environment, Amirkabir University of Technology, Tehran, Iran0000-0002-9397-6269Journal Article20211219The analysis of the basic parameters of the foreshocks is one of the most applied researches for risk reduction of earthquakes. Because identification of behavioral pattern of foreshocks can help researchers in detection of the active fault conditions in different areas. Also accurate analysis of these parameters help to study of earthquake prediction as more effective. In this study, we study about behavioral pattern of foreshocks in different tectonic zons in Iran. This research was conducted for prediction of probability the earthquakes with M>5 in Iran. According to this research, accurate analysis of the basic seismic parameters of foreshocks (including: relationship between depth and magnitude of foreshocks) is studied with target for the prediction of big earthquake in various zons for a ten-year period (from 2007 until 2017). The results of this research suggest that there are certain similarities in the magnitude-depth models for the one zone and also different for various zones. Therefore, this can be used as a precursor in earthquake prediction with Magnitude>5 for different zones in Iran. The important results presented in this article can be presented in the following cases:- Investigation of the information of seismicity parameters of foreshocks regarding the relationship between the focal depth of the main earthquake and the frequency of the foreshocks that used in some parts of the world as a precursor of earthquake suggested that main shocks with M>5 and shallow depth have more foreshocks abundances (Fig 2). - Due to the relationship between the type of fault with the occurrence and non-occurrence of aftershocks in different parts of the world, in the case of earthquakes greater than 5 in Iran, in earthquakes with reverse faults have relatively more aftershocks recorded compared to strike-slip faults.- The results of the statistical study conducted in this study show that for earthquakes with reverse fault, the frequency of foreshocks increases with magnitude. However, we do not see such conditions for earthquakes with faults of strike-slip.- The result of this study shows that more earthquake especially in Zagros zone and near salt domes happened without foreshocks. The reason for this is related to effect of salt dome on movement fault from slide to creep. The creep is a gradual movement and it is not usually accompanied by rapid movement such as slides that lead to large and recordable earthquakes.- Based on the present study on earthquakes, for the Zagros (especially in the northern and central part) and Central Iran and Sanandaj-Sirjan, can be used more confidently as a precursor of earthquake because in this zones earthquakes happened with more foreshocks.- In Zagros and Iran Markazi zone the relationship between variations of the depth and magnitude of foreshocks is fruitful for predicting of the main shocks.- For other zones we need to have more complete data bank that has earthquakes with higher frequency of foreshocks. Based on this data bank we can present suitable relations and models for the study of foreshock with the aim of predicting the big earthquakes.The analysis of the basic parameters of the foreshocks is one of the most applied researches for risk reduction of earthquakes. Because identification of behavioral pattern of foreshocks can help researchers in detection of the active fault conditions in different areas. Also accurate analysis of these parameters help to study of earthquake prediction as more effective. In this study, we study about behavioral pattern of foreshocks in different tectonic zons in Iran. This research was conducted for prediction of probability the earthquakes with M>5 in Iran. According to this research, accurate analysis of the basic seismic parameters of foreshocks (including: relationship between depth and magnitude of foreshocks) is studied with target for the prediction of big earthquake in various zons for a ten-year period (from 2007 until 2017). The results of this research suggest that there are certain similarities in the magnitude-depth models for the one zone and also different for various zones. Therefore, this can be used as a precursor in earthquake prediction with Magnitude>5 for different zones in Iran. The important results presented in this article can be presented in the following cases:- Investigation of the information of seismicity parameters of foreshocks regarding the relationship between the focal depth of the main earthquake and the frequency of the foreshocks that used in some parts of the world as a precursor of earthquake suggested that main shocks with M>5 and shallow depth have more foreshocks abundances (Fig 2). - Due to the relationship between the type of fault with the occurrence and non-occurrence of aftershocks in different parts of the world, in the case of earthquakes greater than 5 in Iran, in earthquakes with reverse faults have relatively more aftershocks recorded compared to strike-slip faults.- The results of the statistical study conducted in this study show that for earthquakes with reverse fault, the frequency of foreshocks increases with magnitude. However, we do not see such conditions for earthquakes with faults of strike-slip.- The result of this study shows that more earthquake especially in Zagros zone and near salt domes happened without foreshocks. The reason for this is related to effect of salt dome on movement fault from slide to creep. The creep is a gradual movement and it is not usually accompanied by rapid movement such as slides that lead to large and recordable earthquakes.- Based on the present study on earthquakes, for the Zagros (especially in the northern and central part) and Central Iran and Sanandaj-Sirjan, can be used more confidently as a precursor of earthquake because in this zones earthquakes happened with more foreshocks.- In Zagros and Iran Markazi zone the relationship between variations of the depth and magnitude of foreshocks is fruitful for predicting of the main shocks.- For other zones we need to have more complete data bank that has earthquakes with higher frequency of foreshocks. Based on this data bank we can present suitable relations and models for the study of foreshock with the aim of predicting the big earthquakes.https://jesphys.ut.ac.ir/article_86921_0bace7b448d342c3b637a7e12eb979ea.pdfInstitute of Geophysics, University of TehranJournal of the Earth and Space Physics2538-371X48220220823Observing of Pre-flare Very Long-period Pulsations, for 12 Solar Flares, as a Sign of Flare’s OnsetObserving of Pre-flare Very Long-period Pulsations, for 12 Solar Flares, as a Sign of Flare’s Onset3473598690810.22059/jesphys.2022.328165.1007344FAMaliheJalali RadM.Sc. Student, Department of Physics, Payame Noor University, Tehran, Iran0000-0002-0815-4974NargesFathalianAssistant Professor, Department of Physics, Payame Noor University, Tehran, Iran (P.O.Box 19395-3697)0000-0003-2824-1773Journal Article20210830Solar flares are sudden bursts in the solar atmosphere, which have emissions, from radio wavelengths up to gamma rays, and according to their energy are classified into different classes (A, B, C, M, and X, respectively). The process of releasing magnetic energy in flares is done by magnetic reconnection, which is often created by a complex magnetic field. Flares accelerate many electrons and ions, raising their energy to the limit of relative energy. These accelerating particles play a very important role in the release of large solar flare energies. Considering the fact that flares emit radiation when they explode, most of them create light spectrum and sometimes X-rays and ultraviolet rays, which are emitted mainly by the photosphere and chromosphere into concentrated sources called footpoints and ribbons. These radiations and emissions occur when the lower layers of the sun's atmosphere heat up during a flare, and this heating due to the collision of particles probably plays an important role in the occurrence of the flare. In addition, they emit high-energy radiation such as hard X-rays (HXR) from electrons and gamma rays from ions. The main part of these emissions is in the form of electromagnetic emission (soft X-rays) and energetic particles. Emissions radiated from a large flare or a solar mass eruption (with an energy more than J), when reaching the earth, can have destructive effects on the Earth's atmosphere, as well as the orbits of satellites or magnetic and electrical facilities of devices like ships and airplanes. Therefore, predicting the time of the flare occurrence and determining its class type can help reduce these destructive effects.<br />One of the observable structures that can be seen before a flare occurs, are oscillations with very long period pulsations (VLPs) of the order of 8-30 minutes, which occur about one to two hours before the flare onset, and were first reported by Tan et al. (2016) in the pre-flare phase. MHD oscillations and longitudinal electric current in flare loops can be appropriate candidates to explain the formation of VLPs. Investigating pre-flare VLPs can also help us in understanding the origin of flares. With the help of observational data of X-ray radiation (SXR), onboard the GOES satellite, during the pre-flare phase, these pulses can be observed at similar time scales during flare processes.<br />In this paper, using the abovementioned data, we selected eighteen flares for the study of which 6 flares are in class C and 12 flares are in class M. Of these, twelve had typical VLPs before flare-onset, which were all in the M class, with the exception of one. The periodicity that we calculated for the VLPs of these flares, with the help of the Fast Fourier Transform is 14 to 28.9 minutes, which is in agreement with the results of Tan et al. (2016). The number of pulses observed in each pre-flare is between 3 and 7. For the other six remaining flares of our selection, no typical pre-flare VLP was observed, which all but one of them, were in class C.Solar flares are sudden bursts in the solar atmosphere, which have emissions, from radio wavelengths up to gamma rays, and according to their energy are classified into different classes (A, B, C, M, and X, respectively). The process of releasing magnetic energy in flares is done by magnetic reconnection, which is often created by a complex magnetic field. Flares accelerate many electrons and ions, raising their energy to the limit of relative energy. These accelerating particles play a very important role in the release of large solar flare energies. Considering the fact that flares emit radiation when they explode, most of them create light spectrum and sometimes X-rays and ultraviolet rays, which are emitted mainly by the photosphere and chromosphere into concentrated sources called footpoints and ribbons. These radiations and emissions occur when the lower layers of the sun's atmosphere heat up during a flare, and this heating due to the collision of particles probably plays an important role in the occurrence of the flare. In addition, they emit high-energy radiation such as hard X-rays (HXR) from electrons and gamma rays from ions. The main part of these emissions is in the form of electromagnetic emission (soft X-rays) and energetic particles. Emissions radiated from a large flare or a solar mass eruption (with an energy more than J), when reaching the earth, can have destructive effects on the Earth's atmosphere, as well as the orbits of satellites or magnetic and electrical facilities of devices like ships and airplanes. Therefore, predicting the time of the flare occurrence and determining its class type can help reduce these destructive effects.<br />One of the observable structures that can be seen before a flare occurs, are oscillations with very long period pulsations (VLPs) of the order of 8-30 minutes, which occur about one to two hours before the flare onset, and were first reported by Tan et al. (2016) in the pre-flare phase. MHD oscillations and longitudinal electric current in flare loops can be appropriate candidates to explain the formation of VLPs. Investigating pre-flare VLPs can also help us in understanding the origin of flares. With the help of observational data of X-ray radiation (SXR), onboard the GOES satellite, during the pre-flare phase, these pulses can be observed at similar time scales during flare processes.<br />In this paper, using the abovementioned data, we selected eighteen flares for the study of which 6 flares are in class C and 12 flares are in class M. Of these, twelve had typical VLPs before flare-onset, which were all in the M class, with the exception of one. The periodicity that we calculated for the VLPs of these flares, with the help of the Fast Fourier Transform is 14 to 28.9 minutes, which is in agreement with the results of Tan et al. (2016). The number of pulses observed in each pre-flare is between 3 and 7. For the other six remaining flares of our selection, no typical pre-flare VLP was observed, which all but one of them, were in class C.https://jesphys.ut.ac.ir/article_86908_dc346f33f08b9935042f5bb3d394ef09.pdfInstitute of Geophysics, University of TehranJournal of the Earth and Space Physics2538-371X48220220823On the design and implementation of digital filters to process meteorological signalsOn the design and implementation of digital filters to process meteorological signals3613808546010.22059/jesphys.2022.329174.1007349FAAbolfazlNeyestaniAssistant Professor, Department of Physics, Razi University, Kermanshah, Iran0000-0002-0771-1548Journal Article20210826Separation of different frequency bands in the complex and combined signals related to meteorological variables and also climatic indices requires the use of digital filtering methods. In this way, the information on different frequency bands can be organized and used. Given that these signals generally exhibit complex and nonlinear behavior, the use of mathematical filtering methods to identify their stochastic and periodic components leads to a better understanding of their behavior and helps modeling them as well. Therefore, the use of digital filters in order to recognize regular variabilities and facilitate statistical forecasting is one of the main goals in this field.<br />The design and implementation of these filters are possible both in time and frequency domains. In the frequency space, this process is performed based on the Fourier transform of the signals on the basis of the Fast Fourier Transform Algorithm (FFT), in which the variances of the desired signal can be extracted based on spectral analysis in different frequencies. By employing different types of non-recursive and recursive digital filters, which they can be implemented as low-pass, high-pass, band-pass, and band-stop, the related signal inthe time domain for each state can be constructed, and the corresponding spectrum can be studied. The isolated spectrum can be related to the effect of a special phenomenon that influences the main signal. In addition, it is possible to remove high frequency components from the original signal, which include noises and may not contain important information. Moreover, the process of optimal smoothing the original signal can also be carried out.<br />In this study, different digital filters have been designed and then applied to meteorological data such as monthly surface temperature and precipitation. Two synoptic stations over Iran are selected and the related discrete monthly signals are constructed for 504 months during 1979-2021. Then, the moving average (MA) filter is used as a main filter, because it is the most common filter in digital signal processing (DSP), and also it is the easiest digital filter to understand and use. In spite of its simplicity, the moving average filter is optimal<em> </em>for a common task such as reducing random noise while retaining a sharp step response. This makes it the premier filter for time domain encoded signals. The filtering process in this study is conducted to denoise the original signals, and also to examine seasonal, annual, and inter-annual components of the original signals. Since employed filters are digital, they must be applied to the initial discrete signal in the form of convolution with the finite impulse response (FIR) of the filter in the time domain, or they can be applied in the form of multiplication in the frequency domain based on discrete Fourier transform and then using of the inverse Fourier transform to recover the desired signal.<br />The results of this study show the importance of using digital filters in analyzing the spectral contents of meteorological signals. Furthermore, the Hamming filter, which is defined based on the cosine truncation and windowing, shows better performance in attenuating Gibbs oscillations in the lateral sidelobes of the filter frequency response than the simple moving average (MA) filter. In addition, the correlation analysis is carried out separately to indicate the linear relationships between different frequency components of the signals. The higher correlations are observed in annual frequency bands of the temperature and precipitation signals for the selected stations. It shows the effect of external climate forcing on both temperature and precipitation that is stemmed from the earth’s motion around the sun during a year. Obviously, choosing more weights in the design of a filter can improve the filtering performance, but it should be avoided to use more weights than necessary.Separation of different frequency bands in the complex and combined signals related to meteorological variables and also climatic indices requires the use of digital filtering methods. In this way, the information on different frequency bands can be organized and used. Given that these signals generally exhibit complex and nonlinear behavior, the use of mathematical filtering methods to identify their stochastic and periodic components leads to a better understanding of their behavior and helps modeling them as well. Therefore, the use of digital filters in order to recognize regular variabilities and facilitate statistical forecasting is one of the main goals in this field.<br />The design and implementation of these filters are possible both in time and frequency domains. In the frequency space, this process is performed based on the Fourier transform of the signals on the basis of the Fast Fourier Transform Algorithm (FFT), in which the variances of the desired signal can be extracted based on spectral analysis in different frequencies. By employing different types of non-recursive and recursive digital filters, which they can be implemented as low-pass, high-pass, band-pass, and band-stop, the related signal inthe time domain for each state can be constructed, and the corresponding spectrum can be studied. The isolated spectrum can be related to the effect of a special phenomenon that influences the main signal. In addition, it is possible to remove high frequency components from the original signal, which include noises and may not contain important information. Moreover, the process of optimal smoothing the original signal can also be carried out.<br />In this study, different digital filters have been designed and then applied to meteorological data such as monthly surface temperature and precipitation. Two synoptic stations over Iran are selected and the related discrete monthly signals are constructed for 504 months during 1979-2021. Then, the moving average (MA) filter is used as a main filter, because it is the most common filter in digital signal processing (DSP), and also it is the easiest digital filter to understand and use. In spite of its simplicity, the moving average filter is optimal<em> </em>for a common task such as reducing random noise while retaining a sharp step response. This makes it the premier filter for time domain encoded signals. The filtering process in this study is conducted to denoise the original signals, and also to examine seasonal, annual, and inter-annual components of the original signals. Since employed filters are digital, they must be applied to the initial discrete signal in the form of convolution with the finite impulse response (FIR) of the filter in the time domain, or they can be applied in the form of multiplication in the frequency domain based on discrete Fourier transform and then using of the inverse Fourier transform to recover the desired signal.<br />The results of this study show the importance of using digital filters in analyzing the spectral contents of meteorological signals. Furthermore, the Hamming filter, which is defined based on the cosine truncation and windowing, shows better performance in attenuating Gibbs oscillations in the lateral sidelobes of the filter frequency response than the simple moving average (MA) filter. In addition, the correlation analysis is carried out separately to indicate the linear relationships between different frequency components of the signals. The higher correlations are observed in annual frequency bands of the temperature and precipitation signals for the selected stations. It shows the effect of external climate forcing on both temperature and precipitation that is stemmed from the earth’s motion around the sun during a year. Obviously, choosing more weights in the design of a filter can improve the filtering performance, but it should be avoided to use more weights than necessary.https://jesphys.ut.ac.ir/article_85460_1518e4744a7db7e801b583f0225c1142.pdfInstitute of Geophysics, University of TehranJournal of the Earth and Space Physics2538-371X48220220823A study of clear air turbulence by spontaneous imbalance theoryA study of clear air turbulence by spontaneous imbalance theory3814018545010.22059/jesphys.2022.329288.1007351FASaraMahmoodvandM.Sc. Graduated, Department of Space Physics, Institute of Geophysics, University of Tehran, Tehran, Iran0000-0002-4467-8353MohammadMirzaeiAssociate Professor, Department of Space Physics, Institute of Geophysics, University of Tehran, Tehran, Iran0000-0003-0813-3994Ali RezaMohebalhojehProfessor, Department of Space Physics, Institute of Geophysics, University of Tehran, Tehran, Iran0000-0002-5906-8486Journal Article20210826Emission of inertia–gravity waves (IGWs) through imbalance is a well-known cause of clear air turbulence (CAT) in the upper troposphere. IGWs may initiate CAT by locally modifying the environmental characteristics of the meteorological quantities like static stability and wind shear. CAT is a micro-scale phenomenon for which there are also mechanisms other than IGWs. Accurate forecasting methods using numerical models and CAT diagnostic indices are still being studied and developed (Sharman and Lane, 2016). Following Knox et al. (2008) (hereafter KMW), the current study is focused on detecting CAT by spontaneous imbalance theory and the effect of IGWs on the flow.<br />For this purpose, the lifecycle of the baroclinic waves, including their phases of growth, overturn and decay as well as the generation and propagation of IGWs are investigated by numerical simulation using the Weather Research and Forecasting (WRF) model in a channel of 4000 km length, 10000 km width and 22 km height in respectively the zonal, meridional and vertical directions on the f plane, with a horizontal resolution of 25 km and vertical resolution of 0.25 km. Based on the wave–vortex decomposition (WVD) method, the unbalanced flow, and the dimensional and non-dimensional IGW amplitude have been estimated. In the next step, the non-dimensional wave amplitude has been alternatively determined for reference, based on the Lighthill–Ford theory of spontaneous imbalance in KMW method. Then the turbulent kinetic energy (TKE) dissipation and eddy dissipation rate (EDR) have been calculated to determine the intensity and location of CAT.<br />The results showed that KMW method uses a proportionality constant to make the non-dimensional wave amplitude as order of the Rossby number and determines the constant empirically by matching distributions of pilot reports of turbulence to the pattern of TKE dissipation. For this reason, the EDR has the best fit with the location of observed CAT and the minimum value of Richardson number. This is while most values of the non-dimensional wave amplitudes calculated by the WVD and harmonic divergence analysis are less than unity and have values of the order of the Rossby number itself. On day 8, when the baroclinic wave and IGWs are at their peak of activity, the pattern of distribution of EDR by WVD indicates that there is moderate turbulence all around the jet stream region, and the maximum values of EDR are located below the jet core and in the jet-exit region, which is similar to the location of wave activity and location of CAT in previous studies. Also minimum values of Richardson number are at the jet-exit region where the maxima of EDR reveal moderate turbulence there. The distribution of EDR by KMW, unlike the distribution of EDR by WVD, shows that in most areas of the flow, there is no sign of turbulence except in a few patchy places near the jet region, where moderate turbulence is predicted. Thus making use of an optimal WDV could improve the accuracy of detecting unbalanced parts of the flow and predicting areas of CAT in the upper troposphere in the vicinity of the jet stream.Emission of inertia–gravity waves (IGWs) through imbalance is a well-known cause of clear air turbulence (CAT) in the upper troposphere. IGWs may initiate CAT by locally modifying the environmental characteristics of the meteorological quantities like static stability and wind shear. CAT is a micro-scale phenomenon for which there are also mechanisms other than IGWs. Accurate forecasting methods using numerical models and CAT diagnostic indices are still being studied and developed (Sharman and Lane, 2016). Following Knox et al. (2008) (hereafter KMW), the current study is focused on detecting CAT by spontaneous imbalance theory and the effect of IGWs on the flow.<br />For this purpose, the lifecycle of the baroclinic waves, including their phases of growth, overturn and decay as well as the generation and propagation of IGWs are investigated by numerical simulation using the Weather Research and Forecasting (WRF) model in a channel of 4000 km length, 10000 km width and 22 km height in respectively the zonal, meridional and vertical directions on the f plane, with a horizontal resolution of 25 km and vertical resolution of 0.25 km. Based on the wave–vortex decomposition (WVD) method, the unbalanced flow, and the dimensional and non-dimensional IGW amplitude have been estimated. In the next step, the non-dimensional wave amplitude has been alternatively determined for reference, based on the Lighthill–Ford theory of spontaneous imbalance in KMW method. Then the turbulent kinetic energy (TKE) dissipation and eddy dissipation rate (EDR) have been calculated to determine the intensity and location of CAT.<br />The results showed that KMW method uses a proportionality constant to make the non-dimensional wave amplitude as order of the Rossby number and determines the constant empirically by matching distributions of pilot reports of turbulence to the pattern of TKE dissipation. For this reason, the EDR has the best fit with the location of observed CAT and the minimum value of Richardson number. This is while most values of the non-dimensional wave amplitudes calculated by the WVD and harmonic divergence analysis are less than unity and have values of the order of the Rossby number itself. On day 8, when the baroclinic wave and IGWs are at their peak of activity, the pattern of distribution of EDR by WVD indicates that there is moderate turbulence all around the jet stream region, and the maximum values of EDR are located below the jet core and in the jet-exit region, which is similar to the location of wave activity and location of CAT in previous studies. Also minimum values of Richardson number are at the jet-exit region where the maxima of EDR reveal moderate turbulence there. The distribution of EDR by KMW, unlike the distribution of EDR by WVD, shows that in most areas of the flow, there is no sign of turbulence except in a few patchy places near the jet region, where moderate turbulence is predicted. Thus making use of an optimal WDV could improve the accuracy of detecting unbalanced parts of the flow and predicting areas of CAT in the upper troposphere in the vicinity of the jet stream.https://jesphys.ut.ac.ir/article_85450_7e1b701e996ab75b21d7781c86fceb90.pdfInstitute of Geophysics, University of TehranJournal of the Earth and Space Physics2538-371X48220220823Deterministic and Fuzzy Evaluation of Human and Climate Contributions in Changing Hydrologic Regime: A Case Study of the Gorganrood Watershed at Tamar River Hydrometric StationDeterministic and Fuzzy Evaluation of Human and Climate Contributions in Changing Hydrologic Regime: A Case Study of the Gorganrood Watershed at Tamar River Hydrometric Station4034208708110.22059/jesphys.2022.329382.1007352FAMohammad MasoudMohammadpour KhoieM.Sc. Student, School of Civil Engineering, College of Engineering, University of Tehran, Tehran, Iran0000-0001-5101-6713MohsenNasseriAssociate Professor, School of Civil Engineering, College of Engineering, University of Tehran, Tehran, Iran0000-0002-7584-7631Seyyed Mohammad AliBanihashemiAssociate Professor, School of Civil Engineering, College of Engineering, University of Tehran, Tehran, Iran0000-0001-5399-5329Journal Article20210828Human and climate are two major scio-hydrologic drivers that determine hydrological regimes and patterns. In this regards, Land Use and Land Cover (LULC) changes, agricultural development, etc., on global and regional scales, hydro-climatological components have influenced these regimes. The effects of each driver on the variation of hydrological components have been assessed in different studies, but these approaches are not accurate enough at watershed-scales that experience the simultaneous impacts of climate dynamics and LULC changes. Different studies have considered both climate and human altertions in the hydrological cycle, and quantified their contributions in such basin. Results of these researches can help decision makers in water management of the pros and cons of water and land use policies. The Gorganrood watershed is an important basin in the northern part of Iran, especially from the agricultural point of view, which has considerably experienced hydrological and extreme events changes. While the consequence of each climate change and LULC changes have been assessed in the watershed, there is no study, which considers the complicated interactions of these drivers. In this paper, the authors firstly evaluated the contributions of LULC and climate change on the variation of streamflow. Secondly, the modified fuzzy arithmetic method has been used to achieve their fuzzy contributions. To this purpose, the computational period was firstly divided into two different temporal spans known as the reference and affected periods. The reference period is the first temporal span in which climate controls the hydrological responses. Then, the statistical behavior of the time-serries changes due to human activities, and the affected periods. Two hydrological models, Soil and Water Assessment Tool (SWAT) and a black box Artificial Neural Networks (ANNs), were used to simulate the streamflow in the watershed. However, the results of the hydrological models showed their general acceptable performance to simulate the recorded streamflow at Tamar hydrometric station, but the results of the conceptual model (SWAT) showed that the performance of the model in the dry season is not as good as in the wet season. In the next step, the contributions of human and climate activities were assessed via two different methods. The first method is simple differential method, which compares the projection of the calibrated model in the second period with observations in both periods. The second set of contribution rates was calculated using the climate elasticity method via recorded monthly data and implemented derivation rules. In the first method, the contribution rate of human activities is significantly higher than the rate of climate change, and the result of the second method is a reverse. Because of differences in the methods’ concepts, the calculated contributions rates are different. To assess the uncertainty grouped with the estimations, a novel approach was developed using fuzzy mathematics. The uncertain version of the contribution rates showed that in each <em>α</em>-cut (fuzzy uncertainty level), the contribution of human alternation (LULC change) as the most important human interventions is more significant than climate direvers. In other words, during the simulation period, the effect of LULC change on the flow was very noteworthy, while climate change had relatively less effect on the behavioral change of the flow.Human and climate are two major scio-hydrologic drivers that determine hydrological regimes and patterns. In this regards, Land Use and Land Cover (LULC) changes, agricultural development, etc., on global and regional scales, hydro-climatological components have influenced these regimes. The effects of each driver on the variation of hydrological components have been assessed in different studies, but these approaches are not accurate enough at watershed-scales that experience the simultaneous impacts of climate dynamics and LULC changes. Different studies have considered both climate and human altertions in the hydrological cycle, and quantified their contributions in such basin. Results of these researches can help decision makers in water management of the pros and cons of water and land use policies. The Gorganrood watershed is an important basin in the northern part of Iran, especially from the agricultural point of view, which has considerably experienced hydrological and extreme events changes. While the consequence of each climate change and LULC changes have been assessed in the watershed, there is no study, which considers the complicated interactions of these drivers. In this paper, the authors firstly evaluated the contributions of LULC and climate change on the variation of streamflow. Secondly, the modified fuzzy arithmetic method has been used to achieve their fuzzy contributions. To this purpose, the computational period was firstly divided into two different temporal spans known as the reference and affected periods. The reference period is the first temporal span in which climate controls the hydrological responses. Then, the statistical behavior of the time-serries changes due to human activities, and the affected periods. Two hydrological models, Soil and Water Assessment Tool (SWAT) and a black box Artificial Neural Networks (ANNs), were used to simulate the streamflow in the watershed. However, the results of the hydrological models showed their general acceptable performance to simulate the recorded streamflow at Tamar hydrometric station, but the results of the conceptual model (SWAT) showed that the performance of the model in the dry season is not as good as in the wet season. In the next step, the contributions of human and climate activities were assessed via two different methods. The first method is simple differential method, which compares the projection of the calibrated model in the second period with observations in both periods. The second set of contribution rates was calculated using the climate elasticity method via recorded monthly data and implemented derivation rules. In the first method, the contribution rate of human activities is significantly higher than the rate of climate change, and the result of the second method is a reverse. Because of differences in the methods’ concepts, the calculated contributions rates are different. To assess the uncertainty grouped with the estimations, a novel approach was developed using fuzzy mathematics. The uncertain version of the contribution rates showed that in each <em>α</em>-cut (fuzzy uncertainty level), the contribution of human alternation (LULC change) as the most important human interventions is more significant than climate direvers. In other words, during the simulation period, the effect of LULC change on the flow was very noteworthy, while climate change had relatively less effect on the behavioral change of the flow.https://jesphys.ut.ac.ir/article_87081_9470f00a94d118d2de6807b7d3c54305.pdfInstitute of Geophysics, University of TehranJournal of the Earth and Space Physics2538-371X48220220823Investigation of Seasonal dust in northeastern Iran and numerical simulation of extreme dust events using WRF-CHEM modelInvestigation of Seasonal dust in northeastern Iran and numerical simulation of extreme dust events using WRF-CHEM model4214408691010.22059/jesphys.2022.330319.1007361FAAzarZarrinAssociate Professor, Department of Geography, Ferdowsi University of Mashhad, Mashhad, Iran0000-0002-4542-3176NargesSalehabadiM.Sc. Graduated, Department of Geography, Ferdowsi University of Mashhad, Mashhad, Iran0000-0002-8322-5277AbbasMofidiAssistant Professor, Department of Geography, Ferdowsi University of Mashhad, Mashhad, Iran0000-0002-5050-0033Abbas AliDadashi-RoudbariPost-Doc Researcher, Department of Geography, Ferdowsi University of Mashhad, Mashhad, Iran0000-0002-9308-1019Journal Article20210915In recent years, dust storm has become a serious environmental concern and has attracted a lot of attention among atmospheric scientists. Northeast of Iran is a large and strategic population area. Due to its proximity to large arid regions in Central Asia, this region has a high risk of experiensing dust events. In recent years, it has faced many problems regarding dust phenomena. This study is conducted to investigate seasonal dust events in northeastern Iran. To achieve this goal, a combination of station data, reanalysis, satellite and output of the WRF-Chem numerical model have been used simultaneously to improve our understanding of the dust seasonal cycle in northeast of Iran. Accordingly, this research was organized in two parts: monitoring and modeling of dust phenomenon. The results of this study may be useful for forecasting dust storms as well as spatial planning.<br />To investigate the dust events seasonal variabilities, the dust surface mass concentration of MERRA-2 dataset, aerosol optical depth (AOD) of the combined Dark Target (DT) and Deep Blue (DB) algorithms of MODIS sensor of Terra and Aqua satellites were examined during the long-term period (2004-2018).<br />Since the emission of dust is highly dependent on biophysical components, it is necessary to use numerical models. The WRF-Chem numerical model was used for this purpose. The study area includes northeastern Iran and parts of Central Asia. The horizontal resolution of the child domin of 30 km model was performed with 32 vertical levels. The NCEP / FNL is used as boundary conditions with 3-hourly time step and 1-degree horizontal resolution for the model configuration. Four extreme dust events were selected to investigate the transport of dust to northeastern Iran. The selected dust events occurred on November 13, 2007, May 29, 2008, June 8, 2015, and October 17, 2017 in northeastern Iran. Therefore, case events were simulated with a time step of 180 seconds and output every three hours using GOCART, AFWA, UoC_S01 and UoC_S11 schemes.<br />The results showed that the maximum dust activity occurred in spring with AOD value equal to 0.59 and dust surface mass concentration is 645.2 µg m <sup>-3</sup>. The summer is in the next place. Seasonal analysis of AOD and dust using satellite and reanalysis data, showed that Aralkum, Kyzylkum, Karakum and Kara-Bogaz-Gol are the main dust sources in Central Asia that are active in all seasons.<br />Comparison of dust simulation results for PM<sub>2.5 </sub>and PM<sub>10</sub> variable with observational data of air quality control stations in Mashhad showed that GOCART scheme can well depict dust events and provide a low bias towards station data. Also, the study of correlation coefficient between simulation and observation showed that the GOCART scheme explains nearly 90% of the variance of data. The root mean square error (RMSE) for the GOCART scheme is less than 20 micrograms per cubic meter for PM<sub>2.5</sub>. Accordingly, the GOCART scheme is a suitable scheme for dust study in Northeast of Iran and the WRF-Chem model can be used to operationally forecast dust storms. The dust detection algorithm (DDA) of the AIRS sensor and the aerosol optical depth (AOD) of the MODIS sensor confirm the contribution of the mentioned sources of dust in transferring dust to the northeast of Iran. The results showed that three of the case studies occurred as a result of the passage of an extratropical Rossby wave and the deepening of the trough on the territory of Turkmenistan. In contrast, the summer case study is the result of the establishment of a summer circulation pattern that has occurred with the simultaneous establishment of an anticyclonic circulation in the southern part of Turkmenistan and the northeastern parts of Iran and a cyclonic circulation in the Sistan plain and southeastern parts of the country.In recent years, dust storm has become a serious environmental concern and has attracted a lot of attention among atmospheric scientists. Northeast of Iran is a large and strategic population area. Due to its proximity to large arid regions in Central Asia, this region has a high risk of experiensing dust events. In recent years, it has faced many problems regarding dust phenomena. This study is conducted to investigate seasonal dust events in northeastern Iran. To achieve this goal, a combination of station data, reanalysis, satellite and output of the WRF-Chem numerical model have been used simultaneously to improve our understanding of the dust seasonal cycle in northeast of Iran. Accordingly, this research was organized in two parts: monitoring and modeling of dust phenomenon. The results of this study may be useful for forecasting dust storms as well as spatial planning.<br />To investigate the dust events seasonal variabilities, the dust surface mass concentration of MERRA-2 dataset, aerosol optical depth (AOD) of the combined Dark Target (DT) and Deep Blue (DB) algorithms of MODIS sensor of Terra and Aqua satellites were examined during the long-term period (2004-2018).<br />Since the emission of dust is highly dependent on biophysical components, it is necessary to use numerical models. The WRF-Chem numerical model was used for this purpose. The study area includes northeastern Iran and parts of Central Asia. The horizontal resolution of the child domin of 30 km model was performed with 32 vertical levels. The NCEP / FNL is used as boundary conditions with 3-hourly time step and 1-degree horizontal resolution for the model configuration. Four extreme dust events were selected to investigate the transport of dust to northeastern Iran. The selected dust events occurred on November 13, 2007, May 29, 2008, June 8, 2015, and October 17, 2017 in northeastern Iran. Therefore, case events were simulated with a time step of 180 seconds and output every three hours using GOCART, AFWA, UoC_S01 and UoC_S11 schemes.<br />The results showed that the maximum dust activity occurred in spring with AOD value equal to 0.59 and dust surface mass concentration is 645.2 µg m <sup>-3</sup>. The summer is in the next place. Seasonal analysis of AOD and dust using satellite and reanalysis data, showed that Aralkum, Kyzylkum, Karakum and Kara-Bogaz-Gol are the main dust sources in Central Asia that are active in all seasons.<br />Comparison of dust simulation results for PM<sub>2.5 </sub>and PM<sub>10</sub> variable with observational data of air quality control stations in Mashhad showed that GOCART scheme can well depict dust events and provide a low bias towards station data. Also, the study of correlation coefficient between simulation and observation showed that the GOCART scheme explains nearly 90% of the variance of data. The root mean square error (RMSE) for the GOCART scheme is less than 20 micrograms per cubic meter for PM<sub>2.5</sub>. Accordingly, the GOCART scheme is a suitable scheme for dust study in Northeast of Iran and the WRF-Chem model can be used to operationally forecast dust storms. The dust detection algorithm (DDA) of the AIRS sensor and the aerosol optical depth (AOD) of the MODIS sensor confirm the contribution of the mentioned sources of dust in transferring dust to the northeast of Iran. The results showed that three of the case studies occurred as a result of the passage of an extratropical Rossby wave and the deepening of the trough on the territory of Turkmenistan. In contrast, the summer case study is the result of the establishment of a summer circulation pattern that has occurred with the simultaneous establishment of an anticyclonic circulation in the southern part of Turkmenistan and the northeastern parts of Iran and a cyclonic circulation in the Sistan plain and southeastern parts of the country.https://jesphys.ut.ac.ir/article_86910_7a4f678baea270f84789e9dc5e18a8c4.pdfInstitute of Geophysics, University of TehranJournal of the Earth and Space Physics2538-371X48220220823Statistical modeling of the mean annual temperature at Mehrabad station, TehranStatistical modeling of the mean annual temperature at Mehrabad station, Tehran4414528544710.22059/jesphys.2022.332720.1007372FAArmanJahediPh.D. Student, Department of Geography, Faculty of Humanities, University of Zanjan, Zanjan, Iran0000-0002-9534-9058Journal Article20211027Regarding climate changes and global warming, it seems that the behavior of climate elements in the future should be predicted and known. Therefore, in this study, using modeling by a set of ARIMA statistical models, models on the time series of the mean annual temperature at Mehrabad station in Tehran during 1951-2015 were fitted to investigate a significant model by trial and error in order to identify the most appropriate model. Since the time series of the observations had a normal distribution, modeling was performed on the time series without applying Box Cox transformation. First, for static and non-static investigations, the time series of annual mean temperature observations was plotted simply in diagrams. In addition, the first and second order regression line equations were used to further ensure the type of time series behavior of the mean annual temperature. The results showed that the time series behavior of temperature at this station is linear. Since the time series behavior was linear, the order d = 1 was determined. Second, the first-order differentiation was performed on the time series. In the third step, the order p and q were determined using autocorrelation and partial autocorrelation of the differentiated values ( ). After investigating the significance of the order of the components of each of the models, the following models were selected as significant models, respectively:<br /><em>1) ARIMA(0,1,1</em><br /><em>2) ARIMA</em><em>(2,1,0</em><br />Since the first significant model was observed with suspicion, as a result each of the components (p, d, q) of the above two models were tested up to the 3<sup>th</sup> order. Finally, these two models were selected as significant models. Also, Akaike information criterion (AIC) was considered to determine the most appropriate model among the above two models. ARIMA model (0,1,1 had the minimum value of AIC compared to the other model. As a result, using this model, the temperature time series at this station was predicted from the end of the period to ¼ of the first time series. Given the concept of uncertainty, which underlies descriptive and inferential statistics, as a result, it seems that uncertainties should be expressed with high statistical certainty. In this regard, we used statistical tests of autocorrelation, Pearson correlation coefficient, standard normal homogeneity, cumulative deviations, milestones, sign on the time series of ARIMA model residues (0,1,1 , and drawing methods for residual normality, residual independence, constant residual variance and <em>portmanteau </em>test to consider further criteria to increase the statistical reliability of the applied model. The results of all statistical tests showed the random residual time series of the model. These tests showed that the best model for modeling the time series of the mean annual temperature at Mehrabad station, Tehran is ARIMA model (0,1,1 . Since the upper and lower limits of the predicted series as well as the predicted observations show the same behavior of the temperature time series at Mehrabad station, it can be said that the estimation of the predicted numerical values is still appropriate for this model to predict the temperature variable at this station. Finally, the results showed that the mean temperature of the predicted series is likely to be 17.742 ͦ C, and the mean annual temperature will increase by 0.038 ͦ C compared to the previous year.Regarding climate changes and global warming, it seems that the behavior of climate elements in the future should be predicted and known. Therefore, in this study, using modeling by a set of ARIMA statistical models, models on the time series of the mean annual temperature at Mehrabad station in Tehran during 1951-2015 were fitted to investigate a significant model by trial and error in order to identify the most appropriate model. Since the time series of the observations had a normal distribution, modeling was performed on the time series without applying Box Cox transformation. First, for static and non-static investigations, the time series of annual mean temperature observations was plotted simply in diagrams. In addition, the first and second order regression line equations were used to further ensure the type of time series behavior of the mean annual temperature. The results showed that the time series behavior of temperature at this station is linear. Since the time series behavior was linear, the order d = 1 was determined. Second, the first-order differentiation was performed on the time series. In the third step, the order p and q were determined using autocorrelation and partial autocorrelation of the differentiated values ( ). After investigating the significance of the order of the components of each of the models, the following models were selected as significant models, respectively:<br /><em>1) ARIMA(0,1,1</em><br /><em>2) ARIMA</em><em>(2,1,0</em><br />Since the first significant model was observed with suspicion, as a result each of the components (p, d, q) of the above two models were tested up to the 3<sup>th</sup> order. Finally, these two models were selected as significant models. Also, Akaike information criterion (AIC) was considered to determine the most appropriate model among the above two models. ARIMA model (0,1,1 had the minimum value of AIC compared to the other model. As a result, using this model, the temperature time series at this station was predicted from the end of the period to ¼ of the first time series. Given the concept of uncertainty, which underlies descriptive and inferential statistics, as a result, it seems that uncertainties should be expressed with high statistical certainty. In this regard, we used statistical tests of autocorrelation, Pearson correlation coefficient, standard normal homogeneity, cumulative deviations, milestones, sign on the time series of ARIMA model residues (0,1,1 , and drawing methods for residual normality, residual independence, constant residual variance and <em>portmanteau </em>test to consider further criteria to increase the statistical reliability of the applied model. The results of all statistical tests showed the random residual time series of the model. These tests showed that the best model for modeling the time series of the mean annual temperature at Mehrabad station, Tehran is ARIMA model (0,1,1 . Since the upper and lower limits of the predicted series as well as the predicted observations show the same behavior of the temperature time series at Mehrabad station, it can be said that the estimation of the predicted numerical values is still appropriate for this model to predict the temperature variable at this station. Finally, the results showed that the mean temperature of the predicted series is likely to be 17.742 ͦ C, and the mean annual temperature will increase by 0.038 ͦ C compared to the previous year.https://jesphys.ut.ac.ir/article_85447_4d16bc049ab29bade76c1c5986b0f171.pdfInstitute of Geophysics, University of TehranJournal of the Earth and Space Physics2538-371X48220220823Effect of non-thermal and trapped electrons on solitary waves and chaos in auroral acceleration regionsEffect of non-thermal and trapped electrons on solitary waves and chaos in auroral acceleration regions4534758690210.22059/jesphys.2022.329937.1007373FAMojtabaHashemzadeh DehaghaniAssistant Professor, Department of Plasma Physics and Basic Particle, Faculty of Physics, Shahrood University of Technology, Shahrood, Iran0000-0001-5580-563XJournal Article20211030In this paper, using the reductive perturbation method, the propagation of nonlinear solitary waves and chaos phenomenon and its stability were studied in auroral acceleration regions in the presence of electrons with the Cairns-Gurevich distribution function. Using the continuity, momentum transfer, and Poisson equations, and considering the density of electrons as the Cairns-Gurovich distribution function, and using two different models, Korteweg–De Vries (KdV) and modified KdV equations were obtained. It was shown that the solutions of these equations are in the form of solitary waves. The effect of non-thermal and trapped electrons and wave velocity on these waves were studied. In the next section, pseudo-potentials and total mechanical energy are obtained. Considering a quasi-periodic factor, KdV and modified KdV equations were reviewed and the chaos and its stability were studied in the auroral acceleration regions. Results showed that by increasing the wave velocity and non-thermal and trapped parameters, the size of the field increased, and the depth of the potential well was also increased. These results confirmed each other. It was indicated that in the case of b=0, this distribution function became as the Maxwellian distribution function. In the case b>0, in addition to free particles, the trapped and non-thermal particles also affect the distribution function. In this case, the width of the distribution function became larger, which indicated that the more energetic electrons existed in this case. It is also concluded that for both nonlinear equations, the solutions can exist in the form of rarefactive and compressive solitons. Three-dimensional graphs of total mechanical energy were also plotted for different values of the wave velocity and non-thermal and trapped parameters. Results for this case also showed that for the total energy of E<sub>1</sub>, by increasing the b parameter, the energy deviated from the uniform function and reached the saddle state. It was also shown that the wave velocity was similar to the b parameter. It was found that for different values of U and b parameters, the behavior of the total energy of E<sub>2</sub> was different from the total energy diagram of E<sub>1</sub>. Poincaré return mapping diagrams confirmed the existence of a closed cycle indicating chaos in these plasmas. Results of this section also showed that for solitons with function ψ<sub>1</sub>, by increasing the U parameter, the Poincaré return mapping cycle region increased. Poincaré return mapping lines were also more focused in this case. For solitons with ψ<sub>1</sub> functions, by increasing the wave velocity, Poincaré's return map goes from a quasi-stable state to a stable state. By increasing the quasi-periodic frequency, the Poincaré return map goes from steady-state to quasi-steady state so that a cycle converts to two cycles with a certain overlap. Finally, it was concluded that using real parameters, the wave velocity was in the interval 13km/s<v'<52km/s and the electric field was approximately 5mV/m and the Debye length became 15m. It was also concluded that the results of the recent work were in good agreement with the results obtained from the Viking, Freja and S3-3 satellites.In this paper, using the reductive perturbation method, the propagation of nonlinear solitary waves and chaos phenomenon and its stability were studied in auroral acceleration regions in the presence of electrons with the Cairns-Gurevich distribution function. Using the continuity, momentum transfer, and Poisson equations, and considering the density of electrons as the Cairns-Gurovich distribution function, and using two different models, Korteweg–De Vries (KdV) and modified KdV equations were obtained. It was shown that the solutions of these equations are in the form of solitary waves. The effect of non-thermal and trapped electrons and wave velocity on these waves were studied. In the next section, pseudo-potentials and total mechanical energy are obtained. Considering a quasi-periodic factor, KdV and modified KdV equations were reviewed and the chaos and its stability were studied in the auroral acceleration regions. Results showed that by increasing the wave velocity and non-thermal and trapped parameters, the size of the field increased, and the depth of the potential well was also increased. These results confirmed each other. It was indicated that in the case of b=0, this distribution function became as the Maxwellian distribution function. In the case b>0, in addition to free particles, the trapped and non-thermal particles also affect the distribution function. In this case, the width of the distribution function became larger, which indicated that the more energetic electrons existed in this case. It is also concluded that for both nonlinear equations, the solutions can exist in the form of rarefactive and compressive solitons. Three-dimensional graphs of total mechanical energy were also plotted for different values of the wave velocity and non-thermal and trapped parameters. Results for this case also showed that for the total energy of E<sub>1</sub>, by increasing the b parameter, the energy deviated from the uniform function and reached the saddle state. It was also shown that the wave velocity was similar to the b parameter. It was found that for different values of U and b parameters, the behavior of the total energy of E<sub>2</sub> was different from the total energy diagram of E<sub>1</sub>. Poincaré return mapping diagrams confirmed the existence of a closed cycle indicating chaos in these plasmas. Results of this section also showed that for solitons with function ψ<sub>1</sub>, by increasing the U parameter, the Poincaré return mapping cycle region increased. Poincaré return mapping lines were also more focused in this case. For solitons with ψ<sub>1</sub> functions, by increasing the wave velocity, Poincaré's return map goes from a quasi-stable state to a stable state. By increasing the quasi-periodic frequency, the Poincaré return map goes from steady-state to quasi-steady state so that a cycle converts to two cycles with a certain overlap. Finally, it was concluded that using real parameters, the wave velocity was in the interval 13km/s<v'<52km/s and the electric field was approximately 5mV/m and the Debye length became 15m. It was also concluded that the results of the recent work were in good agreement with the results obtained from the Viking, Freja and S3-3 satellites.https://jesphys.ut.ac.ir/article_86902_086f612930fc5761f43a66c2b2540418.pdfInstitute of Geophysics, University of TehranJournal of the Earth and Space Physics2538-371X48220220823Numerical solution of two-layer shallow water equations using mode splitting methodNumerical solution of two-layer shallow water equations using mode splitting method4774948690710.22059/jesphys.2022.334561.1007384FAHakimGolshahyAssistant Professor, Department of Physics, Shoushtar Branch, Islamic Azad University, Shoushtar, Iran0000-0002-3555-2791Journal Article20211205In the numerical models that use iterative methods to solve the momentum equations by applying the rigid-lid approximation, the number of iterations increases for high resolution, therefore processing time increases. An alternative method is applying a free surface and splitting equations to barotropic and baroclinic modes. The surface gravity waves that are faster than slow moving internal gravity waves; impose a limitation on the time steps with the CFL condition. Thus, mode splitting method is computationally efficient to handle the multiple time steps separating the barotropic and baroclinic mode equations. In this method, the barotropic mode equations are solved at small time steps consistent with the fast surface gravity wave speeds and the baroclinic mode equations are solved at larger time steps consistent with the slow internal gravity wave speeds. This is used in most of the ocean circulation models and is an unavoidable choice to high resolution models.<br />In this study, we considered the shallow water equations for two-layer basin with vorticity-divergence formulation using mode splitting method by a small time step of barotropic mode within a larger time step of baroclinic mode. The primary systems of equations that contain both upper and lower layer variables, are rewritten in terms of new (barotropic and baroclinic) variables without any variation or more approximation of primary systems. This procedure can be extended to multi-layer systems so that primary N-layer system of equations is changed to 1 system of barotropic mode equations and N-1 systems of baroclinic mode equations coupled together.<br />For numerical experiments, a fully baroclinic (non-barotropic) initial condition is considered in a constant depth rectangular domain with 64, 128 and 256 grid points in each direction and periodic boundaries. For spatial differencing, second order centered scheme with low computational cost and fourth-order compact scheme with high computational cost are used. For time integration, a semi-implicit descretization based on leapfrog scheme is implemented with the Robert-Asselin time filter for both barotropic and baroclinic systems of equations, similarly.<br />Mode splitting method may presents numerical instabilities on the larger baroclinic time steps, in spite of time step limitation based on CFL condition coming from each system of barotropic and baroclinic mode equations taken individually. Here, it is controlled by increasing the coefficient of time filter to some extent.<br />First, we solve the baroclinic mode equations to derive all baroclinic variables that are necessary to solve barotropic mode equations during a baroclinic time step. In this case, these variables can be taken to be constant up to the next baroclinic time level or determined by time interpolation between two successive baroclinic time levels.<br />To assess the performance of the numerical method, relative error of energy conservation is calculated. Results show that for the ratio of baroclinic time step to up to 20 times of that of barotropic time step, time evolution of the barotropic and baroclinic variables have appropriate correspondence to the basic state, in which the barotropic mode has the same time step as the baroclinic mode. When this ratio increases, the differences of errors from basic state are presented more clearly. These errors are increased on fourth order compact method insofar as it leads to numerical instabilities so the time filter coefficient had to be increased, while second order scheme is not sensitive and stays stable with small coefficient. Moreover, taking constant for baroclinic variables to solve barotropic mode equations makes the solution on fourth order compact scheme for large baroclinic time step unstable, but on the other hand time interpolation provides more stable condition and has a good performance almost on both spatial schemes.In the numerical models that use iterative methods to solve the momentum equations by applying the rigid-lid approximation, the number of iterations increases for high resolution, therefore processing time increases. An alternative method is applying a free surface and splitting equations to barotropic and baroclinic modes. The surface gravity waves that are faster than slow moving internal gravity waves; impose a limitation on the time steps with the CFL condition. Thus, mode splitting method is computationally efficient to handle the multiple time steps separating the barotropic and baroclinic mode equations. In this method, the barotropic mode equations are solved at small time steps consistent with the fast surface gravity wave speeds and the baroclinic mode equations are solved at larger time steps consistent with the slow internal gravity wave speeds. This is used in most of the ocean circulation models and is an unavoidable choice to high resolution models.<br />In this study, we considered the shallow water equations for two-layer basin with vorticity-divergence formulation using mode splitting method by a small time step of barotropic mode within a larger time step of baroclinic mode. The primary systems of equations that contain both upper and lower layer variables, are rewritten in terms of new (barotropic and baroclinic) variables without any variation or more approximation of primary systems. This procedure can be extended to multi-layer systems so that primary N-layer system of equations is changed to 1 system of barotropic mode equations and N-1 systems of baroclinic mode equations coupled together.<br />For numerical experiments, a fully baroclinic (non-barotropic) initial condition is considered in a constant depth rectangular domain with 64, 128 and 256 grid points in each direction and periodic boundaries. For spatial differencing, second order centered scheme with low computational cost and fourth-order compact scheme with high computational cost are used. For time integration, a semi-implicit descretization based on leapfrog scheme is implemented with the Robert-Asselin time filter for both barotropic and baroclinic systems of equations, similarly.<br />Mode splitting method may presents numerical instabilities on the larger baroclinic time steps, in spite of time step limitation based on CFL condition coming from each system of barotropic and baroclinic mode equations taken individually. Here, it is controlled by increasing the coefficient of time filter to some extent.<br />First, we solve the baroclinic mode equations to derive all baroclinic variables that are necessary to solve barotropic mode equations during a baroclinic time step. In this case, these variables can be taken to be constant up to the next baroclinic time level or determined by time interpolation between two successive baroclinic time levels.<br />To assess the performance of the numerical method, relative error of energy conservation is calculated. Results show that for the ratio of baroclinic time step to up to 20 times of that of barotropic time step, time evolution of the barotropic and baroclinic variables have appropriate correspondence to the basic state, in which the barotropic mode has the same time step as the baroclinic mode. When this ratio increases, the differences of errors from basic state are presented more clearly. These errors are increased on fourth order compact method insofar as it leads to numerical instabilities so the time filter coefficient had to be increased, while second order scheme is not sensitive and stays stable with small coefficient. Moreover, taking constant for baroclinic variables to solve barotropic mode equations makes the solution on fourth order compact scheme for large baroclinic time step unstable, but on the other hand time interpolation provides more stable condition and has a good performance almost on both spatial schemes.https://jesphys.ut.ac.ir/article_86907_34b156895712684a83bcfe5ed08a694d.pdfInstitute of Geophysics, University of TehranJournal of the Earth and Space Physics2538-371X48220220823Assessment of the performance of cumulus and boundary layer schemes in the WRF-NMM model in simulation of heavy rainfalls over the Bushehr Province during 2000-2020Assessment of the performance of cumulus and boundary layer schemes in the WRF-NMM model in simulation of heavy rainfalls over the Bushehr Province during 2000-20204955148690110.22059/jesphys.2022.334737.1007387FANafisehPegahfarAssistant Professor, Atmospheric Science Center, Iranian National Institute for Oceanography and Atmospheric Science, Tehran, Iran0000-0003-4885-7428Journal Article20211205The mesoscale numerical weather prediction system of Weather Research and Forecasting (WRF), with two cores of ARW and NMM, has been used for atmospheric research, operational forecasting, and dynamical downscaling of Global Climate Models. Many parameterizations for each physics option can be accessed in this model. It is noteworthy that the performance of the model depends on the selected configuration and varies in different areas. Therefore, choosing a configuration with the lowest error for each terrain is mandatory. Here, the performances of various physics schemes, including cumulus and boundary layer schemes of the WRF-NMM model, were examined to simulate twelve heaviest extreme rainfall events in the southwest of Iran, the Bushehr Province, during 2000-2020. These events lasted for eighteen days. Three domains with 27, 9, and 3 km resolution were used in the model configurations, with no cumulus option for the smallest one. The initial and boundary conditions were used from the European Centre for Medium-Range Weather Forecasts (ECMWF) Reanalysis v5 (ERA5) datasets. One hundred and eight simulations were done using six cumulus schemes of KF, BMJ, SAS, oldSAS, NSAS, and TiedTKE, and seventy-two runs were done to evaluate the boundary layer schemes of MRF, MYJ, QNSE, and YSU. The simulated precipitation patterns were assessed using two observational data sets, including (I) in-situ measured data from eleven automatic weather stations and (II) grid point data from Global Precipitation Measurement (GPM) satellite with 0.1-degree horizontal resolution. Four statistic indices of Root Mean Square Error, Correlation Coefficient, Standard Deviation, and Bias were applied in the evaluation process. The evaluation process with the data measured at 11 automatic weather stations was done using outputs of the third domain. The outputs of the second domain were used for evaluation basis on GPM data at grid points. For a comprehensive analysis, the assessment process was performed separately for rainfall events (March-April and November-December events) in coastal and non-coastal stations. Comparison of precipitation from simulations of various cumulus schemes with the eleven in-situ data showed that the schemes from SAS family well performed at March-April events at coastal and noncoastal stations. While, the KF scheme produced the least error at coastal and noncoastal stations during the November-December events. The precipitation data from 1271 GPM grid-point data revealed that the oldSAS scheme generated the least error for the March-April and November-December events. According to the number of GPM grid-point data, the oldSAS scheme opted as the cumulus option for the next runs. Evaluation of WRF-NMM simulations using different boundary layer physics with the in-situ data indicated that MRF scheme produced the minor error at coastal and noncoastal stations for both March-April and November-December events. Using the 1271 GPM grid-point data illustrated that the QNSE and MRF (MYJ and MRF) options did the best performance for March-April (November-December) events. In conclusion, based on the number of GPM grid-point data compared with in-situ measured data, it is suggested that the oldSAS cumulus scheme and MRF boundary layer scheme can be chosen with some robustness in predicting the amount and pattern of the heavy rainfall precipitation in Bushehr Province of Iran. It is also notable that the default options introduced by the model for cumulus scheme and boundary layer scheme in the WRF-NMM model produce the largest error and are not appropriate for the selected area. This reveals the importance of adequately selecting physics options for this area.The mesoscale numerical weather prediction system of Weather Research and Forecasting (WRF), with two cores of ARW and NMM, has been used for atmospheric research, operational forecasting, and dynamical downscaling of Global Climate Models. Many parameterizations for each physics option can be accessed in this model. It is noteworthy that the performance of the model depends on the selected configuration and varies in different areas. Therefore, choosing a configuration with the lowest error for each terrain is mandatory. Here, the performances of various physics schemes, including cumulus and boundary layer schemes of the WRF-NMM model, were examined to simulate twelve heaviest extreme rainfall events in the southwest of Iran, the Bushehr Province, during 2000-2020. These events lasted for eighteen days. Three domains with 27, 9, and 3 km resolution were used in the model configurations, with no cumulus option for the smallest one. The initial and boundary conditions were used from the European Centre for Medium-Range Weather Forecasts (ECMWF) Reanalysis v5 (ERA5) datasets. One hundred and eight simulations were done using six cumulus schemes of KF, BMJ, SAS, oldSAS, NSAS, and TiedTKE, and seventy-two runs were done to evaluate the boundary layer schemes of MRF, MYJ, QNSE, and YSU. The simulated precipitation patterns were assessed using two observational data sets, including (I) in-situ measured data from eleven automatic weather stations and (II) grid point data from Global Precipitation Measurement (GPM) satellite with 0.1-degree horizontal resolution. Four statistic indices of Root Mean Square Error, Correlation Coefficient, Standard Deviation, and Bias were applied in the evaluation process. The evaluation process with the data measured at 11 automatic weather stations was done using outputs of the third domain. The outputs of the second domain were used for evaluation basis on GPM data at grid points. For a comprehensive analysis, the assessment process was performed separately for rainfall events (March-April and November-December events) in coastal and non-coastal stations. Comparison of precipitation from simulations of various cumulus schemes with the eleven in-situ data showed that the schemes from SAS family well performed at March-April events at coastal and noncoastal stations. While, the KF scheme produced the least error at coastal and noncoastal stations during the November-December events. The precipitation data from 1271 GPM grid-point data revealed that the oldSAS scheme generated the least error for the March-April and November-December events. According to the number of GPM grid-point data, the oldSAS scheme opted as the cumulus option for the next runs. Evaluation of WRF-NMM simulations using different boundary layer physics with the in-situ data indicated that MRF scheme produced the minor error at coastal and noncoastal stations for both March-April and November-December events. Using the 1271 GPM grid-point data illustrated that the QNSE and MRF (MYJ and MRF) options did the best performance for March-April (November-December) events. In conclusion, based on the number of GPM grid-point data compared with in-situ measured data, it is suggested that the oldSAS cumulus scheme and MRF boundary layer scheme can be chosen with some robustness in predicting the amount and pattern of the heavy rainfall precipitation in Bushehr Province of Iran. It is also notable that the default options introduced by the model for cumulus scheme and boundary layer scheme in the WRF-NMM model produce the largest error and are not appropriate for the selected area. This reveals the importance of adequately selecting physics options for this area.https://jesphys.ut.ac.ir/article_86901_4f3fff5af299daba9dbfe426938b93dd.pdf