Finite element programs allow for an enhancement in the computational capabilities of earthquake site response models. However, many finite element programs involve complicated constitutive models that are difficult for the end user to implement. We present a methodology for modeling earthquake site response within a general finite element framework, using an overlay model to represent nonlinear soil behavior. Using parallel load-carrying elements with varying stiffness and yield stress, the behavior of any given backbone stress-strain relation can be replicated, along with hysteretic unloading-reloading (Masing) behavior. Our finite element modeling methodology makes use of existing conventional elastoplastic material models available in any general finite element program, without requiring the specification of any complicated constitutive models. To represent overlay elements in a finite element model, the user defines a number of finite elements and assigns each of them identical node numbers. The only necessary input parameters are density, elastic constants, and the backbone curve of a stress-strain relation. This paper focuses on the application of one-dimensional total-stress site response, but the framework could be easily extended to model cyclic hardening and softening, three-dimensional wave propagation, and soil-structure interaction.