A boundary integral formulation is presented and applied to model the ground motion on alluvial valleys under incident P, S and Rayleigh waves. It is based on integral representations for the diffracted and the refracted elastic waves using single‐layer boundary sources. This approach is called indirect BEM in the literature as the sources' strengths should be obtained as an intermediate step. Boundary conditions lead to a system of integral equations for boundary sources. A discretization scheme based on the numerical and analytical integration of exact Green's functions for displacements and tractions is used. Various examples are given for two‐dimensional problems of diffraction of elastic waves by soft elastic inclusion models of alluvial deposits in an elastic half‐space. Results are displayed in both frequency and time domains. These results show the significant influence of locally generated surface waves in seismic response and suggest approximations of practical interest. For shallow alluvial valleys the response and its resonant frequencies are controlled by a coupling mechanism that involves both the simple one‐dimensional shear beam model and the propagation of surface waves.