We perform high-resolution Monte Carlo simulations to study steady-state flow in three-dimensional leaky heterogeneous aquitards. The vertical extend of the aquitards is finite, and constant-mean head boundaries apply in the planar direction. The loghydraulic conductivity, Y(x), is considered to be Gaussian, statistically homogeneous, and described by an exponential auto-covariance function. We have performed numerical simulations to investigate the effect of (1) the dimensionless ratio L/λ, L being the aquitards' vertical extend and λ being the integral scale of Y(x); and (2) the variance of Y(x). The turning bands and sequential Gaussian simulation methods were used to generate the heterogeneous Y(x) fields, and a block-centered, seven-point finite difference scheme was utilized to solve the steady-state, three-dimensional, saturated flow equation. We have found that for small ratios of L/λ, the arithmetic mean of hydraulic conductivities serves as a good approximation to equivalent homogeneous aquitards, whereas the behavior at large ratios (L/λ > 4) approaches that of equivalent infinite aquitards. Our results, for various finite ratios L/λ and for variances of Y(x) to 7, support the analytical expressions developed in 1996 by Paleologos, Neuman, and Tartakovsky for bounded, strongly heterogeneous media. Our numerical results indicate that for statistically isotropic media and for infinite and/or finite flow domains, the exponential generalization of theoretical expressions developed through first-order analyses provides an accurate prediction of ensemble mean hydraulic behavior.