Direct numerical simulations of granular beds consisting of uniformly sized spherical particles being eroded by a shear flow of Newtonian liquid have been performed. The lattice-Boltzmann method has been used for resolving the flow of the interstitial liquid. Fluid and solid dynamics are fully coupled with the particles having finite size and undergoing hard-sphere collisions. Only laminar flow has been considered with particle-based Reynolds numbers in the range 0.02 to 0.6. The parameter range of the simulations covers the transition between static and mobilized beds. The transition occurs for 0.10<θ<0.150.10<θ<0.15 with θθ the Shields number. The transition is insensitive of the Reynolds number and the solid-over-liquid density ratio. Incipient bed motion has been interpreted in terms of the probability density functions of the hydrodynamic forces acting on the spheres in the upper layer of the bed.