Arrival times from an ensemble of discrete earthquakes independently contain information on hypocenter locations and jointly provide information on the velocity model. A properly formulated least squares estimation procedure can be used to determine simultaneously both hypocenter and velocity model parameters. By means of this procedure a sufficiently well distributed set of earthquakes provides information on crustal structure without some of the inherent difficulties associated with explosion refraction interpretation. Useful error estimates and resolving information are provided, and tests on artificial data indicate that successful inversions are obtainable, even for models with low-velocity zones, by using only P wave first arrivals at a number of stations. The method may have application for detecting P wave low-velocity regions resulting from dilatancy and is directly extendable to modeling shear wave velocity structure when S arrivals are used in place of P arrivals.