Vegetation acts as a natural engineering approach that can enhance the stability of shallow soil slopes. However, the intrinsic uncertainties of vegetation, particularly the randomness of root features, have rarely been considered in vegetated slope stability analysis. To address the issue, a modified random field model is proposed to simulate the spatial variability of root parameters (e.g., root volume ratio, Rv) with different root orientations and architectures. The proposed random field model, combined with the finite element method (FEM) and limit equilibrium method (LEM), is utilised to investigate the effects of root architecture, orientation, and depth on vegetated slope reliability, considering spatial variability of soil and root properties simultaneously under rainfall infiltration. It is found that the proposed model effectively characterises the spatial variability of root parameters with different root orientations and architectures, as validated against theoretical and field-measured data. Ignoring root orientation during simulation can lead to a misestimation of the slope reliability index by up to nearly 17.1%. Moreover, without considering the spatial variability of Rv within the root zone or spatially variable soil properties (i.e., void ratio, saturated water permeability, and soil water retention curve) beyond the root zone, the vegetated slope reliability index can be overestimated by up to 4.8 times. It is recommended to adopt a triangular root system, concentrating most of the roots in the shallow depth zone and extending horizontally, which can reduce slope stability uncertainty and improve overall slope reliability.