We present a seismic waveform inversion methodology based on the Gauss-Newton method from pre-stack seismic data. The inversion employs a staggered-grid finite difference solution of the 2-D elastic wave equation in the time domain, allowing accurate simulation of all possible waves in elastic media. The partial derivatives for the Gauss-Newton method are obtained from the differential equation of the wave equation in terms of model parameters. The resulting wave equation and virtual sources from the reciprocity principle allow us to apply the Gauss-Newton method to seismic waveform inversion. The partial derivative wavefields are explicitly computed by convolution of forward wavefields propagated from each source with reciprocal wavefields from each receiver. The Gauss-Newton method for seismic waveform inversion was proposed in the 1980s but has rarely been studied. Extensive computational and memory requirements have been principal difficulties which are addressed in this work. We used different sizes of grids for the inversion, temporal windowing, approximation of virtual sources, and parallelizing computations. With numerical experiments, we show that the Gauss-Newton method has significantly higher resolving power and convergence rate over the gradient method, and demonstrate potential applications to real seismic data.