To tackle the arising problem of existing weak shale beddings in the slope stability analysis in the oil sand deposits, it’s important to understand the stress-strain behavior of such weak shale beddings in different formats of layering. Stimulation of 3D sand/shale formation is set out using the triaxial test. Peak strength, friction angle, and the plastic strain of the specimen are validated using a simple setting numerical model with the control of the mesh effect for the optimum computational time. Pygeostat will then be introduced to generate heterogeneous oil sand realizations mimicking HIS in the McMurray formation with four carious shale configurations, including the volume fractions of shale, the horizontal and vertical range, and the inclination angle of shale beddings. Different combinations of volume fraction of shale will influence the stress in both the x and y directions. The following machine learning models are proposed to predict the upscaling geomechanics parameters of oil sands. Linear Regression Model of Ridge and Polynomial features, the Artificial neural network (ANN) model and Random Forest will assist in the ease of the calculation process.