A three-dimensional numerical method of moments has been developed for solute flux through nonstationary flows in porous media. The solute flux is described as a space-time process where time refers to the solute flux breakthrough and space refers to the transverse displacement distribution at a control plane. Flow nonstationarity may stem from various sources, such as the medium's conductivity nonstationarity and complex hydraulic boundary conditions. The first two statistics of solute flux are derived using a Lagrangian framework and are expressed in terms of the probability density functions (PDFs). These PDFs are given in terms of one-and two-parcel moments of travel time and transverse locations, and these moments are related to the Eulerian velocity moments. The moment equations obtained analytically for flow and transport are so complex that numerical techniques are used to obtain solutions. In this study, we investigate the influence of various factors, such as the grid resolution relative to correlation length and the number of solute parcels comprising a source, on the accuracy of the calculation results. It has been found that for the computation of means and variances using the developed moment equations, hydraulic head requires at least one numerical grid element per correlation length scale. At least two grid elements are required for velocity, and 1-2 grid elements for the solute flux variance. Five parcels are required per correlation length scale to approximate the initial solute source distribution. The effects of boundary and hydraulic nductivity nonstationarity on flow and transport are also considered. Flow nonstationarity caused by either hydraulic boundary condition or conductivity nonstationarity significantly influences the transport process. The calculation results of numerical method of moments are compared with Monte Carlo simulations. The comparison indicates that the two methods are consistent with each other for head variance, velocity covariance in longitudinal direction, and mean and variance of total solute flux, but numerical method of moment underestimates the velocity variance in transverse direction. The method is applied to an environmental project for predicting the solute flux in the saturated zone below the Yucca Mountain project area, demonstrating the applicability of the method to complex subsurface environments.