In general, stochastic partial differential equations and/or boundary conditions describing flow in phreatic aquifers are nonlinear. This paper uses second-order perturbation to approximate the response of a phreatic aquifer system subject to multiple random inputs. The approach is equivalent to assuming quadratic nonlinearities in the original nonlinear phreatic aquifer system. Using frequency domain representations, the resulting model is decomposed into linear and nonlinear (quadratic) frequency response functions (FRF). For Gaussian input, the response of a nonlinear phreatic aquifer system is non-Gaussian and is typically skewed. Thus the variance spectrum is insufficient to describe the response process. The bispectrum, defined as the Fourier transform of the third-order cumulant sequence, represents the distribution of skewness with frequency and further serves as a measure of the nonlinearity of the system. Given input processes of aquifer recharge and stream stage fluctuations, the linear and quadratic FRF are used to construct the spectrum and bispectrum of the response process. The role of nonlinearity in a stream-aquifer system is examined. ¿ American Geophysical Union 1994