AbstractNumerical models are often used for simulating groundwater flow. Written in state‐space form, the dimension of these models is of the order of the number of grid cells used and can be very high (more than a million). As a result, these models are computationally very demanding, especially if many different scenarios have to be simulated. In this paper we introduce a model reduction approach to develop an approximate model with a significantly reduced dimension. The reduction method is based upon several simulations of the large‐scale numerical model. By computing the covariance matrix of the model results, we obtain insight into the variability of the model behaviour. Moreover, by selecting the leading eigenvectors of this covariance matrix, we obtain the spatial patterns that represent the directions in state space where the model variability is dominant. These patterns are also called empirical orthogonal functions. We can project the original numerical model onto those dominant spatial patterns. The result is a low‐dimensional model that is still able to reproduce the dominant model behaviour. Copyright © 2004 John Wiley & Sons, Ltd.