Abstract Slip rate is a critical parameter for describing geologic and earthquake rates of known active faults. Although faults are inherently three-dimensional surfaces, the paucity of data allows for estimating only the slip rate at the ground surface and often only few values for an entire fault. These values are frequently assumed as proxies or as some average of slip rate at depth. Evidence of geological offset and single earthquake displacement, as well as mechanical requirements, show that fault slip varies significantly with depth. Slip rate should thus vary in a presumably similar way, yet these variations are rarely considered. In this work, we tackle the determination of slip rate depth distributions by applying the finite element method on a 2D vertical section, with stratification and faults, across the central Apennines, Italy. In a first step, we perform a plane-stress analysis assuming visco-elasto-plastic rheology and then search throughout a large range of values to minimize the RMS deviation between the model and the interseismic GPS velocities. Using a parametric analysis, we assess the accuracy of the best model and the sensitivity of its parameters. In a second step, we unlock the faults and let the model simulate 10 kyr of deformation to estimate the fault long-term slip rates. The overall average slip rate at depth is approximately 1.1 mm/yr for normal faults and 0.2 mm/yr for thrust faults. A maximum value of about 2 mm/yr characterizes the Avezzano fault that caused the 1915, Mw 7.0 earthquake. The slip rate depth distribution varies significantly from fault to fault and even between neighbouring faults, with maxima and minima located at different depths. We found uniform distributions only occasionally. We suggest that these findings can strongly influence the forecasting of cumulative earthquake depth distributions based on long-term fault slip rates.