In order to understand physical processes of megathrust earthquakes, and thus to analyze hazards of the regions, proper information of the rupture processes is essential. In this study, we develop a technique for determination of large spatial-scale rupture processes on megathrust faults in subduction zones as soon as possible, based on teleseismic broadband full-waveform inversions using realistic models of finite fault geometries and 3D global velocity structures. Finite fault geometry models are constructed for eleven subduction zones of the circum-Pacific area, utilizing the 3D plate interface model Slab 1.0 of USGS. Green’s functions for available pairs of subfaults and globally distributed stations to be used in the inversions are computed in the global 3D velocity model using spectral element method rather than in a 1D velocity model. This scheme accommodates synthetic waves with lateral velocity variations especially for surface waves, which is useful in estimations of large spatial-scale slip distributions. The data utilized are teleseismic broadband full-waveform records including body waves and surface waves at the epicentral distance of 30o ~ 80o. The observed data are inverted for fault rupture processes, based on the non-negative least-square technique with a multiple source time window scheme for each subfault. The applicability of the method is tested for the 2011 Mw 9.0 Tohoku-Oki earthquake, and the result shows that the proposed scheme of the finite fault inversion produces reasonable results of slip distributions.