Estimating the slip rate in the North Tabriz Fault using focal mechanism data and GPS velocity field

Milad Salmanian; Asghar Rastbood; Masoud Mashhadi Hossainali

doi:10.1515/jogs-2022-0167

Article Open Access

Estimating the slip rate in the North Tabriz Fault using focal mechanism data and GPS velocity field

Milad Salmanian , Asghar Rastbood and Masoud Mashhadi Hossainali

Published/Copyright: March 19, 2024

Published by

Become an author with De Gruyter Brill

Author Information

From the journal Journal of Geodetic Science Volume 14 Issue 1

Abstract

This study delves into slip distribution on the North Tabriz Fault (NTF), a critical aspect of seismic hazard analysis due to its proximity to the Tabriz metropolis. The study operates within a uniform elastic half-space, maintaining constant values for fault geometry and regional rheological parameters throughout the research. To calculate strain boundary conditions, permanent and periodic global positioning system (GPS) data from the northwest region were utilized. The fault was constrained perpendicularly while allowing tangential movement, facilitating the determination of its annual slip rate using the boundary element method, with the Okada analytical model serving as the fundamental solution. The findings underscore the intricate relationship between the fault’s slip rate and boundary conditions, revealing a predominant right-lateral strike-slip motion. The study offers two slip rate assessments, obtained through earthquake focal mechanisms and GPS velocity data, yielding values of 5 and 5.5 mm / year , respectively. Importantly, the alignment of these calculated slip rates with paleo-seismological data underscores the credibility of the results generated via the boundary element method, distinguishing it as a reliable approach when compared to other numerical and analytical techniques. This research provides valuable insights into the behavior and slip dynamics of the NTF, which is pivotal for assessing seismic risks.

Keywords: boundary element method; focal mechanism; GPS velocity field; North Tabriz fault; slip rate

1 Introduction

The slip rate of the fault is a crucial factor in assessing the seismic risk within the study area. As the North Tabriz fault (NTF) holds paramount importance among the faults in Iran’s northwestern region and is located near the populous city of Tabriz, accurately determining its slip rate becomes highly significant in seismic risk estimation. This study aimed to employ an optimal and efficient inversion method for slip rate determination. By utilizing the most up-to-date data, the study performs the inversion operation using two distinct boundary conditions: the strain field obtained from global positioning system (GPS) velocity field data and the stress field obtained from earthquake focal mechanism data. The convergence of estimates obtained from both boundary conditions and their consistency with estimates derived from paleoseismological methods validate the results of this study. The use of contemporary techniques to determine the deformation of the Earth’s crust is of great importance in the field of tectonic research. By integrating geological, geodetic, and seismic data, fault models can be tailored to comprehend the characteristics and behavior of faults over time in crustal deformation studies. Such models have the potential to forecast future deformations and seismic hazards, which are especially critical for cities located in close proximity to active seismic areas, such as Tabriz, as per previous observations.

The youthful crust in northwestern Iran has been influenced by the northward motion of the Arabian plate, resulting in the deformation of the region. Since the Late Miocene, this motion has generated north-south compressive forces and east-west tensile forces, which have regulated the distribution of volcanoes and led to the predominance of strike-slip faults. The collision between the Arabian plate and Eurasia at the boundary of the Eocene-Oligocene is the crucial factor in the regulation of active deformation in Iran (Agard et al. 2011, Allen et al. 2004, Dewey et al. 1986, Jackson and McKenzie 1984, McKenzie 1972, Mouthereau et al. 2012). The antecedent scientific investigations concerning the GPS have furnished a thorough scrutiny of the dispersion of deformation in Iran (Khorrami et al. 2019, Rahmani et al. 2023).

According to Pollard (1987), the slip rate distribution along a fault should be symmetrical and elliptical if no external factors are taken into account. However, studies have revealed that several factors such as heterogeneous elastic properties, inelastic deformation process, the mechanical interaction between adjacent faults, and varying stress conditions in remote areas from the fault can significantly affect the slip rate along a fault. Bürgmann et al. (1994) found that these external factors can cause an increase in the slip rate along a fault, resulting in a shift toward a symmetric distribution.

The low loading rate of faults located in intracontinental deformation zones can result in seismicity transfer between different sections of the fault, even at considerable distances. These occurrences, which may transpire long after a seismic quiet period, can ultimately lead to the emergence of significant earthquakes. Consequently, the examination of mechanical interaction and the determination of slip rate based on this interaction are critical for investigating these areas.

The Tabriz fault’s northern side can be easily detected and delineated on both satellite images and aerial photographs. Such images reveal alterations in the late Quaternary alluvial cones and watercourses, brought about by right-lateral displacement. By restoring the images to their pre-displacement state, it can be observed that the canals dug along the late quaternary geomorphology exhibit a right-lateral horizontal shift (Rizza et al. 2013).

Karakhanian et al. (2004) conducted an investigation of the displacement patterns resulting from seismic activity in the north of Tabriz. They observed maximum horizontal displacements of 100 m in large valleys near the village of Zabarlo and minimum displacements 3–5 m in secondary waterways near the village of Khajeh Marjan, with at least 2–3 m of vertical displacement. The vertical displacements in the western part of the fault north of Tabriz were reported to be two to seven times less than the corresponding horizontal displacements, with a maximum vertical displacement of 10 m . Based on their observations of 8 m of horizontal displacement of small waterways along the main branch of the fault (with a maximum age of late Psistocene), total slip rate of 2 mm / year was estimated for this part of the fault.

Rizza et al. (2013) employed two geodetic observation techniques to derive the slip rate values for the NTF. Specifically, the slip rate values were estimated to be 7.3 mm / year based on GPS stations observations and 6.0 mm / year through radar interferometric analysis. In a separate study, Djamour et al. (2011) utilized the block model to estimate the slip rate of the NTF, which stood at 7.3 mm / year . In addition, Khodaverdian et al. (2015) developed a kinematic finite element modeling approach to investigate land surface deformations in the Iranian plateau. The model incorporated key inputs such as the geometry of faults, the slip rate of faults based on geological and geomorphological studies, the GPS velocity field during large earthquake intervals, and main stress directions and velocity values at the boundary points of the modeling area. The study estimated the sliding rate for the NTF at 5.8 mm / year .

Karimzadeh et al. (2013) utilized radar interferometry observations to model the simple elastic failure of short baseline time series and proposed an average slip rate of 8.7 mm / year for the NTF. In contrast, Hessami et al. (2003b) employed paleo-seismological methods in the northwestern part of the fault located in the west of Tabriz city and reported a sliding rate of 3.1 to 6.4 mm / year for this section of the fault based on various techniques.

The slip rates of the NTF were evaluated by Haji-Aghajany et al. (2019), employing troposphere-corrected InSAR measurements and the 3D ray tracing technique. The resulting slip rates were estimated to be 7.6 and 5.6 mm / year . Figure 1 presents a summary of the previous estimations regarding the slip rate for the NTF.

Figure 1

Slip rate values for the NTF obtained from different studies.

The proximity of the Tabriz metropolis to the NTF highlights the significance of investigating the mechanical interaction occurring on this fault, particularly for assessing seismic risk in the city. The present study employed the boundary element method to determine the slip rate on the NTF. Two approaches were utilized, namely, the first method based on the strain rate deduced from GPS observations and the second method employing the stress tensor derived from the focal mechanism data of earthquakes. The results from each method were compared with each other.

2 Tectonic setting

Iran, situated within the oblique collision zone of the Arabian and Eurasian tectonic plates, is experiencing intracontinental deformation. The topography of Iran has been shaped by this collision, which occurs at a rate of 22 ± 2 mm / year , with the northwestern region of Iran contributing 8 ± 2 mm / year (Vernant et al. 2004). This convergence is distributed among various processes, including shortening in the Zagros Mountains, internal deformation caused by predominantly strike-slip faults in central Iran, and shortening in the Alborz Mountains (Figure 2). The northwestern region of Iran, which includes the North Anatolian Fault, the Eastern Anatolian Fault, and the Caucasus Mountains, is part of a complex system due to the interaction between the Arabian plate, the Anatolian Plateau, and the Eurasian Plate. This complex fault system transfers some of the northward movement of the Arabian Plate to the Anatolian Plateau, while the tilting of the collision zone in the Zagros Mountains results in the release of movement between the shortening in the Caucasus and the occurrence of right-lateral strike-slip in the NTF (Jackson 1992).

Figure 2

The present study focuses on the region encompassing the NTF or the northwestern region of Iran. The study area has been demarcated based on the set scale for this research.

The intricate NTF is situated in the northwestern region of Iran and spans nearly 150 km (Figure 2), with the inclusion of the North Mishu fault as its northwestern sequence pushing its length to around 210 km. It commences from the southern area of Marand with a northwest-southeast trend and passes through the northern border of Tabriz city, extending near Bostanabad. The fault is linked to an overturned fault zone from the northwest side, which rotates to the west-southwest in the northern territory of Urmia Lake, known as the Sufian and Tasuj faults. Furthermore, the southeastern extension of the NTF connects with the North and South Bozgosh, Dozdozan, and South Sarab overturned fault zones, causing this complex to shift toward the east-northeast (Berberian and Yeats 1999). In initial published seismic reports, the NTF was categorized as a steep reverse fault (Berberian and Arshadi 1976). However, through an examination of aerial photographs, evidence of right-lateral strike-slip displacement along the fault was discovered (Berberian and Arshadi 1976). Subsequently, other researchers provided stronger evidence of the clockwise movement of waterways and other quaternary geomorphological features along the fault by documenting them on the ground (Hessami et al. 2003b, Karakhanian et al. 2004).

Berberian (1997) conducted a division of the northern fault of Tabriz into various segments, which collectively span a total distance of 210 km. The delineation of these segments was predicated upon an analysis of surface ruptures that were observed during the earthquakes that occurred in the years 1780, 1721, and 1786. According to studies conducted by Karakhanian et al. (2004), the geometry of the Tabriz fault comprises misaligned segments arranged in a right lateral manner. The estimation of earthquake risk is particularly significant for two primary sections, namely, the eastern and western regions of Tabriz city. The western area extends from Sofian to Tabriz city, and it exhibits clear indications of waterway movements, indicating the right-lateral strike-slip. Furthermore, morpho-geostructural complexities and paleoseismological investigations reveal conspicuous evidence of surface faulting caused by previous earthquakes (Hessami et al. 2003b).

The NTF exhibits a northwest sequence passing across the northern slopes of Misho Mountain, which has been identified as a distinct fault, commonly referred to as the North Misho Fault (Berberian and Yeats 1999). The fault’s eastern segment, situated north of Tabriz, originates from the southern end of its western component and stretches from Tabriz city to Basmanj village, with a tension basin lying between the two pieces. Young branches that have separated from both parts have created a depression along the western fault line. Notably, these young branches possess evident cliffs that can be traced in the northwest of Tabriz city, likely formed during the 1780 earthquake, which had a magnitude of 7.4 ( M w = 7.4 ), and previous seismic events.

3 Data

The present study employed the boundary element method to determine the slip rate on the northern fault of Tabriz. Two methods were utilized, namely, the first method based on the strain rate deduced from GPS observations, utilizing the GPS velocity field data obtained from Khorrami et al. (2019), and the second method utilizing the data of the focal mechanism of earthquakes. The earthquake focal mechanism data were collected from two dependable sources, namely, the Global Centered Moment Tensor (GCMT) and the Iranian Seismological Center (IRSC), within the longitude range of 45 – 47 ° and the latitude range of 37 – 39 ° in the study area.

Table 1 displays the focal mechanism data utilized in this investigation.

Table 1

The data gathered from two reputable sources, namely, GCMT and IRSC

Number	Origin time and location parameters						Nodal Plane 1			Nodal Plane 2			Reference
Number	Date	Time	Latitude (°N)	Longitude (°E)	Depth (km)	Mw	Strike (°)	Dip (°)	Rake (°)	Strike (°)	Dip (°)	Rake (°)	Reference
1	20050926	18:57:11.9	37.36	47.77	19.8	5.2	194	43	55	57	56	118	CMT
2	20120811	12:23:15	38.393	46.806	9.0	6.5	267	81	−175	176	85	−9	IRSC
3	20120811	22:24: 6.1	38.35	46.73	27.4	5.2	345	59	1	254	89	149	CMT
4	20120814	14:02:25	38.503	46.810	7.4	5.1	92	83	−176	1	86	−7	IRSC
5	20120815	17:49:04	38.440	46.670	4.0	4.9	80	70	165	175	76	20	IRSC
6	20120815	17:49: 8.8	38.39	46.71	13.2	5.0	246	50	133	10	56	51	CMT
7	20120816	17:14:14	38.540	46.770	10.0	4.8	256	84	165	348	76	6	IRSC
8	20120816	17:14:17.8	38.35	46.81	25.2	4.8	263	78	172	355	82	12	CMT
9	20121107	6:26:31	38.458	46.565	10.0	5.7	272	75	−173	181	84	−15	IRSC
10	20121107	06:26:33.3	38.40	46.61	15.0	5.6	183	83	7	92	83	173	CMT
11	20121223	6:38:57	38.487	44.934	12.0	5.0	70	68	149	172	62	25	IRSC
12	20130126	15:10:52.8	38.37	46.87	21.6	4.8	269	66	163	6	74	25	CMT

The table presents the earthquake data with corresponding column headers indicating the date, time, latitude, longitude, focal mechanism depth, earthquake magnitude, as well as the strike, dip, and rake parameters associated with each focal mechanism.

The earthquake focal mechanisms and GPS velocities in the study area are depicted in Figure 3. The boundary element method is the fundamental approach of this study, in which the GPS velocity field is employed to estimate the strain rate as the boundary element of the first method, while the focal mechanism data of earthquakes are used as the boundary element of the second method to estimate the stress.

Figure 3

Tectonic map of active faults in the northwestern region of Iran.

Figure 3 depicts the active faults in the northwestern region of Iran, revealing the tectonic deformation of the area through GPS velocity vectors represented by red arrows. The GPS velocity field data presented in this figure was sourced from Khorrami et al. (2019), while the earthquake focal mechanism data utilized were extracted from reliable sources such as the Global Centroid Moment Tensor (GCMT) and Iranian Seismological Center (IRSC). The fault lines shown in the illustration are indicated by black lines and were obtained from the studies by Ambraseys and Melville (2005), Berberian and Yeats (1999), and Berberian (1994).

Incorporating the NTF into the modeling required subsurface information about the fault, including its position, extension, slope, and length. To obtain this information, the active faults map of Iran (Hessami et al. 2003a) was utilized (Table 2). In addition, the locking depth of the fault was determined to be 15 km by averaging the values obtained in prior research (Djamour et al. 2011, Haji-Aghajany et al. 2019, Karimzadeh et al. 2013).

Table 2

In the simplest scenario, the NTF was characterized using average values

Section Number	X UTM ( m )	Y UTM ( m )	Strike (°)	Dip (°)	Length ( km )	Width ( km )
1	528,493	4,255,466	120	90	200	10,000

Table 2 presents the fault segment name in the first column, UTM coordinates of the fault segment’s starting point in the second and third columns, azimuth and slope of the fault in the fourth and fifth columns, and fault length and width in the sixth and seventh columns.

4 Methodology

In this investigation, the estimation of the slip rate along the NTF has been conducted employing the boundary element approach. The slip estimation was carried out by employing two distinct boundary conditions: strain boundary conditions and stress boundary conditions. The stress boundary conditions were acquired from earthquake focal mechanisms, while the strain boundary conditions were obtained from the GPS velocity field data. To perform the inversion and ascertain the slip rate, Okada’s fundamental solution was employed.

4.1 Calculation of strain boundary conditions

The deformation is characterized by the gradient of the displacement field. When three distinct points ( X ) within the initial coordinate system experience displacement via three nonparallel vectors ( u ) and are subsequently transformed into the final coordinate system ( x ) , the connection between these displacement vectors and the reference state is established by assuming a state of uniform deformation, where initially parallel lines maintain their parallelism in the final state. This relationship is mathematically formulated as follows:

(1) u i = t i + G ij X j , G ij = ∂ u i ∂ X j .

Equation (1) involves a constant value, t i , which represents the displacement of a point with respect to the origin, and G ij denotes the displacement gradients in the reference state. The Lagrangian displacement gradient tensor is denoted as G (Means 2012).

Equation (1) demonstrates that when solving the system of equations in a two-dimensional mode, six unknown parameters must be resolved. The unknowns consist of two components of the transfer vector and four components of the Lagrangian displacement gradient tensor. In three-dimensional cases, there will be a total of 12 unknowns, including three components of the transfer vector and nine components of the Lagrangian displacement gradient tensor. Each GPS station that has either known displacement or displacement speed will offer two equations in two-dimensional space and three equations in three-dimensional space. Therefore, a minimum of three points that are not collinear in two-dimensional space and four points that are not coplanar in three-dimensional space are necessary to decipher the strain tensor or displacement gradient tensor.

Linear algebra methods can be employed to solve the system of linear equations. To implement this, equation (1) must be reformulated as three matrices, whereby two matrices encompass the known quantities, and the other matrix encompasses unknown quantities. In the two-dimensional scenario, the reformulated equations for the reference state are as follows:

(2) u 1 1 u 2 1 u 1 2 u 2 2 ⋯ ⋯ u 1 n u 1 n = 1 0 0 1 1 0 0 1 ⋯ ⋯ ⋯ ⋯ 1 0 0 1 X 1 1 X 2 1 0 1 X 1 2 X 2 2 0 1 ⋯ ⋯ ⋯ ⋯ X 1 n X 2 n 0 1 0 0 X 1 1 X 2 1 0 0 X 1 2 X 2 2 ⋯ ⋯ ⋯ ⋯ 0 0 X 1 n X 2 n t 1 t 2 G 11 G 12 G 21 G 22 .

The structure of the equations in a three-dimensional mode is similar to those presented earlier, except for the inclusion of an additional index. The left-hand side of the equation will have a 3 n × 1 vector, while the right-hand side will contain a 3 n × 12 matrix and a 12 × 1 vector. Equation (2) is not restricted to only three stations in a two-dimensional mode or four stations in a three-dimensional mode, but can be applied to general cases with n stations. If there are more than three stations in a two-dimensional mode or more than four stations in a three-dimensional mode, the number of equations will exceed the number of unknown parameters. In such cases, supplementary information can be utilized to assess the accuracy of the calculated parameters.

The classical inversion theory is utilized to solve equation (2), with a particular focus on solving the linear least squares problem. The problem can be defined as follows:

(3) b → = M a → .

The problem involves a vector b → , which represents the known displacements or displacement velocities, a design matrix denoted by M that includes the initial position of the stations, and an unknown vector a → . To resolve this problem and acquire the vector a → , the product of vector b → and the inverse of matrix M is computed.

(4) a → = M − 1 b → .

If the stations are not weighted according to their distance from the point of strain tensor computation, the general linear quadratic problem can be resolved by utilizing the singular value decomposition (SVD) method (Press et al. 1992). Although there are faster methods available, SVD has an advantage in producing a more stable solution when the normal equations are close to being singular. This situation arises when the stations are nearly aligned in the two-dimensional mode or almost co-planar in the three-dimensional mode. In the event that the stations are weighted based on their distance from the strain tensor calculation point, the linear square problem can be solved with the weighted least square method, though with less stability (Menke 2018). In both scenarios, the solution comprises of (1) the unknown parameters of the model a , (2) the variances or the square of the standard deviation of the error of parameters a , which comprise the diagonal components of the covariance matrix, and (3) a statistical evaluation of the goodness of fit ( χ 2 ) (Press et al. 1992).

Once the unknowns are calculated, that is, the vector a , the displacement gradient (the last four elements in the two-dimensional mode or the last nine elements in the three-dimensional mode of vector a ) can be computed. If the reference position coordinates are known, the Lagrangian strain tensor can be expressed as follows:

(5) E ij = 1 2 [ G ij + G ji + G ki G kj ] .

4.1.1 Calculation of displacement gradient tensor using GPS velocity field

Various techniques exist for computing deformation parameters. While some methods rely on observations for parameter determination, geodesy methods utilize adjustment coordinates. To calculate displacement values along the coordinate axes, at least two observation series must be evaluated in two time epochs at each point. GPS observations can serve as an example. Several methods, including Delaunay triangulation and differential methods, are available to compute displacement gradient tensor using GPS observations. However, due to the nonuniform allocation of GPS stations, none of the aforementioned techniques provides an ideal solution.

The differential method consists of two approaches, namely, nearest points and weighted distances. Both of these approaches involve generating a uniform grid across the area under analysis and computing the velocity gradient at the central points of the grids. The key difference between the two methods lies in their approach to establish a relationship between the GPS station velocities and the analysis of central points of each grid cell.

The nearest point method involves computing gradients by selecting a fixed number of the closest stations to each central point of the grid cell. The degree of spatial sensitivity is directly associated with the GPS station density and varies in different regions of the network.

The weighted distance method involves utilizing all available stations in the network for the calculations, but with each station’s contribution weighted based on its distance from the central point of the cell. This is achieved by assigning a constant α , which indicates how the influence of a station diminishes with increased distance from the central point of the cell. The W factor is applied to each distance to adjust the weights accordingly:

(6) W = exp − d 2 2 α 2 .

Regarding the aforementioned formula, the variable d represents the distance between a GPS station and the central point of the grid. As per this formula, stations located at distances 1 α and 2 α contribute 67 and 34 % , respectively, to the adjustment calculations. Conversely, those stations located beyond a distance of 3 α have a contribution of less than 1 % . When presented in a matrix format, the W factor is depicted as a diagonal matrix and is incorporated in the linear adjustment process with the following relationship, as stated by Menke (2018):

(7) m = [ G T WG ] − 1 G T Wd .

Equation (7) involves several variables, including the design matrix G , which is a 2 n × 6 matrix positioned on the right-hand side of equation (2). The column vector of velocities on the left-hand side of equation (2) is denoted by d . In addition, the column vector of velocity and transport gradients found on the right-hand side of equation (2) is represented by m .

In spite of having more than three stations and accounting for the known errors associated with GPS velocity vectors, it is possible to estimate the accuracy of unknown values by assuming homogeneous strain between the stations. However, results obtained from calculations in grids where the absolute value of the computed unknown value is less than the 1 Σ error are not included. It is noteworthy that the calculations at the central point of each grid are performed independently of those at other nodes.

During the calculation process, a constant value of α is selected for the entire network such that the change in the desired parameter is maximized, and the number of excluded nodes where the calculated value is lower than the estimated accuracy is minimized. To achieve this, a trial-and-error process is employed to determine the appropriate α value for the entire network. It was found that an α value of 31.79 satisfied the required conditions for the entire network.

Delaunay’s triangulation method involves constructing displacement models for individual triangles using neighboring points, followed by the calculation of strain parameters at the triangles’ center of gravity. However, this method assumes that strain within each triangle is homogeneous, which is a disadvantage because the parameters of shape change vary from one point to another. Under the Delaunay triangulation approach, the strain tensor is computed at the center of gravity of each triangle.

4.2 Calculation of stress boundary conditions using earthquake focal mechanism data

Angelier (1979) proposed a method to infer the orientation of principal stresses by reversing the solution of the earthquake focal plane. However, determining the stress tensor in the Earth’s crust using only a single earthquake’s focal mechanism is not reliable (McKenzie 1969). To address this issue, a robust inversion algorithm is required to obtain the stress tensor from a significant number of earthquakes, as suggested by McKenzie (1969).

The issue of fault slip within a stress field was first introduced by Wallace (1951) and Bott (1959). Bott proposed that slip occurs on any fault plane in accordance with the maximum shear stress. Moreover, it was shown that the orientation of the shear stress is influenced by the fault plane’s orientation in the stress field and the shape ratio ( R ) , which is defined as follows:

(8) R = σ 1 − σ 2 σ 1 − σ 3 .

Regarding (8), the principal stresses axes are σ 1 , σ 2 , and σ 3 , and consistently σ 3 ≤ σ 2 ≤ σ 1 .

Carey (1974) utilized Bott’s stress criterion for inversion and added the proposition that the movements on all strike-slip faults stem from a common tensor.

Gephart and Forsyth (1984)’s stress inversion method introduced the fault instability condition, which was later utilized by Lund and Slunga (1999) to improve its effectiveness. Subsequently, in Vavryčuk (2014)’s iterative joint inversion technique, this same condition was implemented in Michael (1984)’s approach. Michael’s linear inversion approach accurately determines the principal stress direction but is associated with significant R ratio (shape ratio) estimation errors. Unlike Gephart and Forsyth (1984)’s method, Michael’s approach is linear and necessitates multiple iterations to resolve stress reversal once the fault instability condition is applied.

The sensitivity of the shape ratio ( R ) heavily relies on the accurate selection of faults, with the substitution of faults by auxiliary nodal planes leading to notable errors. Through conducting numerical tests, we demonstrate the remarkable robustness of the iterative joint inversion method in determining stress and fault orientations, resulting in significantly improved accuracy in calculating the shape ratio.

Michael utilized a method that involves utilizing the shear tension ( τ ) and normal stress ( σ n ) acting on the fault in the following manner:

(9) σ n = T i n i = τ ij n i n j ,

(10) τ N i = T i − σ n n i = τ ij n j − τ jk n j n k n i = τ kj n j ( δ ik − n i n k ) .

In this method, δ ik is the Kronecker delta operator, T represents the tension acting on the fault, n denotes the fault’s normal vector, and N signifies the unit vector direction of the shear stress acting along the fault. i , j , and k are the unit vectors in the directions x , y , and z , respectively. Equation (10) can be rewritten as follows:

(11) τ kj n j ( δ ik − n i n k ) = τ N i .

Wallace (1951)’s study demonstrates that Mohr’s stress plane graphical solution is a more efficient method for determining the intensity on multiple planes within a stress system, compared to using the formula. The graphical solution involves plotting shearing stress on the ordinate and normal stress on the abscissa, which results in arcs of circles that depict the relationship between these variables. These arcs of circles consistently represent the varying normal stress and shearing stress within the stress system. Moreover, the three principal stresses ( n 1 , n 2 , n 3 ) are depicted as normal stress on three planes perpendicular to the principal stresses. It is crucial to note that these three planes do not undergo any shearing stress. In addition, Wallace’s findings reveal that shear stress is aligned at right angles to the normal stress lines or is “down gradient” of the normal stress.

To obtain the right-hand side of equation (10), Michael (1984) utilized the assumption made by Wallace (1951) to determine the orientation of shear stress and slip direction along the fault. It was also suggested that the magnitude of shear stress was constant for all earthquakes analyzed. As this approach could not determine the absolute stress value, the variable τ was normalized to 1 in equation (11). This permitted the representation of equation (11) as a matrix in the following manner:

(12) A t = s ,

where t is the vector of stress components.

(13) t = [ τ 11 τ 12 τ 13 τ 22 τ 23 ] ,

and A is a matrix in terms of the normal vector n :

(14) A = n 1 ( n 2 2 + 2 n 2 2 ) n 2 ( 1 − 2 n 1 2 ) n 3 ( 1 − 2 n 1 2 ) n 2 ( − n 1 2 + n 3 2 ) n 1 ( 1 − 2 n 2 2 ) − 2 n 1 n 2 n 3 n 3 ( − 2 n 1 2 − n 2 2 ) − 2 n 1 n 2 n 3 n 1 ( 1 − 2 n 3 2 ) n 1 ( − n 2 2 + n 3 2 ) − 2 n 1 n 2 n 3 n 2 ( n 1 2 + 2 n 3 2 ) n 3 ( 1 − 2 n 2 2 ) n 3 ( − n 1 2 − 2 n 2 2 ) n 2 ( 1 − 2 n 3 2 ) .

In addition, the symbol s denotes the orientation of the slip vector. To ascertain the stress tensor’s five unknown components, the focal mechanism equation (14) employs the slip direction and normal vector of k earthquakes and produces 3 k linear equations. The system of equations is solved using the generalized linear inversion technique (Lay and Wallace 1995) by inserting it into equation (9).

(15) t = A − 1 s .

Equations (16)–(18) are utilized to derive the vector s in the following manner:

(16) s 1 = cos ( rake ) × cos ( strike ) + cos ( dip ) × sin ( rake ) × sin ( strike )

(17) s 2 = cos ( rake ) × sin ( strike ) − cos ( dip ) × sin ( rake ) × cos ( strike ) ,

(18) s 3 = − sin ( rake ) × sin ( dip ) .

Furthermore, the components of the vector n can be expressed in the following manner:

(19) n 1 = − sin ( dip ) × sin ( strike ) ,

(20) n 2 = sin ( dip ) × cos ( strike ) ,

(21) n 3 = − cos ( dip ) .

To determine the most optimal deviatoric stress tensor, it is crucial to assess two measures. The first measure, denoted as β , determines the level of alignment between the projected tangential tension and the fault plane slip direction on each plane. While an angle of zero is ideal, it may not always be achievable. The second measure is the magnitude of τ or τ → , which is ideally equal to one, but this may not always be the scenario.

The objective of the linear inversion method is to minimize the difference between the calculated shear stress vector and the slip vector that is observed.

In general, the slip plane is expected to exhibit higher instability. In addition, based on the stress diagram of Mohr’s circle, it is possible to identify the plane closest to Mohr’s envelope (Jaeger et al. 2009). In three dimensions, the relative magnitude of the principal stresses ( R ) can create a Mohr’s circle without requiring a scale, as affirmed by Gephart and Forsyth (1984). The spatial arrangement of the nodal planes in relation to the principal stress field determines their location on the Mohr diagram. By employing the Mohr diagram obtained from the four tensor components generated by the earthquake focal mechanism inversion process, it is feasible to determine the nodal plane with the greatest instability for all sliding friction coefficients. According to Vavryčuk et al. (2013), the instability relationship utilized in the iterative joint inversion technique from the focal mechanism of earthquakes can be expressed as follows:

(22) I = τ − μ ( σ − 1 ) μ + 1 + μ 2 ,

where:

(23) σ = n 1 2 + ( 1 − 2 R ) n 2 2 − n 3 2 ,

(24) τ = n 1 2 + ( 1 − 2 R ) 2 n 2 2 + n 3 2 − [ n 1 2 + ( 1 − 2 R ) 2 n 2 2 + n 3 2 ] 2 .

The measure of instability is communicated as a numeric quantity ranging from zero, signifying the utmost level of stability, to one, connoting the lowest level of stability.

The iterative joint inversion technique used to ascertain the focal mechanisms of earthquakes entails the initial application of Michael’s method in a standard manner, without any preconditions or knowledge concerning the orientation of fault plates. Subsequently, the principal stress direction and R ratio are determined, and their values are used to assess the instability (as per equation (22)) of the nodal plane for all inverted focal mechanisms. The fault planes are identified as those nodal planes exhibiting the highest level of instability. The fault plane orientations obtained from the initial iteration are used in subsequent iterations via Michael’s method. This process is repeated until the stress values attain an optimal level.

In the current study, the optimal stress values were achieved after six iterations. To determine fault instability using equation (22), it is crucial to know the friction coefficient value μ , which typically ranges between 0.2 and 0.8 within the fault region. However, the precise value for the NTF is unknown. In the inversion process, a typical approach is to assign an average value, such as 0.6, to the friction coefficient. Alternatively, multiple values can be used during the inversion, and the value that produces the highest instability can estimate the optimal friction coefficient for the region. Vavryčuk (2014) conducted numerical experiments that demonstrated the iterative joint inversion method for stress to be rapid, precise, and superior to the conventional linear inversion method.

4.3 Solving the inverse problem using the boundary element method

Typically, the deformation occurring at a given observation point, such as ( x , y , z ), is computed as the sum of the deformation caused by sliding on each of the planar separations (faults) and the uniform deformation field present within the studied area. The displacement vector completely characterizes the resultant deformation field.

(25) u → = u x x ˆ + u y y ˆ + u z z ˆ = u d d ˆ + u s s ˆ + u n n ˆ .

The displacement gradient tensor is expressed as follows:

(26) d u x d x d u y d x d u y d x d u x d y d u y d y d u y d y d u x d z d u y d z d u y d z .

The physical characteristics of the surrounding environment, namely, Poisson’s ratio and Young’s modulus, are considered constants. Utilizing these constants alongside the displacement gradient tensor and the governing equation of the environment, it is possible to calculate the stress, strain, and rotation tensors of a rigid body.

For every individual element, the distinct components of separation, namely, fault slip, along the three axes of the strike, dip, and normal to the fault, are as follows:

(27) D s = u s − − u s + , D d = u d − − u d + , D n = u n − − u n + .

In equation (27), the negative index positioned earlier represents the absolute displacement of the footwall, while the positive index indicates the absolute displacement of the hanging wall.

In this study, the three-dimensional boundary element method was utilized to establish the slip rate. When employing this method, there exist three significant techniques for modeling, namely, the virtual stress method, the displacement discontinuity method, and the direct boundary integral method (Crouch et al. 1983). As a fault is akin to a fracture or crack, it comprises two surfaces or boundaries, one of which is effectively aligned with the other. In such cases, traditional boundary element methods, including the direct integration method, prove ineffective in simulation. Therefore, the displacement discontinuity method was developed by Crouch (1976) to address these issues. Consequently, due to the presence of displacement discontinuity in faults, this research employed the displacement discontinuity method to model fault movement.

When using the boundary element method, faults are defined as planar rectangular discontinuities within a uniform elastic half-space, referred to as elements. The displacement or sliding of said elements can be achieved through various techniques, such as applying stress, strain, or displacement gradient tensor, implementing boundary conditions within remote regions, or through the application of displacement or stress on other elements. In addition, the utilization of composite boundary conditions is also feasible.

It is noteworthy that the discontinuity of a fault invariably pertains to the motion of the hanging wall in relation to the footwall. Hence, a negative shear dislocation indicates hanging wall movement in the positive direction of the fault extension or the X P axis of the plane coordinate system, signifying right-lateral movement. Similarly, a positive displacement in the direction of the fault dip implies the direction of reverse movement.

Based on Figure 4, the research examines the NTF by simplifying its geometry to the most basic form. It is treated as an element, and boundary conditions are established at the midpoint of this element. These conditions define three stress or displacement conditions in the azimuthal, sloping, and perpendicular directions to the element. Each direction has a corresponding boundary condition. The boundary conditions can involve stress components, absolute displacements, relative displacements, or any combination of these parameters within the element.

Figure 4

NTF geometry in the simplest possible form and as a single element.

The fault slip rate components, or relative displacement components, belonging to the element can be taken as known and integrated into the model as a boundary condition. Alternatively, these components can be treated as unknown variables and derived from the modeling process. The estimation of relative fault components within an element involves modeling that adheres to the specified initial boundary conditions and minimizes the strain energy within the model.

As for the interaction between an element containing a fault, or a part of it, with other elements, and the consequent modification of the shape of the background region, the process involves solving a series of linear equations, which are further elaborated below.

The process involves identifying a set of boundary conditions, encompassing displacement or stress, at the central point of an element.

It is important to acknowledge that the boundary conditions are solely established at the central point of an element, and they are not extended to the complete plane of the element. While a higher level of element division can lead to more precise outcomes, implementing this approach will concurrently escalate the computational time and memory storage demands.
The succeeding step involves generating a sequence of linear equations, which take the following form:
(28) τ 1 s τ 1 d σ 1 n τ 2 s τ 2 d σ 2 n = A 11 ss A 11 sd A 11 sn A 11 ds A 11 dd A 11 dn A 11 ns A 11 nd A 11 nn A 12 ss A 12 sd A 12 sn A 12 ds A 12 dd A 12 dn A 12 ns A 12 nd A 12 nn A 21 ss A 21 sd A 21 sn A 21 ds A 21 dd A 21 dn A 21 ns A 21 nd A 21 nn A 22 ss A 22 sd A 22 sn A 22 ds A 22 dd A 22 dn A 22 ns A 22 nd A 22 nn D 1 s D 1 d D 1 n D 2 s D 2 d D 2 n .

Or:
(29) τ s ( x i . y i . z i ) = ∑ j = 1 J ( A ij ss D j s + A ij sd D j d + A ij sn D j n ) + τ s b τ d ( x i . y i . z i ) = ∑ j = 1 J ( A ij ds D j s + A ij dd D j d + A ij dn D j n ) + τ d b σ n ( x i . y i . z i ) = ∑ j = 1 J ( A ij ns D j s + A ij nd D j d + A ij nn D j n ) + σ n b .

These equations consist of τ s , τ d , and σ n , which represent the stress boundary conditions along the strike, dip, and perpendicular to the fault, respectively. The coefficients A are denoted as effect coefficients or Green’s functions, which are computed using Okada’s fundamental solution (Okada 1985). Green’s functions are functions that connect the deformation field (displacement and its gradient) to discontinuities or faults in a uniform half-space (bounded by a free surface in a semi-infinite environment). The estimation of Green’s functions is subject to the availability of an analytical solution. An analytical solution is derived by considering the analytic solution of the problem and the conversion relation between the coordinate systems. By employing Green’s functions appropriate for displacement-related boundary conditions, the aforementioned system of equations can be formulated.
Upon resolving the set of linear equations, the unknown dislocation components ( D ) can be acquired.
Through the computation of the relative displacements (fault components), it is feasible to generate the analytical solution of the deformation field at any site within the environment, using the Okada analytical model.

4.3.1 Fundamental solution

To model the displacements that occur from the slip rate of faults, we employed the analytical model proposed by Okada (1985) in this investigation. The Okada model is based on the principle of dislocation theory. Initially, the deformation field that arises from a single force (point source) is modeled in this approach. Subsequently, the deformation field that arises from a rectangular source (fault plane) is formulated by integrating the point source relationships. The input parameters of the Okada model are generally classified into two categories: physical and geometrical parameters.

The physical parameters utilized in the Okada model are the Lame coefficients μ and λ of the studied region, which are estimated approximately. To this end, global average values can be implemented based on the results of the sensitivity analysis of the Okada model. The Okada model’s geometrical parameters involve the fault’s length, width, locking depth, dip, strike, separation rate or slip rate, initial fault point coordinates, and observation point coordinates. Essentially, this model transforms the fault displacement or slip rate into a displacement or velocity field resulting from it based on the fault’s geometry and the physical properties of the studied region.

5 Results

The initial approach utilized in this study involved the application of strain boundary conditions to the boundary element model, using the strain rate obtained from periodic and permanent GPS station observations, as reported by Khorrami et al. (2019). Initially, we employed the weighted distance method to compute the displacement gradient tensor in the study area, which corresponded best with all the observations in the region. This was followed by the calculation of the strain tensor via equation (5). By incorporating the strain boundary condition into the model, we were able to estimate the slipping rate in the NTF as 5.5 mm/year.

The dip-slip rate for the fault was estimated using two methods, with the first method yielding a rate of 10 − 8 mm/year, and the second method estimating a rate of 10 − 3 mm/year. Therefore, the dip-slip rate is practically zero and differs significantly from the slip rate of strike-slip. These results serve as evidence of the fault strike-slip in the northern region of Tabriz.

The outcomes of this investigation concerning the rate of slip in the NTF compared with past studies are demonstrated in Figure 5.

Figure 5

The results of this study and previous studies. The initial two lines of the presented outcomes depict the findings of this research conducted under distinct boundary conditions.

Table 3 presents the results of this study.

Table 3

Slip rate along strike slip and dip slip estimated in this study

Boundary condition	Strike-slip component (mm/year)	Dip-slip component (mm/year)
GPS velocity field	5.516	10 − 8
Focal mechanism	4.965	10 − 3

Based on the findings presented in Table 3, this study provides results for stress and strain boundary conditions. The strain boundary condition was determined through the utilization of GPS velocity field data, while the stress boundary condition was obtained using earthquake focal mechanism data. The slip rate of both the strike-slip component and the dip-slip were computed. The observed substantial value of the strike-slip component, in contrast to the negligible value of the dip-slip component, offers compelling evidence that the Tabriz fault is primarily characterized as a strike-slip fault.

Figures 6 and 7 depict the values and directions of tension and compression, which were obtained from the strain tensor analysis conducted within the studied region. The directions and values of maximum shear strain are presented in Figure 8. The results indicated that one of the planes of the maximum shear strain is in close agreement with the primary rupture of the fault, thereby confirming its strike-slip of the fault. It is noteworthy that the orientation of the principal axes of the strain tensor was considered constant throughout the entire studied region.

Figure 6

The direction and values of compression in the study region are shown by blue lines.

Figure 7

The direction and values of tension in the study region are shown by red lines.

Figure 8

The directions and values of maximum shear strain in the study region.

The second approach utilized in this study involved utilizing the stress tensor derived from earthquake focal mechanism data collected from two trusted sites, GCMT and IRSC, within the study area to determine the stress boundary condition for the boundary element model. This method involved calculating the optimal stress tensor in the study area using the iterative joint inversion method based on the focal mechanism data, which produced the best agreement with all observations within the study area. By applying the stress boundary condition to the boundary element model, the slip rate for the NTF was estimated to be 5 mm/year.

Figure 9 illustrates the principal stress axes present within the NTF region, while Figure 10 presents the associated uncertainty levels related to the calculation of these stress axes. In addition, Figure 11 depicts the Mohr’s circle diagram utilizing the method proposed by Vavryčuk (2014) for the data analyzed during this investigation within the NTF zone.

$Figure 9 The principal stress axes within the NTF. Each occurrence within the figure is represented by circles corresponding to values of σ 1 {\sigma }_{1} , a multiplication sign denoting σ 2 {\sigma }_{2} , and a plus sign representing σ 3 {\sigma }_{3} . The derivation of these stresses is indicated by the color green.$

Figure 9

The principal stress axes within the NTF. Each occurrence within the figure is represented by circles corresponding to values of σ 1 , a multiplication sign denoting σ 2 , and a plus sign representing σ 3 . The derivation of these stresses is indicated by the color green.

$Figure 10 The presented figure depicts the associated uncertainty levels in determining the principal stresses. The colors red, green, and blue indicate the confidence levels corresponding to the σ 1 {\sigma }_{1} , σ 2 {\sigma }_{2} , and σ 3 {\sigma }_{3} axes, respectively.$

Figure 10

The presented figure depicts the associated uncertainty levels in determining the principal stresses. The colors red, green, and blue indicate the confidence levels corresponding to the σ 1 , σ 2 , and σ 3 axes, respectively.

Figure 11

Mohr’s circle diagram for the data analyzed in this study.

In this study, the estimation of slip rate is accomplished by employing two distinct boundary conditions, namely stress and strain. Through the inversion process utilizing both boundary conditions, it was observed that the strike-slip component held significant magnitude, while the dip-slip component was found to be minimal and inconsequential. This inversion and subsequent research serve as concrete evidence confirming fault slippage north of Tabriz. Notably, the outcomes of this investigation exhibit the highest degree of alignment with results derived from paleoseismological methods, surpassing previous endeavors. This outcome serves as empirical validation of the suitability and accuracy of the boundary element method.

6 Discussion

The sensitivity analysis revealed that the change in displacement parameter had the most considerable impact on the outputs of Okada’s analytical model. Therefore, this study fixed all the geometric and physical parameters of the fault and the surrounding region based on the available estimations from geological and geophysical evidence. Then, by utilizing the boundary element method with Okada’s fundamental solution, the slip rate was determined in the NTF.

To perform calculations utilizing Green’s functions in the half-space environment, it is essential to image the area of study in the same environment. Since the NTF is situated entirely north of Tabriz in the UTM’s zone 38, the UTM imaging system was used to convert the environment from spherical to half-space. Furthermore, based on the image of the fault in the half-space environment, the fault’s width was assumed to be infinitely wide.

Geological measurements and geomorphological studies, including Pedrami (1987), Hessami et al. (2006), and Karakhanian et al. (2004), were implemented to analyze waterway movement, morphological indicators, and age of sediments corresponding to climatic events. These methods helped to determine the slip rate for the NTF, which ranges between 2 and 6.4 mm / year . Furthermore, Salmaniyan et al. (2023) employed a physical methodology combined with Markov chain Monte Carlo to calculate the slip rate, resulting in an estimated value of 5.3 mm / year . These estimations align with the slip rate established in this study and are lower than the current slip rates calculated through geodetic observations, such as GPS and radar interferometry, as noted by Djamour et al. (2011), Karimzadeh et al. (2013), Masson et al. (2006), and Vernant et al. (2004). By determining the slip rate via focal mechanism and geophysical data, the difference can be eliminated or lessened, and the results from geodetic methods can be brought closer to geological methods.

When analyzing the rate of landslides, it was discovered that the cause of landslides in the fault is the regional stress calculated using earthquake focal mechanism data and GPS observations, which have been utilized as boundary conditions. The inclination of the main horizontal axes toward the maximum shortening and stretching of the regional stress rate tensor along the northern fault is also responsible for the right-lateral slip rate. As a result, the image of the primary stress components along the fault is greater than the image of the corresponding component perpendicular to them.

Hessami et al. (2003b) conducted research on the NTF using paleo-seismological methods, revealing a slip rate of 1.3–6.4 mm/year. In this study, the slip rate for the same fault was estimated using two different boundary conditions, resulting in values of 5.5 and 5.0 mm/year. These estimates align with the range of slip rates obtained from paleo-seismological methods and demonstrate the superiority of the boundary element method over conventional approaches such as the block model and finite element method, both scientifically and technically. The results also indicate the seismic activity of the fault.

In the present study, the slip rate of strike-slip is determined under two distinct boundary conditions: stress and strain. The slip rate of strike-slip, utilizing the stress boundary condition, is calculated to be 5.0 mm / year , whereas the slip rate, employing the strain boundary condition, is determined to be 5.5 mm / year . Conversely, the magnitude of the slip rate for the dip-slip, considering the stress boundary condition, falls within the order of 10 − 3 mm / year , whereas the strain boundary condition yields a slip rate within the order of 10 − 8 mm / year . Remarkably, the findings of this research exhibit the highest level of agreement with the outcomes obtained through paleoseismological methods.

The modeling process in this study assumed that the Earth is a homogeneous elastic body, but utilizing a layered model that considers the heterogeneity of the Earth can result in enhanced findings (Sun and Okubo 2002, Sun et al. 1996). The effect of the curvature of the fault surface can also be investigated by using triangular elements instead of rectangular elements to determine the slip rates more accurately (Maerten et al. 2005, Marshall et al. 2008, Thomas 1993). Moreover, incorporating the effect of the earth’s gravity as a boundary condition using the model developed by Wang et al. (2006) can lead to a more realistic determination of slip rates.

7 Conclusion

This study utilized the boundary element method with Okada fundamental solution to determine the slip rate on the north fault of Tabriz. Two sets of boundary conditions were employed in modeling, namely, the strain field obtained from GPS observations and the stress field obtained from the focal mechanism data of earthquakes. The findings reveal that the NTF exhibits right-lateral strike-slip characteristics.

This study utilized a strain rate tensor obtained from GPS observations (in the first method) and a stress tensor derived from earthquake focal mechanism data (in the second method) for the entire study area, taking into account the size of the area under investigation. To enhance the findings, it is recommended to divide the study area into sections with various stress conditions and model each area separately with unique stress boundary conditions Yousefi-Bavil and Moayyed (2015).

To determine the fault slip rate, Green’s functions related to the half-space were utilized. However, it is recommended to replace the half-space Green’s functions with spherical Green’s functions to improve the deterministic slip rates by accounting for the effect of the earth’s sphericity.

Conflict of interest: There is no conflict of interest to declare.
Data availability statement: Data and Methodology used in the article are given in the text.

References

Agard, P., J. Omrani, L. Jolivet, H. Whitechurch, B. Vrielynck, W. Spakman, P. Monié, B. Meyer, and R. Wortel. 2011. Zagros orogeny: a subduction-dominated process. Geological Magazine 148(5–6), 692–25.Search in Google Scholar

Allen, M., J. Jackson, and R. Walker. 2004. Late Cenozoic reorganization of the Arabia‐Eurasia collision and the comparison of short‐term and long‐term deformation rates. Tectonics 23(2), 1–16.Search in Google Scholar

Ambraseys, N. N. and C. P. Melville. 2005. A history of Persian earthquakes. Cambridge, England: Cambridge University Press.Search in Google Scholar

Angelier, J. 1979. Determination of the mean principal directions of stresses for a given fault population. Tectonophysics 56(3–4), T17–26.Search in Google Scholar

Berberian, M. 1994. Natural hazards and the first earthquake catalogue of Iran. Volume 1: Historical hazards in Iran prior to 1900, p. 603. Tehran: International Institute of Earthquake Engineering and Seismology.Search in Google Scholar

Berberian, M. 1997. Seismic sources of the Transcaucasian historical earthquakes. Historical and prehistorical earthquakes in the Caucasus, vol. 28, p. 233–11.Search in Google Scholar

Berberian, M. and S. Arshadi. 1976. On the evidence of the youngest activity of the North Tabriz Fault and the seismicity of Tabriz city. Geological Survey Iran Report 39, 397–18.Search in Google Scholar

Berberian, M. and R. S. Yeats. 1999. Patterns of historical earthquake rupture in the Iranian Plateau. Bulletin of the Seismological Society of America 89(1), 120–39.Search in Google Scholar

Bott, M. H. P. 1959. The mechanics of oblique slip faulting. Geological Magazine 96(2), 109–17.Search in Google Scholar

Bürgmann, R., D. D. Pollard, and S. J. Martel. 1994. Slip distributions on faults: effects of stress gradients, inelastic deformation, heterogeneous host-rock stiffness, and fault interaction. Journal of Structural Geology 16(12), 1675–90.Search in Google Scholar

Carey, E. 1974. Analyse theorique et rumerique d’un modele mecanique elementaire applique a l’etude d’une population de failles. Comptes Rendus Hebdomadaires des Seances de l’Academie des Sciences, Serie D: Sciences Naturelles 279(11), 891–4.Search in Google Scholar

Crouch, S. 1976. Solution of plane elasticity problems by the displacement discontinuity method. I. Infinite body solution. International Journal for Numerical Methods in Engineering 10(2), 301–43.Search in Google Scholar

Crouch, S. L., A. M. Starfield, and F. Rizzo. 1983. Boundary element methods in solid mechanics. Journal of Applied Mechanics 50(3), 704–5.Search in Google Scholar

Dewey, J., M. Hempton, W. Kidd, F. Saroglu, and A. Şengör. 1986. Shortening of continental lithosphere: the neotectonics of Eastern Anatolia—a young collision zone. Geological Society, London, Special Publications 19(1), 1–36.Search in Google Scholar

Djamour, Y., P. Vernant, H. R. Nankali, and F. Tavakoli. 2011. NW Iran-eastern Turkey present-day kinematics: results from the Iranian permanent GPS network. Earth and Planetary Science Letters 307(1–2), 27–4.Search in Google Scholar

Gephart, J. W. and D. W. Forsyth. 1984. An improved method for determining the regional stress tensor using earthquake focal mechanism data: Application to the San Fernando earthquake sequence. Journal of Geophysical Research: Solid Earth 89(B11), 9305–20.Search in Google Scholar

Haji-Aghajany, S., B. Voosoghi, and Y. Amerian. 2019, Estimating the slip rate on the north Tabriz fault (Iran) from InSAR measurements with tropospheric correction using 3D ray tracing technique. Advances in Space Research 64(11), 2199–208.Search in Google Scholar

Hessami, K., F. Jamali, and H. Tabassi. 2003a. Major active faults of Iran. Tehran: IIEES.Search in Google Scholar

Hessami, K., F. Nilforoushan, and C. J. Talbot. 2006. Active deformation within the Zagros Mountains deduced from GPS measurements. Journal of the Geological Society 163(1), 143–8.Search in Google Scholar

Hessami, K., D. Pantosti, H. Tabassi, E. Shabanian, M. R. Abbassi, K. Feghhi, and S. Solaymani. 2003b. Paleoearthquakes and slip rates of the North Tabriz Fault, NW Iran: preliminary results. Annals of Geophysics 46, 903–15.Search in Google Scholar

Jackson, J. 1992. Partitioning of strike‐slip and convergent motion between Eurasia and Arabia in eastern Turkey and the Caucasus. Journal of Geophysical Research: Solid Earth 97(B9), 12471–9.Search in Google Scholar

Jackson, J. and D. McKenzie. 1984. Active tectonics of the Alpine—Himalayan Belt between western Turkey and Pakistan. Geophysical Journal International 77(1), 185–64.Search in Google Scholar

Jaeger, J. C., N. G. Cook, and R. Zimmerman. 2009. Fundamentals of rock mechanics. Hoboken, New Jersey: John Wiley & Sons.Search in Google Scholar

Karakhanian, A. S., V. G. Trifonov, H. Philip, A. Avagyan, K. Hssami, F. Jamali, M. S. Bayraktutan, H. Bagdassarian, S. Arakelian, and V. Davtian. 2004. Active faulting and natural hazards in Armenia, eastern Turkey and northwestern Iran. Tectonophysics 380(3–4), 189–19.Search in Google Scholar

Karimzadeh, S., Z. Cakir, B. Osmanoğlu, G. Schmalzle, M. Miyajima, R. Amiraslanzadeh, and Y. Djamour. 2013. Interseismic strain accumulation across the North Tabriz Fault (NW Iran) deduced from InSAR time series. Journal of Geodynamics 66, 53–8.Search in Google Scholar

Khodaverdian, A., H. Zafarani, and M. Rahimian. 2015. Long term fault slip rates, distributed deformation rates and forecast of seismicity in the Iranian Plateau. Tectonics 34(10), 2190–220.Search in Google Scholar

Khorrami, F., P. Vernant, F. Masson, F. Nilfouroushan, Z. Mousavi, H. Nankali, S. A. Saadat, A. Walpersdorf, S. Hosseini, and P. Tavakoli. 2019. An up-to-date crustal deformation map of Iran using integrated campaign-mode and permanent GPS velocities. Geophysical Journal International 217(2), 832–43.Search in Google Scholar

Lay, T., and T. C. Wallace. 1995. Modern global seismology. Amsterdam, Netherlands: Elsevier.Search in Google Scholar

Lund, B. and R. Slunga. 1999. Stress tensor inversion using detailed microearthquake information and stability constraints: Application to Ölfus in southwest Iceland. Journal of Geophysical Research: Solid Earth 104(B7), 14947–64.Search in Google Scholar

Maerten, F., P. Resor, D. Pollard, and L. Maerten. 2005. Inverting for slip on three-dimensional fault surfaces using angular dislocations. Bulletin of the Seismological Society of America 95(5), 1654–65.Search in Google Scholar

Marshall, S. T., M. L. Cooke, and S. E. Owen. 2008. Effects of nonplanar fault topology and mechanical interaction on fault-slip distributions in the Ventura Basin, California. Bulletin of the Seismological Society of America 98(3), 1113–27.Search in Google Scholar

Masson, F., Y. Djamour, S. Van Gorp, J. Chéry, M. Tatar, F. Tavakoli, H. Nankali, and P. Vernant. 2006, Extension in NW Iran driven by the motion of the South Caspian Basin. Earth and Planetary Science Letters 252(1–2), 180–8.Search in Google Scholar

McKenzie, D. 1972. Active tectonics of the Mediterranean region. Geophysical Journal International 30(2), 109–85.Search in Google Scholar

McKenzie, D. P. 1969. The relation between fault plane solutions for earthquakes and the directions of the principal stresses. Bulletin of the Seismological Society of America 59(2), 591–1.Search in Google Scholar

Means, W. D. 2012. Stress and strain: Basic concepts of continuum mechanics for geologists. Berlin, Germany: Springer Science & Business Media.Search in Google Scholar

Menke, W. 2018. Geophysical data analysis: Discrete inverse theory. Cambridge, Massachusetts, United States: Academic Press.Search in Google Scholar

Michael, A. J. 1984. Determination of stress from slip data: faults and folds. Journal of Geophysical Research: Solid Earth 89(B13), 11517–26.Search in Google Scholar

Mouthereau, F., O. Lacombe, and J. Vergés. 2012. Building the Zagros collisional orogen: timing, strain distribution and the dynamics of Arabia/Eurasia plate convergence. Tectonophysics 532, 27–60.Search in Google Scholar

Okada, Y. 1985. Surface deformation due to shear and tensile faults in a half-space. Bulletin of the seismological society of America 75(4), 1135–54.Search in Google Scholar

Pedrami, M. 1987. The timing, rate and nature of the late Quaternary tectonism in Iran. Bull. INQAN.Search in Google Scholar

Pollard, D. D. 1987. Theoretical displacements and stresses near fractures in rock: with applications to faults, joints, veins, dikes, and solution surfaces. Fracture Mechanics of Rock 277–49.Search in Google Scholar

Press, W. H., S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery. 1992. Numerical recipes in C. 2. Cambridge, England: Cambrige University.Search in Google Scholar

Rahmani, M., V. Nafisi, J. Asgari, and A. Nadimi. 2023. Crustal deformation in NW Iran: Insights from different invariant and variant components of geodetic strain rate tensors. ISPRS Annals of the Photogrammetry, Remote Sensing and Spatial Information Sciences 10, 631–7.Search in Google Scholar

Rizza, M., P. Vernant, J. Ritz, M. Peyret, H. Nankali, H. Nazari, Y. Djamour, R. Salamati, F. Tavakoli, and J. Chery. 2013. Morphotectonic and geodetic evidence for a constant slip-rate over the last 45 kyr along the Tabriz fault (Iran). Geophysical Journal International 193(3), 1083–94.Search in Google Scholar

Salmaniyan, M., A. Rastbood, and M. M. Hossainali, 2023, Physics-Based Approach to Deep Interseismic Creep: Implications for North Tabriz Fault Behavior Using MCMC. Engineering Proceedings 56(1), 26.Search in Google Scholar

Sun, W. and Okubo, S. 2002. Effects of earth’s spherical curvature and radial heterogeneity in dislocation studies—for a point dislocation. Geophysical Research Letters 29(12), 46-41–4.Search in Google Scholar

Sun, W., S. Okubo, and P. Vaníček. 1996. Global displacements caused by point dislocations in a realistic Earth model. Journal of Geophysical Research: Solid Earth 101(B4), 8561–77.Search in Google Scholar

Thomas, A. L. 1993. POLU3D: A three-dimensional, polygonal element, displacement discontinuity boundary element computer program with applications to fractures, faults, and cavities in the earth’s crust. Master Thesis at Stanford University.Search in Google Scholar

Vavryčuk, V. 2014. Iterative joint inversion for stress and fault orientations from focal mechanisms. Geophysical Journal International 199(1), 69–7.Search in Google Scholar

Vavryčuk, V., Bouchaala, F., and Fischer, T. 2013. High-resolution fault image from accurate locations and focal mechanisms of the 2008 swarm earthquakes in West Bohemia, Czech Republic. Tectonophysics 590, 189–95.Search in Google Scholar

Vernant, P., F. Nilforoushan, D. Hatzfeld, M. Abbassi, C. Vigny, F. Masson, H. Nankali, J. Martinod, A. Ashtiani, and R. Bayer. 2004. Present-day crustal deformation and plate kinematics in the Middle East constrained by GPS measurements in Iran and northern Oman. Geophysical Journal International 157(1), 381–98.Search in Google Scholar

Wallace, R. E. 1951. Geometry of shearing stress and relation to faulting. The Journal of Geology 59(2), 118–30.Search in Google Scholar

Wang, R., F. Lorenzo-Martín, and F. Roth. 2006. PSGRN/PSCMP—a new code for calculating co-and post-seismic deformation, geoid and gravity changes based on the viscoelastic-gravitational dislocation theory. Computers & Geosciences 32(4), 527–41.Search in Google Scholar

Yousefi-Bavil, A. and M. Moayyed. 2015. Paleo and modern stress regimes of central North Tabriz fault, Eastern Azerbaijan Province, NW Iran. Journal of Earth Science 26, 361–72.Search in Google Scholar

Received: 2023-08-21

Revised: 2023-12-27

Accepted: 2023-12-30

Published Online: 2024-03-19

This work is licensed under the Creative Commons Attribution 4.0 International License.

Articles in the same Issue

https://doi.org/10.1515/jogs-2022-0167

Keywords for this article

boundary element method; focal mechanism; GPS velocity field; North Tabriz fault; slip rate

Creative Commons

BY 4.0