In the research to be presented, we express concentration or flux of solutes as a distribution over their velocity. We then derive an integrodifferential equation that governs the plume evolution and multi-component reactions solute distribution over velocity at given times and locations for a particle ensemble, based on a presumed velocity correlation structure and an ergodic cross-sectional velocity distribution. This way, the spatial evolution of breakthrough curves away from the source and the effective reaction rates are predicted based on cross-sectional velocity distribution and the connectivity, which is expressed by the velocity transition probability density. The transition probability can be specified via a copula function that can help construct a joint distribution with a given correlation and given marginal velocities or or in the special case of Gaussian copula, via an Ornstein-Uhlenbeck process. Using this approach, we analyze the breakthrough curves depending on the velocity distribution and correlation properties.