We develop a three-dimensional numerical model to simulate bed form dynamics in a turbulent boundary layer. In the numerical model, hydrodynamics is solved in a moving generalized boundary-fitted curvilinear coordinate system, such that the domain boundary exactly follows complex time-dependent bed form geometry. The resolved turbulent features are computed via large-eddy simulation, while the subgrid scale turbulent motions are modeled with a dynamic mixed model. A second-order accurate arbitrary Lagrangian-Eulerian method is used to guarantee conservation of sediment mass, while the grid moves arbitrarily due to the motion of the bed. Transport of suspended load is modeled using the Eulerian approach with a pickup function as the bottom boundary condition for sediment entrainment at the bed. Transport of bed load and suspended load are combined in a mass balance equation for the bed, which evolves due to the spatiotemporally varying bed stress induced by the turbulent flow field above the bed and gravity (gravity-induced avalanche flow). Motion of the bed in turn affects the flow field in a coupled hydrodynamic moving bed simulation, in which bed features evolve due to resolved details of the turbulent flow. We compare different bed elevation models and demonstrate the capability of the present model through simulation of sand ripple formation and evolution induced by turbulence in an oscillatory flow. A resolution study demonstrates the need for fine grid resolution to resolve a bulk of the near-wall turbulence, which is essential for bed form initiation.