The wavefield is typically simulated for seismic exploration applications through solving the wave equation for a specific seismic source location. The direct relation between the form (or shape) of the wavefield and the source location can provide insights useful for velocity estimation and interpolation. As a result, I derive partial differential equations that relate changes in the wavefield shape to perturbations in the source location, especially along the Earth's surface. These partial differential equations have the same structure as the wave equation with a source function that depends on the background (original source) wavefield. The similarity in form implies that we can use familiar numerical methods to solve the perturbation equations, including finite difference and downward continuation. In fact, we can use the same Green's function to solve the wave equation and its source perturbations by simply incorporating source functions derived from the background field. The solutions of the perturbation equations represent the coefficients of a Taylor's series type expansion of the wavefield as a function of source location. As a result, we can speed up the wavefield calculation as we approximate the wavefield shape for sources in the vicinity of the original source. The new formula introduces changes to the background wavefield only in the presence of lateral velocity variation or in general terms velocity variations in the perturbation direction. The approach is demonstrated on the smoothed Marmousi model.