Abstract In this paper, we develop a numerical scheme for the space-time fractional parabolic equation, i.e. an equation involving a fractional time derivative and a fractional spatial operator. Both the initial value problem and the non-homogeneous forcing problem (with zero initial data) are considered. The solution operator E ⁢ ( t ) {E(t)} for the initial value problem can be written as a Dunford–Taylor integral involving the Mittag-Leffler function e α , 1 {e_{\alpha,1}} and the resolvent of the underlying (non-fractional) spatial operator over an appropriate integration path in the complex plane. Here α denotes the order of the fractional time derivative. The solution for the non-homogeneous problem can be written as a convolution involving an operator W ⁢ ( t ) {W(t)} and the forcing function F ⁢ ( t ) {F(t)} . We develop and analyze semi-discrete methods based on finite element approximation to the underlying (non-fractional) spatial operator in terms of analogous Dunford–Taylor integrals applied to the discrete operator. The space error is of optimal order up to a logarithm of 1 h {\frac{1}{h}} . The fully discrete method for the initial value problem is developed from the semi-discrete approximation by applying a sinc quadrature technique to approximate the Dunford–Taylor integral of the discrete operator and is free of any time stepping. The sinc quadrature of step size k involves k - 2 {k^{-2}} nodes and results in an additional O ⁢ ( exp ⁡ ( - c k ) ) {O(\exp(-\frac{c}{k}))} error. To approximate the convolution appearing in the semi-discrete approximation to the non-homogeneous problem, we apply a pseudo-midpoint quadrature. This involves the average of W h ⁢ ( s ) {W_{h}(s)} , (the semi-discrete approximation to W ⁢ ( s ) {W(s)} ) over the quadrature interval. This average can also be written as a Dunford–Taylor integral. We first analyze the error between this quadrature and the semi-discrete approximation. To develop a fully discrete method, we then introduce sinc quadrature approximations to the Dunford–Taylor integrals for computing the averages. We show that for a refined grid in time with a mesh of O ⁢ ( 𝒩 ⁢ log ⁡ ( 𝒩 ) ) {O({\mathcal{N}}\log({\mathcal{N}}))} intervals, the error between the semi-discrete and fully discrete approximation is O ⁢ ( 𝒩 - 2 + log ⁡ ( 𝒩 ) ⁢ exp ⁡ ( - c k ) ) {O({\mathcal{N}}^{-2}+\log({\mathcal{N}})\exp(-\frac{c}{k}))} . We also report the results of numerical experiments that are in agreement with the theoretical error estimates.