Suppression of surface-related and internal multiples is an outstanding challenge in seismic data processing. The former is particularly difficult in shallow water, whereas the latter is problematic for targets buried under complex, highly scattering overburdens. We have developed a two-step, amplitude- and phase-preserving, inversion-based workflow that addresses these problems. We apply robust estimation of primaries by sparse inversion (R-EPSI) to solve simultaneously for the surface-related primaries Green’s function and the source wavelet. A significant advantage of the inversion approach of the R-EPSI method is that it does not rely on an adaptive subtraction step that typically limits other demultiple methods such as surface-related multiple elimination. The resulting Green’s function is used as the input to a Marchenko equation-based approach to predict the complex interference pattern of all overburden-generated internal multiples at once. In this approach, no a priori information about the subsurface is needed. In theory, the interbed multiples can be predicted with correct amplitude and phase and, again, no adaptive filters are required. We illustrate this workflow by applying it on an Arabian Gulf field data example. It is crucial that all preprocessing steps are performed in an amplitude-preserving way to restrict any impact on the accuracy of the multiple prediction. In practice, some minor inaccuracies in the processing flow may end up as prediction errors for which corrections will be needed. Hence, we conclude that the use of conservative adaptive filters were necessary to obtain the best results after interbed multiple removal. The obtained results indicate promising suppression of surface-related and interbed multiples.