A method for processing energy at a free-surface reflection relating to an air-water interface, the method comprising: receiving input seismic data recorded with seismic sensors; receiving wave-height data that describes an actual shape of a top surface of a body of water; separating up-going energy at a receiver and down-going energy following a reflection at the sea-surface, using the input seismic data and a linear operator modified to take into account the wave-height data, wherein said modified linear operator is a modified representation of a linear least squares Radon algorithm using a tau-p transform solving a phase shift equation d = L p, where d represents the input data for a 2D shot-point gather of n traces for one frequency slide, p is the tau-p domain model containing m slowness traces and L is the reverse slant stack operator, wherein said modified operator is given by L n,m = e - 2 πifτ pr + Re - 2 πifτ gh , where R is the reflection coefficient at the sea surface and τ pr and τ gh are the modified time delays of the primary and ghost components, respectively, where the modified time delay for the primary component is given by τ pr = ( h n + Δ h nm )s m - Δτ nm and the modified time delay for the ghost is given by z gh = ( h n - Δ h nm - Δh wave ) s m + Δ τ nm + Δ τ wave , where h n is an offset between a source and a seismic sensor, Δh nm is an offset difference, s m is the slowness of the p-trace, Δτ nm is the travel time of the ghost from a horizontal water surface to the seismic sensor for trace n and slowness m, Δh wave is the distance traveled by the down-going wave field along the flat sea surface due to the wave's height and Δτ wave is the time delay introduced by the water's height; and generating an image of the subsurface based on the up-going energy or the down-going energy or a combination of the input seismic data and one of the up-going or down-going energy. receiving input seismic data recorded with seismic sensors; receiving wave-height data that describes an actual shape of a top surface of a body of water; separating up-going energy at a receiver and down-going energy following a reflection at the sea-surface, using the input seismic data and a linear operator modified to take into account the wave-height data, wherein said modified linear operator is a modified representation of a linear least squares Radon algorithm using a tau-p transform solving a phase shift equation d = L p, where d represents the input data for a 2D shot-point gather of n traces for one frequency slide, p is the tau-p domain model containing m slowness traces and L is the reverse slant stack operator, wherein said modified operator is given by L n,m = e - 2 πifτ pr + Re - 2 πifτ gh , where R is the reflection coefficient at the sea surface and τ pr and τ gh are the modified time delays of the primary and ghost components, respectively, where the modified time delay for the primary component is given by τ pr = ( h n + Δ h nm )s m - Δτ nm and the modified time delay for the ghost is given by z gh = ( h n - Δ h nm - Δh wave ) s m + Δ τ nm + Δ τ wave , where h n is an offset between a source and a seismic sensor, Δh nm is an offset difference, s m is the slowness of the p-trace, Δτ nm is the travel time of the ghost from a horizontal water surface to the seismic sensor for trace n and slowness m, Δh wave is the distance traveled by the down-going wave field along the flat sea surface due to the wave's height and Δτ wave is the time delay introduced by the water's height; and generating an image of the subsurface based on the up-going energy or the down-going energy or a combination of the input seismic data and one of the up-going or down-going energy. <EOS>