In this paper the boundary element algorithm which uses the time-convoluted traction kernels is applied to a numerical parametric study on the seismic behavior of three-dimensional Gaussian-shaped hills subjected to vertically propagating incident waves. All calculations were executed in the time-domain and the medium was assumed to have a linear elastic constitutive behavior. Results are discussed in both time and frequency domain with respect to the dimensionless parameters. It was shown that wave length and site geometry, including shape and dimension ratios and, to some extent, wave type are the key independent parameters governing hill amplification behavior. Comparing two- and three-dimensional hills with similar shape ratios, two-dimensional hill had greater characteristic periods, where the three-dimensional hill had greater maximum amplification potential. Three-dimensionality has a strong effect on the seismic responses of the hill; however the rate of seismic response variation with the three-dimensionality factor depends on the shape ratio. It was shown that two-dimensional behavior was dominant in low height three-dimensional hills, however, as the shape ratio increased, three-dimensionality effects appeared and the seismic response of the hill tends toward the axisymmetric three-dimensional hill.