We consider three-dimensional steady state flow towards a well that fully penetrates a randomly heterogeneous aquifer confined between horizontal no-flow boundaries, and bounded laterally by a cylindrical, deterministically prescribed constant head boundary. The well is represented by a line sink that produces water at a deterministically prescribed constant rate Q for unit aquifer thickness. The log hydraulic conductivity, Y = ln K, of the aquifer is multivariate Gaussian, statistically homogeneous with a Gaussian spatial autocorrelation function. We develop an analytical solution for the variance of hydraulic head as a function of dimensionless vertical and horizontal locations within the aquifer, variance σ_Y² of Y and dimensionless ratios between the principal spatial correlation scales. Our analysis is based on the non-local theory first proposed for steady state and transient flows in bounded, randomly heterogeneous media by Neuman and Orr (1993), Neuman et al. (1996), Guadagnini and Neuman (1999a,b) and Tartakovsky and Neuman (1998, 1999). In particular, we develop and solve analytically recursive closure approximations of the governing non-local moment equations to second order in σ_Y by means of an appropriate Green's function. We evaluate our analytical solutions by means of Gaussian quadratures for the special case of the isotropic Y field.