<jats:p> Waves propagating across a vertical seismic profiling (VSP) array may be distinguished by their differing arrival times and linear-moveout velocities. Current methods typically assume that the waves propagate uniformly with an unvarying wavelet shape and amplitude. These assumptions break down in the presence of irregular spatial sampling, event truncations, wavelet variations, and noise. I present a new method that allows each event to independently vary in its amplitude and arrival time as it propagates across the array. The method uses an iterative global nonlinear optimization scheme that consists of several least-squares and two eigenvalue problems at each step. Events are stripped from the data one at a time. As stronger events are predicted and removed, weaker events then become visible and can be modeled in turn. As each new event is approximately modeled, the fit for all previously removed events is then revisited and updated. Iterations continue until no remaining coherent events can be distinguished. As VSP data sets are typically not large, the expense of this method is not a significant limitation. I demonstrate with a real-data example that this iterative approach can lead to a significantly better VSP wavefield separation than that which has been available when using conventional techniques. </jats:p>