Landslides significantly threaten lives and infrastructure, especially in seismically active regions. This study conducts a probabilistic analysis of seismic landslide runout behavior, leveraging a large-deformation finite-element (LDFE) model that accounts for the three-dimensional (3D) spatial variability and cross-correlation in soil strength — a reflection of natural soils' inherent properties. LDFE model results are validated by comparing them against previous studies, followed by an examination of the effects of univariable, uncorrelated bivariable, and cross-correlated bivariable random fields on landslide runout behavior. The study's findings reveal that integrating variability in both friction angle and cohesion within uncorrelated bivariable random fields markedly influences runout distances when compared with univariable random fields. Moreover, the cross-correlation of soil cohesion and friction angle dramatically affects runout behavior, with positive correlations enlarging and negative correlations reducing runout distances. Transitioning from two-dimensional (2D) to 3D analyses, a more realistic representation of sliding surface, landslide velocity, runout distance and final deposit morphology is achieved. The study highlights that 2D random analyses substantially underestimate the mean value and overestimate the variability of runout distance, underscoring the importance of 3D modeling in accurately predicting landslide behavior. Overall, this work emphasizes the essential role of understanding 3D cross-correlation in soil strength for landslide hazard assessment and mitigation strategies.