Using the method of decomposition, new analytical solutions of the nonlinear Boussinesq flow equation and of the exact two-dimensional groundwater flow equation subject to a nonlinear free-surface boundary condition are presented and tested with respect to the linearized Boussinesq equation. It is found that for mild regional gradients and for the range of recharge values usually encountered in the field, the extensively used linearized equation with Dupuit assumptions is a reasonable approximation to the exact solution for the hydraulic heads and the regional flow velocities. Discrepancies occur in the presence of high regional hydraulic gradients, unusually high recharge rates, or regions of low hydraulic conductivity. With this result a new nonperturbation solution of the two-dimensional plan view groundwater flow equation in a heterogeneous (heterogeneity represented as a spatial random field) aquifer is presented, and expressions for the statistical properties of the longitudinal and transverse velocities are derived. The results suggest that scale dependency, or spatial variability, in the flow velocity arises naturally as a result of recharge from rainfall and aquifer heterogeneity and may help explain the scale dependency of aquifer dispersion parameters. ¿ American Geophysical Union 1995