As the primary structure-borne path, bearings play a critical role in gearbox vibration and sound radiation. To analyze the contact between rolling elements and races, and predict the bearing stiffness more accurately, detailed models of rotor-bearing systems are required. In this study, a system level finite element/contact mechanics model of the axle system is proposed, which includes the detailed geometry of hypoid gear pair, shafts, tapered roller bearings and axle housing. A numerical method is applied to calculate the tapered roller bearing stiffness matrix. The proposed method is then compared against the traditional analytical method. A series of parametric studies are also performed to investigate the effects of bearing preloads on the transmission error, bearing stiffness, and bearing contact pattern. The proposed approach demonstrates the capability to calculate fully populated bearing stiffness matrix and time-dependent contact characteristics between rollers and raceways, and thus can be employed in the axle system design and further dynamic analysis.