In this study, we make use of a nonstationary stochastic theory in studying solute flux through spatially nonstationary flows in porous media. The nonstationarity of flow stems from various sources, such as multi-scale, nonstationary medium features and complex hydraulic boundary conditions. These flow nonstationarities are beyond the applicable range of the 'classical' stochastic theory for stationary flow fields, but widely exist in natural media. In this study, the stochastic frames for flow and transport are developed through an analytical analysis while the solutions are obtained with a numerical method. This approach combines the stochastic concept with the flexibility of the numerical method in handling medium nonstationarity and boundary/initial conditions. It provides a practical way for applying stochastic theory to solute transport in complex groundwater environments. This approach is demonstrated through some synthetic cases of solute transport in multi-scale media as well as some hypothetical scenarios of solute transport in the groundwater below the Yucca Mountain project area. It is shown that the spatial variations of mean log-conductivity and correlation function significantly affect the mean and variance of solute flux. Even for a stationary medium, complex hydraulic boundary conditions may result in a nonstationary flow field. Flow nonstationarity and/or nonuniform distribution of initial plume (geometry and/or density) may lead to nonGaussian behaviors (with multiple peaks) for mean and variance of the solute flux. The calculated standard deviation of solute flux is generally larger than its mean value, which implies that real solute fluxes may significantly deviate from the mean predictions.