Add to ORKGWe present the extension of the complete flux scheme to advection-diffusion-reaction systems. For stationary problems, the flux approximation is derived from a local system boundary value problem for the entire system, including the source term vector. Therefore, the numerical flux vector consists of a homogeneous and an inhomogeneous component, corresponding to the advection-diffusion operator and the source term, respectively. For time-dependent systems, the numerical flux is determined from a quasi-stationary boundary value problem containing the time-derivative in the source term. Consequently, the complete flux scheme results in an implicit semidiscretization. The complete flux scheme is validated for several test problems