We present a novel model for the true vertical depth (TVD) positioning error in horizontal wells. The model utilizes a non-stationary Gaussian process known as the integrated Ornstein-Uhlenbeck process. This process is continuous and exhibits the known systematic accumulation of vertical errors with increasing measured depth. The well corrections produced are smooth and maintain the underlying shape of the well path. The smoothness can be adjusted through a correlation range parameter. Using the geological constraints contained in the zonation, we can predict or simulate surface depths and well depths simultaneously. The size of the surface and well displacements are governed by the relative magnitude of the surface and well path TVD uncertainties. The resulting well paths and surfaces stay within their uncertainty envelopes and are consistent with the zonation. The surface/well relationships can be expressed as a highly dimensional truncated multivariate Gaussian distribution. We draw samples from this distribution using an efficient rejection sampling strategy that allows fields with hundreds of horizontal wells to be handled.