A comprehensive numerical study of oscillatory wave boundary layers on spatially varying bottom roughness is presented. The study utilizes a model solving incompressible Reynolds-averaged Navier-Stokes equations coupled with k-w turbulence closure, modified in a simple way to incorporate anisotropy in turbulent normal stresses. The model is first validated via comparison with existing oscillating tunnel measurements involving sudden bottom roughness transitions. It is then used to parametrically study oscillatory boundary layer flows, wherein the bed shear stress amplifications and period-averaged streaming characteristics induced by bottom roughness variations are systematically assessed. The effects of variable roughness ratio, gradual roughness transitions, as well as changing flow orientation in plan are all considered. As part of the latter, roughness-induced secondary flows are predicted to occur as the oscillatory flow becomes oriented parallel to a line of roughness transition. This phenomenon is proposed as a natural transverse grain sorting mechanism for coastal flows over graded sediments. Subsequent model testing demonstrates potential generation of secondary circulation cells having characteristic size the order of the wave boundary layer thickness. Analogy is made to similar features known to develop within steady flows, having characteristic size the order of the flow depth.