The　spatial　variability　of　soil　properties　is　pervasive,　and　can　affect　the　propagation　of　seismic　waves　and　the　dynamic　responses　of　soil-structure　interaction　(SSI)　systems.　This　uncertainty　is　likely　to　increase　the　damage　state　of　a　structure　and　its　risk　of　collapse.　Additionally,　conducting　multiscale　simulations　efficiently　in　the　presence　of　uncertainties　is　a　pressing　concern　that　must　be　addressed.　In　this　work,　a　3D　probabilistic　analysis　framework　for　an　SSI　system　considering　site　effects　and　spatial　variability　of　soil　property　(i.e.,　elastic　modulus,　E)　has　been　proposed.　This　framework　is　based　on　the　random　finite　element　method　(RFEM)　and　domain　reduction　method　(DRM).　A　multiscale　model　of　a　five-story　reinforced　concrete　(RC)　frame　structure　was　developed　on　an　ideal　3D　slope　to　verify　the　effectiveness　of　the　proposed　framework.　The　dynamic　responses　of　the　structure　were　analyzed,　and　the　peak　floor　acceleration　(PFA)　and　peak　interstory　drift　ratio　(PSDR)　were　selected　to　estimate　the　damage　state　of　structures.　It　was　found　that　the　proposed　method　significantly　improves　computational　efficiency　approximately　20　times　compared　with　the　direct　method.　In　the　regional　models,　with　the　increase　of　the　coefficient　of　variation　(COV)　of　E,　the　energy　of　seismic　waves　becomes　more　concentrated　at　the　crest　and　the　response　spectrum　value　of　medium　and　long　periods　increases.　In　the　local　SSI　model,　the　soil　variability　increases　the　mean　value　of　PSDR,　resulting　in　a　more　severe　damage　state　compared　to　the　results　from　the　deterministic　analysis.　Consequently,　this　study　provides　some　suggestions　for　engineering　practice,　and　the　importance　of　probabilistic　analysis　considering　spatially　variable　soils　in　the　SSI　problem　is　highlighted.　©　2024　American　Society　of　Civil　Engineers.