This article describes an interdisciplinary methodology to calculate the probability of failure for bearing axial cracking, the dominant failure mode in the intermediate and high-speed stages of many wind turbine gearboxes. This approach is mainly a physics-domain method with needed inputs from the data domain. The gearbox and bearing design along with operations data and component failure records from a wind power plant provide the input to physics-based models and define axial cracking damage metrics. Furthermore, the physics-domain models predict the bearing loads and sliding velocities, which are the essential elements for quantifying the accumulated frictional energy. Both accumulated frictional energy and electrical energy generation are proposed as damage metrics for bearing axial cracking. A first-order reliability method is then used to compare the proposed damage metrics to failure threshold functions and calculate the probability of failure of each individual bearing. Although the probability of failure for the failed turbines is not separated from the population, a feature engineering analysis shows the potential of frictional energy as a damage metric when combined with roller loads, bearing sliding speed, lubricant type, and terrain features. Through statistical analysis of historical data, the proposed methodology enables reliability assessment of axial cracking in individual wind turbine bearings and connects the reliability forecast with turbine design and operations.