In this paper, we present a solute flux approach for analyzing solute transport statistics in statistically nonstationary, unsaturated flow. Flow nonstationarity in the vadose zone may arise from a number of factors. It is useful to develop a systematic approach that incorporates these factors into an uncertainty analysis. We first derive the general forms for solute flux moments. The solute flux moments are associated with one- and two-particle joint probability distribution functions (JPDF). We illustrate our results for certain forms of one-particle and two-particle JPDFs, in which the particle travel time is assumed to be lognormally distributed and the particle transverse displacement normally distributed. In the numerical examples, the Eulerian velocity moments is obtained by solving the head moment equations numerically using a finite-difference method. Our results show that flow nonstationarity has a significant impact on the statistics of solute fluxes and solute breakthrough curves.