In diffusion transport, the flux at a point is typically modeled in terms of the local gradient of a potential. When heterogeneities are present, this local model can break down and it may be more appropriate to model the diffusion flux as a weighted sum of gradients present throughout the domain. Here a discrete nonlocal flux modelconsistent with control-volume implementations-is developed. This scheme is referred to as the control-volume weighted flux scheme (CVWFS). The key component is the modeling of the diffusion flux at a given control-volume face in terms of a weighted sum of gradients at that face and at faces up- and downstream. Criteria for choosing the weights are proposed. This results in numerical solution schemes in which the coefficient matrix is diagonally dominant, has positive off-diagonal elements, and zero row sums. For a particular power-law weighting scheme it is shown how the CVWFS is related to the definition of the Caputo fractional derivative and the one-shift Grunwald approximation of the Riemann-Liouville fractional derivative. On developing transients and boundary condition treatments, the accuracy and suitability of the CVWFS scheme is demonstrated by solving a number of problems governed by Caputo fractional diffusion equations.