<p>This paper investigates a time-fractional Stokes system in a bounded convex domain, presenting an analysis of both semidiscrete and fully discrete schemes. The spatial discretization is carried out using a finite element Galerkin method while keeping the time variable continuous, and an operator-theoretic approach is employed for the error analysis. Optimal error estimates with respect to data regularity and approximation properties are derived for the velocity and pressure terms under both smooth and nonsmooth initial data. Particular care is taken to avoid imposing nonlocal compatibility conditions on the data, which are often difficult to verify in practice. In addition, two fully discrete schemes based on the backward Euler and second-order backward difference methods are analyzed using a convolution quadrature technique, and corresponding error estimates are established. The analysis is further extended to problems posed on nonconvex domains, to distributed-order time-fractional models, and nonconforming lowest-order Crouzeix–Raviart elements. Finally, a unified framework is developed to extend the results to linearized generalized viscoelastic two-grade fluid and viscoelastic Oldroyd-B models. Numerical results are also presented to support the theoretical findings.</p>