We investigate a generalization of cubic splines to Riemannian manifolds. Spline curves are defined as minimizers of the spline energy-a combination of the Riemannian path energy and the time integral of the squared covariant derivative of the path velocity-under suitable interpolation conditions. A variational time discretization for the spline energy leads to a constrained optimization problem over discrete paths on the manifold. Existence of continuous and discrete spline curves is established using the direct method in the calculus of variations. Furthermore, the convergence of discrete spline paths to a continuous spline curve follows from the Γ-convergence of the discrete to the continuous spline energy. Finally, selected example settings are discussed, including splines on embedded finite-dimensional manifolds, on a high-dimensional manifold of discrete shells with applications in surface processing, and on the infinite-dimensional shape manifold of viscous rods.