The purpose of this study is to develop and test a procedure for simulating acceleration time histories of large subduction earthquakes. The ground motions of the large event are obtained by summing contributions from fault elements to simulate the propagation of rupture over the fault surface. The radiation from the fault elements is represented by empirical source functions derived from near-source strong motion recordings of Mw = 7 aftershocks, and gross aspects of wave propagation are modeled by calculated Green's functions. The procedure has been tested against the recorded strong ground motions of the Mw = 8.0 Michoacán, Mexico, and Valparaiso, Chile, earthquakes of 1985. We find that models of heterogeneous slip in these events derived by other investigators from the analysis of teleseismic and near-source velocity seismograms also explain the shorter period motions of the recorded accelerograms. The response spectra of the simulated motions generally have little or no significant bias in the period range of 0.05 to 2 sec for both earthquakes. Also, the peak accelerations, durations, and envelope shapes of the time histories are in good agreement with the recorded motions. The uncertainty associated with the modeling procedure, estimated from the misfit between recorded and simulated response spectra, shows a dependence on period and site conditions. This uncertainty estimate includes the uncertainty involved in using empirical source functions derived from one subduction zone to simulate strong ground motions in another subduction zone. The procedure is applied in a companion paper (Cohee et al., 1991) to estimate strong ground motion characteristics in the Pacific Northwest region of the United States from hypothesized Mw = 8 subduction earthquakes on the Cascadia plate interface.