Fault slip rate distribution plays an important role in earthquake studies. Faults are loaded at very slow rates in continental interiors. So, interaction among faults and resulting slip distribution can give rise to earthquakes on other faults after a long period of quiescence and seismicity that can migrate from one fault to the onother one. NW Iran-Eastern Turkey is a region of active deformation as a result of oblique collision of Arabia-Eurasia tectonic plates. In northwestof Iran, deformation between the Central Iranian block and the Caucasus domain is accommodated by a fault system and mainly by right lateral strike-slip on the North Tabriz fault. In the current study, we did slip rate partitioning in the fault system of northwest Iranian plateau using the concepts of dislocation theory. Modelling approach is described by Gomberg and Ellis (1994), Flerit et al., (2003) and Armijo et al., (2004) and it differs from rigid block models (Reilinger et al., 2006; Djamour et al., 2011) in which dislocation conditions at the boundaries of blocks are often incompatible with geological evidences. In the alternative method of Flerit et al., (2003), slip everywhere has a direction of motion consistent with geological constraints. The dislocations do not divide the region into closed rigid blocks and slip can vary along strike as observed geologically. Finally we obtain a tectonic model for NW Iran-Eastern Turkey that is more realistic than rigid block model (Reilinger et al., 2006; Djamour et al., 2011) or models based on seismic or geologic strain rates (Haines, 1982; Haines and Holt, 1993; Jackson et al., 1995; Masson et al., 2005). For this purpose we use a three dimensional boundary elements method. First, we consider an elastic and homogeneous half-space for the study area. Then geometric data of fault system are collected from geological and geophysical sources including fault length, width, dip, and locking depth. For Lame coefficients, we use global average values. Both mentioned geometrical and physical data are kept fixed in the modeling process. Then, strain tensor that best fits the GPS data is estimated for the study area using least squares method. Then, stress rate tensor is estimated using generalized Hook's law. Geomerical chracteristics of faults, physical characteristics of crust and stress rate tensor act as boundary conditions in the model. Faults are locked in normal direction but they are allowed to slip freely in strike and dip directions under the influence of boundary conditions. Regarding the strike changes of faults, the fault surfaces are divided by different segments in strike direction with constant strikes and dips. Then fault segment surfaces are divided into 1km elements. Finally, we have free slipping elements in strike and dip directions as inputs for modeling. Our model is fitted to the fault traces data set of NW Iran-eastern Turkey. The results indicate the dependency of the partitioned slip rate on the boundary conditions and confirm the existence of interaction among faults. Also, partitioned slip rates show that the Chalderan, Guilato-Siahcheshmeh-Khoy, Nakhchivan, North Tabriz and Pambak-Sevan-Sunik faults are right-lateral strike slip in all cases. Also, the slip rate in these faults is almost symmetric and reaches its maximum value around the center of the faults. We show that the maximum value of slip rate in the fault plane is reduced by partitioning, which it will be definitely closer to reality. According to the gridding for slip rate partitioning in the fault system, the highest value of slip rate is always related to the North Tabriz Fault. Previous studies show that the geological slip rate estimates are lower than the present-day GPS-derived slip-rates along the North Tabriz fault. We show that slip rate partitioning solves this discrepancy by considering the mechanical interaction among faults. Our partitioned slip rates for North Tabriz Fault are lower than geodetic rates and are more consistent with geological rates. Finally, we present a model that fits best with the geological constraints. The proximity of the partitioned slip rate to the paleo-seismic values indicates the closeness of the partitioning results to reality with the Boundary Elements Method, compared to other analytical and numerical methods. This research may open new research direction to handle the differene between geologic and geodetic slip rates values in the Iranian Plateau. The boundary elements method is both faster and more accurate for modeling compared to the finite element method used by Khodaverdian et al. (2015). Considering the effect of topography and sphericity of Earth, using the Galerkin boundary element method developed by Thompson (2019) is proposed to get more realistic results. The coefficients matrix in the of Boundary Elements Method is fully populated and in high dimensions it takes a lot of time to solve the resulting system of equations. Sparsing of the coefficient matrix using wavelet transforms is suggested (Ebrahimnejad et al., 2010) in this study. The use of iterative computational methods along with parallel processing will also reduce the computational time (Thompson and Meade, 2019). [ABSTRACT FROM AUTHOR]