The sea-surface interface between ocean and air is time varying and can be spatially rough as a result of wind, tides, and currents; the shape of this interface changes over time under the influence of wind, tides, etc. As a result, waves impinging on the sea surface are continuously scattered. In the case of marine seismic, the multiple-scattered waves propagate downward into the underwater formation and result in complex seismic responses. To understand the structure of the responses, we have adopted a multistage algorithm for computing the scattered waves at the sea surface. Specifically, we first extrapolate the upgoing incident waves stepwise using thin-slab approximation from the scattering theory based on the De Wolf approximation of the Lippmann-Schwinger equation. Then, we implement the air-water boundary condition at the sea surface. Finally, we use the irregular boundary processing technique to compute the time-varying undulating sea-surface scattered waves from different scattering stages. To overcome the angular limitation of the original parabolic approximation, we introduce a multidirectional parabolic approximation based on computational electromagnetics. Numerical tests indicate that the multistage algorithm presented here can accurately calculate sea-surface scattered waves and should be useful in investigating the structure of marine seismic responses.