In engineering practice, a two-dimensional (2D) geotechnical cross-section is often used to assess soil liquefaction potential and its consequences on civil structures in a specific site for a given earthquake scenario. However, selection of a representative 2D cross-section for a real site is a challenging task because subsurface soils often exhibit spatial variability (e.g., spatially varying soil stratigraphy and engineering properties) in a three-dimensional (3D) subsurface space. Therefore, liquefaction assessment results interpreted from selected 2D cross-sections for representation of a 3D subsurface space may be unconservative, or even biased, leading to significant risk to civil infrastructure. Furthermore, assessment of 3D subsurface models demands higher computational effort than 2D cross-sections, presenting an additional obstacle. To tackle these issues, this study develops a novel and computationally efficient method for 3D liquefaction assessment with appropriate consideration of soil stratigraphy and property spatial variability within a 3D subsurface space. The proposed method integrates the cone penetration test (CPT)-based simplified liquefaction assessment methods with a novel 3D soil stratigraphy and engineering properties modelling method to characterize the spatially varying soil liquefaction potential and liquefaction-induced settlement for a given earthquake scenario in 3D. The proposed 3D method is illustrated using exploration data from Christchurch, New Zealand. An illustrative example indicates that the proposed method can properly assess significantly non-uniform and heterogeneous soils, and evaluate liquefaction potential and liquefaction-induced settlement of soils within the selected methodologies in a 3D subsurface space.