Oscillatory flow produced near a seabed by linear surface gravity waves is analyzed. Using the boundary layer approximation, the problem is reduced to that of a boundary layer near a plate harmonically oscillating in a fluid at rest. Saffman's turbulence model is employed. Results on the ensemble averaged quantities are presented for Reynolds numbers ranging from laminar to fully developed turbulent regime, and for small, intermediate and large values of the roughness parameter. Based on the above model, a numerical fully nonlinear stability analysis of Stokes' solution is performed, and the unsteady transition from laminar to turbulent regime is investigated. The reliability of the analytical results is supported by a comparison with experimental data collected by others.