Abstract: A new indirect boundary integral equation method for numerically solving the scattering of seismic waves by a two-dimensional alluvial valley of arbitrary shape in poroelastic layered half-space is proposed based on the single layer potential approach. The Green's functions of compressional and shear wave sources in poroelastic layered half-space are derived firstly, then the scattered waves in alluvial valley and layered half-space are constructed by using the fictitious wave sources close to the interface between the valley and the half-space respectively, and the magnitude of the fictitious wave sources is determined through continuous boundary conditions, and the total response can be obtained by the superposition of the free field and the scattered field. The precision of the proposed method is verified by the satisfaction extent of boundary conditions and the comparison between the degenerated solutions of single-phased half-space and the well-known solutions. The diffraction of plane SV waves around an alluvial valley in poroelastic layered half-space is studied with a typical example. It is illustrated that there exists significant difference between the scattering of waves around the poroelastic alluvial valley and the dry soil valley, and the drainage condition has large impact on earthquake ground motion in the valley. The wave scattering in layered sites and that in homogenous half-space are different in nature, and the total response characteristics combine the resonant behavior of the soil layers and the characteristics of scattering waves around the alluvial valley. Hence, to simulate the propagation of seismic waves more accurately, it is necessary to consider the multi-phased and layered characteristics of soil medium in reality.