Nonstationarity of flow fields due to pumping wells and its impact on advective transport is of particular interest in well capture zone design and wellhead protection. However, techniques based on Monte Carlo methods to characterize the associated capture zone uncertainty are time consuming and cumbersome. This paper introduces an alternative approach. The mean and covariance of system state variables (i.e., head, pore water velocity, and particle trajectory) are approximated using a first-order Taylor's series with sensitivity coefficients estimated from the adjoint operator for a system of discrete equations. The approach allows nonstationarity due to several sources (e.g., transmissivity, pumping, boundary conditions) to be treated. By employing numerical solution methods, it is able to handle irregular geometry, varying boundary conditions, complicated sink/source terms, and different covariance functions, all of which are important factors for real-world applications. A comparison of results for the Taylor's series approximation with those from Monte Carlo analysis showed, in general, good agreement for most of the tested particles. Particle trajectory variance calculated using Taylor's series approximation is then used to predict well capture zone probabilities under the assumption of normality of the mass transport's state variables. Verification of this assumption showed that not all particle trajectories (depending on their starting location) are normally or log-normally distributed. However, the risk of using the first-order method to delineate the confidence interval of a well capture zone is minimal since it marginally overestimates the 2.5% probability contour. Furthermore, this should be balanced against its greater computation efficiency over the Monte Carlo approach.