Based on a three-dimensional heterogeneous aquifer model exhibiting non-stationary, statistically anisotropic correlation, three hydrostratigraphic models (HSMs) are created within a sedimentary hierarchy. A geostatistical analysis of natural log conductivity (lnK) is conducted for the units of the HSMs. Hydraulic conductivity is then upscaled using numerical and analytical methods. Increasing lnK variances are evaluated. Results suggest that for the aquifer model tested: (1) the numerical method is capable of upscaling irregular domains with reasonable accuracy for a lnK variance up to 7.0. (2) Accuracy of the upscaled equivalent conductivities (K*) and associated performance of the HSMs are sensitive to homogenization level, heterogeneity variance, and boundary condition. Variance is found to be the most significant factor impacting the accuracy of the HSMs. (3) Diagonal tensor appears a good approximation for the full-tensor K*. (4) For the HSM units, when the variance is low (less than 1.0), all analytical methods are nearly equally accurate; however, when variance becomes higher, analytical methods generally are less accurate.