Seismic　wave　scattering　by　multiple　structures　in　layered　half　space　is　a　classically　multiple　scattering　problem　in　infinite　domain.　The　computationally　efficient　substructure　method　is　commonly　adopted　to　solve　this　problem,　because　the　structures　can　be　modelled　by　the　finite　element　scheme　with　the　combination　of　the　high-accuracy　boundary　condition　to　provide　the　dynamic　stiffness　of　the　truncated　infinite　rock　or　soil　media.　However,　it　is　difficult　to　establish　a　high-accuracy　boundary　condition　for　multiple　scattering　problem　in　layered　half　space.　In　this　paper,　a　combined　zigzag-paraxial　boundary　condition　is　presented,　and　it　is　coupled　with　the　seismic　wave　input　method,　as　well　as　with　the　finite　element　scheme,　to　form　the　substructure　method　for　solving　the　multiple　scattering　problem　in　layered　half　space.　To　verify　the　proposed　substructure　method,　it　is　firstly　applied　to　obliquely　incident　P-SV　wave　scatterings　by　a　canyon,　a　cavity　and　an　alluvial　basin　in　homogeneous　half　space.　The　results　are　validated　by　comparing　with　several　existing　analytical　and　numerical　methods.　To　demonstrate　the　capability　of　the　proposed　substructure　method,　it　is　finally　applied　to　obliquely　incident　P-SV　wave　scattering　by　multiple　tunnels　in　layered　half　space,　and　the　effects　of　twin　tunnels,　soft　soil　interlayer　and　a　V-shaped　canyon　on　the　structural　response　are　discussed.　©　2021　Elsevier　Ltd