We have proposed a two-stage model-based interferometric interpolation method for filling in gaps in marine seismic data. The first stage is creating virtual traces and the second stage is utilizing modified virtual traces for interpolation. There are three steps for creating a virtual trace. First, set the virtual trace's source and receiver coordinates and retrieve two common-receiver gathers (CRG) from existing data. At least one CRG should be obtained, otherwise, our method cannot be performed. Second, if only one CRG is obtained, then a synthetic CRG should be generated using known velocity model (e.g., water-layer model). Third, the virtual trace is created by summing correlated trace pairs over all identical source positions in two CRGs. When all virtual traces are created, a least-squares matching filter or a nonstationary prediction error filter (PEF) should be used together with virtual traces to do the interpolation. If the trace gap is small or the subsurface geology is simple, the local matching filter is preferred for the improvement of signal-to-noise ratio, otherwise, we need to estimate a PEF from virtual traces and then use this PEF for interpolation. A simple velocity model and Sigsbee2b velocity model are used to validate the effectiveness of our method. Results show that the trace gaps can be kinematically interpolated, but there are still artifacts that are due to the approximations made in our theory.