Flow of uniform mean velocity U takes place in a formation of spatially variable, random conductivity K(x). Advective transport of a plume of an inert solute is investigated by the Lagrangean approach. The aim of the study is to determine the spatial moments of the plume, i.e., of fluid particle trajectories, for highly heterogeneous aquifers, for which σY > 1, where Y = ln K. A multi-indicator model of the permeability structure, which is different from the common multi-Gaussian one, is proposed: the formation is modeled as a collection of N blocks of different K(j). The structure is defined by the distribution of K(j), the blocks' shape, and the coordinates of their centroids. The following simplifications are adopted: the blocks are inclusions of a regular shape (circles, spheres for isotropic media investigated here) defined by the radius A, and the inclusions are not overlapping, and their centroids are distributed uniformly and independently in space. At the continuous limit the model is characterized by the joint pdf f(Y, A). The model is shown to be quite general and to comprise binary, bimodal, indicator variograms and unimodal distributions of Y as particular cases. The study is focused on the latter case, with Y normal N <⟨Y⟩, σY2> and stationary covariance of given integral scale IY; these are the parameters commonly estimated for sedimentary formations. This leaves still freedom in selecting the pdf f(A). The simple model selected for semianalytical and numerical analysis is that of inclusions of radius R and volume fraction n, submerged in a matrix of effective conductivity Kef. The latter represents the effect of inclusions of much smaller radius, which appear as a nugget in the log conductivity two-point covariance. An approximate analytical solution of the flow is obtained by using a self-consistent approximation, while a fully numerical one is derived in part 3 <Jankovic et al., 2003a>. Transport is solved by particle tracking, and the time-dependent spatial moments (trajectories variance, skewness, kurtosis) are presented in part 2 <Fiori et al., 2003>. In the self-consistent approximation the asymptotic longitudinal macrodispersivity αL, which is a function of Y, shows strong nonlinear effects: inclusions of large positive Y lead to a finite αL, whereas αL grows unbounded for those of negative Y. This effect is not captured by the common first-order approximation in σY, which is symmetrical and overestimates αL for Y > 0 and underestimates it for Y < 0. As a result, the second spatial moment is predicted accurately by the first-order approximation, by cancellation of errors, provided that f(Y) is symmetrical. However, the transient regime and higher-order moments are not captured by the first-order approximation.