Water resources systems management often requires complex mathematical models whose use may be computationally infeasible for many advanced analyses, e.g., optimization, data assimilation, model uncertainty, etc. The computational demand of these analyses can be reduced by approximating the model with a simpler reduced model. Proper Orthogonal Decomposition (POD) is an efficient model reduction technique based on the projection of the original model onto a subspace generated by full-model snapshots. In order to implement this method, an appropriate number of snapshots of the full model must be taken at the appropriate times such that the resulting reduced model is as accurate as possible. Since confined aquifers reach steady state in an exponential manner, we present a simple exponential function that can be used to select snapshots for these types of models. This selection method is then employed to determine the optimal snapshot set for a unit, dimensionless model. The optimal snapshot set is found by maximizing the minimum eigenvalue of the snapshot covariance matrix, a criterion similar to those used in experimental design. The resulting snapshot set can then be translated to any complex, real world model based on a simple, approximate relationship between dimensionless and real-world times. This translation is illustrated using a basin scale model of Central Veneto, Italy. We show that this method is very easy to implement and produces an accurate reduced model that, in the case of Central Veneto, Italy, runs approximately 1,000 times faster than the full model.