Northeastern part of Iran islocated in the Kopeh Dagh seismotectonic province. This province is part of the Alpine-Himalayan orogenic belt. In the last two decades, a number of earthquakes occurred in this region. Therefore, investigating and relocating these earthquakes with adequate precise play a great role in seismicity research of the study area. Adequate crustal structure velocity has a key role in relocating the earthquakes. Since the earthquakes in this region are shallow earthquakes, the upper part of crustal structure can be acquired. So in this study, Mashhad and Quchan seismological networks data were used to determine the crustal structure velocity. In order to determine the velocity structure, V_p/V_s was first determined using the methods of Wadati and travel time of1.72 .Then by applying various filters and selection of several earthquakes 1D velocity structure of the arrival times was determined. Results reveal two velocity discontinuities approximately at 4 and 14 km depths. The resulting model consists of a 4-km-thick layer with V_P = 5.9 kms⁻¹ overlying a 10-km-thick layer with V_P = 6.17 kms⁻¹. Then the initial 1D velocity model was used to determine the three-dimensional velocity structure using LOTOS code. The results show significant changes in the vicinity of Mashhad. Two zones with higher velocity are observed in the north and south of Mashhad due to Kashaf Roud and Binaloud Faults. In the north of Mashhad, a change in velocity is observed, which might be due to the Tous fault. In the south of Mashhad, a sharp velocity discontinuity at the depth of 4 km is probably due to an alluvium cover above a bedrock of varying thickness. Also, the tomography results indicate that the area of Shandiz fault system in the south of Mashhad is associated with a lower velocity zone compared to its adjacent areas.