Joint inversion of multiple geophysical data-sets is promising to reduce uncertainties in independently inverted models. Here, we present an iterative joint inversion approach for P and S traveltime data using cross-gradient function as constraint term. This type of joint inversion scheme links independent inversions through iterations and the cross-gradient function. The primary advantage of this joint inversion strategy is to avoid determining relative weighting of different data-sets. To investigate the performance of this method, we test our algorithm in synthetic examples of P and S traveltime data and field data acquired in west Texas. The results of synthetic example show that the joint inversion significantly reduces the ambiguities of inverted models and improves the identification of boundaries. In results of field data, jointly inverted S velocities have better correlation with P velocities. Moreover, lithologies delineated from Vp/Vs map by joint inversion matches log data very well and also shows clearly a dipping structure below reservoir that was not shown in previous crosswell tomography results.