Desiccation crack is a prevalent natural phenomenon that plays a significant role in the stability of soil slopes. In this study, a hydromechanical coupling model incorporating a layer of stochastic cracks is developed for analyzing cracked soil slopes. To properly consider the anisotropy and spatial variability of desiccation cracks, three crack indices are generated through cross-correlated random fields via Cholesky decomposition. The seepage and mechanical behavior of a cracked slope are analyzed by adjusting stochastic parameters and rainfall conditions. Applied to the Ningzhen Mountains area in China, the model investigates the stability of slopes under various annual meteorological conditions. The results indicate that neglecting the spatial variability of cracked layer properties can lead to inaccurate assessments of instability risks at the base and water accumulation at the top of slopes. During heavy rainfall, slopes with deeper (up to 5 m) and weaker cracked layers often show a roughly planar sliding morphology. Moreover, the uncertainties in crack depth have the most pronounced influence on the uncertainties of the slope stability, more than horizontal permeability or crack aperture. The average crack aperture’s influence on slope stability depends on the relationship between crack infiltration rate and rainfall intensity.