This article reports on the development of a modular two-step, finite-element methodology for modeling earthquake ground motion in highly heterogeneous localized regions with large contrasts in wavelengths. We target complex geological structures such as sedimentary basins and ridges that are some distance away from the earthquake source. We overcome the problem of multiple physical scales by subdividing the original problem into two simpler ones. The first is an auxiliary problem that simulates the earthquake source and propagation path effects with a model that encompasses the source and a background structure from which the localized feature has been removed. The second problem models local site effects. Its input is a set of equivalent localized forces derived from the first step. These forces act only within a single layer of elements adjacent to the interface between the exterior region and the geological feature of interest. This enables us to reduce the domain size in the second step. If the background subsurface structure is simple, one can replace the finite-element method in the first step with an alternative efficient method. The methodology is illustrated in a companion paper (Yoshimura et al., 2003) for several 3D problems of increasing physical and computational complexity. We consider first a flat-layered, stratigraphic system. For this simple case, the first step can be carried out by means of 3D Green's function evaluations. The extension to more general problems is illustrated by two examples: a basin and a hill, with the same background stratigraphy. To verify the two-step procedure with a problem for which the finite-element method is used throughout, we model ground motion in a small region of the Los Angeles Basin, using both the two-step domain-reduction method and the traditional approach in which the computational domain contains both the source and the geological region of interest.