Traveltime data are used to determine upper and lower bounds on velocity variations in the earth by an iterative method. In effect, the range of models is found which is consistent with the data, rather than a single model which is “best fitting” in some sense. The algorithm used, a variant of the row‐action algorithms commonly applied to tomographic inversions, requires little in‐core memory and has proven feasible for data sets of the order of hundreds of thousands to a million traveltimes. Any inequality constraint, such as that the velocity is always positive, may be incorporated into the formulation. Data errors can be included both locally (as strict constraints on each datum) and globally (as a constraint on a total error measure). The method may also be used to derive the velocity structure which results in the minimum 𝓁1 norm of the residual misfit. Data from the Grimsel crosshole experiment are used to map confidence bounds of 0.02 ms on 1521 traveltime residuals into upper and lower bounds on seismic velocities. There is great variation in the widths of these bounds as a function of subsurface position from 0.1 to 0.9 km/s. The distributions of the bounds agree with the parameter resolution values found from a singular value decomposition (SVD) and suggest that a low‐velocity mylonitic zone, seen in tomographic inversions of the traveltime data, is adequately imaged. Though the data were corrected for seismic anisotropy, significant alternating positive and negative velocity perturbations in poorly constrained quadrants of the crosshole region suggest that some residual anisotropy effects are still present.