The widely used analytical models of flow and transport through aquifers of spatially, variable conductivity are based oil a first-order approximation in sigma(2)(Y), and the variance of the log conductivity Y = In K. With Y modeled as random, a few flow (effective conductivity, head and velocity covariances) and transport (macrodispersivity) statistical moments depend oil the mean and two-point covariance of Y only. Thus, the structure is essentially characterized by three parameters: K-G (the geometric mean), sigma(2)(Y) and I-Y (integral scale, with two different values for horizontal and vertical directions for anisotropic aquifers). Numerical modeling requires knowledge of the entire statistical structure of Y and the traditional approach is to regard it as multi-Gaussian, i.e. the vector of Y values at different points is a, multivariate normal vector. Comparison with first-order results showed that the latter are quite robust and apply to sigma(2)(Y) as large as unity. Recently, there is considerable interest in behavior of highly heterogeneous formations, Which have Y values larger than unity. A few numerical studies were conducted primarily for 2D (two-dimensional) configurations. The multi-Gaussian model was questioned and alternative ones were considered. A brief review is presented, with emphasis on the multi-indicator model developed recently by the author in collaboration with A. Flori and I. Jankovic. The main aim of the presentation is to discuss a few conceptual problems arising in modeling of flow and transport in highly heterogeneous aquifers such as: the nonuniqueness of structure models in view of scarcity of field data; robustness of first-order approximation and their limitations; 2D versus 3D models; and macrodispersivity, its meaningfulness, anomalous and non-Gaussian behavior.