A 3D non-homogenous finite element model (FEM) is developed to investigate the spatial variations of interseismic deformation for the Kazerun Fault System (KFS) in the Zagros Mountains of Iran. The model includes 19 fault segments that were extracted from geological maps and previous studies, and the average slips in the dip and strike directions on these segments were computed. The contemporary surface deformation is simulated using a free horizontal detachment surface. The dip angles of the faults in the model are varied at 90°, 70°, 50° and 30° to simulate different 3D representations of the fault systems. Tectonic loading at the boundaries of the region is applied using predicted GPS velocity vectors to the north (southern part of the Central Iran Block) and south (southern region of the Zagros mountain belt), which were obtained by solving inverse and forward problems. Where possible, the fault slip rates that are obtained using our non-homogeneous finite element model are validated using the long-term geologic and instantaneous GPS slip rates. The model is then used to estimate the dip- and strike-slip rates of the fault segments of the KFS for which no a priori information was available. We derive an upper bound of 1 mm/year for the average dip-slip rate in the region, which is consistent with estimates from geomorphologic observations. The modeling results show that in addition to the 4 main faults (Dena, Kazerun, Kareh Bas and Main Recent), other faults, such as the Zagros Front, Main Front, High Zagros and Mishan faults, accommodate up to 2.5 mm/year of the differential movement between the North and Central Zagros. We also investigated the contrast in rigidity between the southern and northern areas of the Zagros mountain belt and found that a rigidity contrast of 2 best explains the GPS data of contemporary surface deformation. Neglecting to account for the rigidity contrast in the model can lead to biased estimates of the fault slip rate of up to 50 %.