To generate multidimensional Gaussian random fields over a regular sampling grid, hydrogeologists can call upon essentially two approaches. The first approach covers methods that are exact but computationally expensive, e.g., matrix factorization. The second covers methods that are approximate but that have only modest computational requirements, e.g., the spectral and turning bands methods. In this paper, we present a new approach that is both exact and computationally very efficient. The approach is based on embedding the random field correlation matrix R in a matrix S that has a circulant/block circulant structure. We then construct products of the square root S1/2 with white noise random vectors. Appropriate subvectors of this product have correlation matrix R, and so are realizations of the desired random field. The only conditions that must be satisfied for the method to be valid are that (1) the mesh of the sampling grid be rectangular, (2) the correlation function be invariant under translation, and (3) the embedding matrix S be nonnegative definite. These conditions are mild and turn out to be satisfied in most practical hydrogeological problems. Implementation of the method requires only knowledge of the desired correlation function. Furthermore, if the sampling grid is a d-dimensional rectangular mesh containing n points in total and the correlation between points on opposite sides of the rectangular is vanishingly small, the computational requirements are only those of a fast Fourier transform (FFT) of a vector of dimension 2dn per realization. Thus the cost of our approach is comparable with that of a spectral method also implemented using the FFT. In summary, the method is simple to understand, easy to implement, and is fast. ¿ American Geophysical Union 1993