In marine seismic experiments it is important to model the reflection from the sea surface. The sea surface is normally assumed to be flat, so the reflection coefficient is (minus) unity and the time delay is constant. In fact the sea surface is rough, causing three effects: the time delay of the reflection will vary; the amplitude of the reflection will vary; and a scattered coda will be generated. It is useful to model these effects in order to understand their importance and possibly to correct for them during processing. Three different numerical methods to model the sea surface interaction in a marine seismic reflection experiment (the frequencies considered are in the band 10-100Hz) are compared: the finite-difference method, the spectral element method and the Kirchhoff method. A plane wave is incident at angles of 0 deg and 30 deg with respect to the vertical on a rough Pierson-Moskowitz surface with 2or 4m significant waveheight and the response is synthesized at 6m, 10m and 50m below the average height of the sea-surface. All three methods display an excellent agreement for the main reflected arrival. The finite-difference and spectral element methods also agree very well all through the scattered coda. The Kirchhoff method shows some discrepancy, particularly in terms of amplitudes.