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 considering 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 propose a multistage algorithm for computing the scattered waves at the sea surface. Specifically, we first extrapolate the upgoing incident waves stepwise using the thin-slab approximation from the scattering theory based on the De Wolf approximation of the LippmannSchwinger 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 multi-directional parabolic approximation based on computational electromagnetics. Numerical tests show that the multistage algorithm presented here can accurately calculate the sea surface scattered waves and should be useful in investigating the structure of marine seismic responses.