This study presents a hybrid spectral method (HSM) to estimate flow uncertainty in large-scale highly nonstationary groundwater systems. Taking advantages of spectral theories in solving unmodeled small-scale variability in hydraulic conductivity, the proposed HSM integrates analytical and numerical spectral solutions in the calculation procedures to estimate flow uncertainty. More specifically, the HSM involves two major computational steps after the mean flow equation is solved. The first step is to apply an analytical-based approximate spectral method (ASM) to predict nonstationary flow variances for entire modeling area. The perturbation-based numerical method, nonstationary spectral method (NSM), is then employed in the second step to correct the regional solution in local areas, where the variance dynamics is considered to be highly nonstationary (e.g., around inner boundaries or strong sources/sinks). The boundary conditions for the localized numerical solutions are based on the ASM closed form solutions at boundary nodes. Since the regional closed form solution is instantaneous and the more expensive perturbation-based numerical analysis is only applied locally around the strong stresses, the proposed HSM can be very efficient, making it possible to model strongly nonstationary variance dynamics with complex flow situations in large-scale groundwater systems. In this study the analytical-based ASM solutions was first assessed to quantify the solution accuracy under transient and inner boundary flow conditions. This study then illustrated the HSM accuracy and effectiveness with two synthetic examples. The HSM solutions were systematically compared with the corresponding numerical solutions of NSM and Monte Carlo simulation (MCS), and the analytical-based solutions of ASM. The simulation results have revealed that the HSM is computationally efficient and can provide accurate variance estimations for highly nonstationary large-scale groundwater flow problems. (C) 2009 Elsevier B.V. All rights reserved.