Combined simulation and optimization models, which are helpful for long-term groundwater planning of complex nonlinear aquifer systems, are developed using alternative modelling approaches. The models incorporate a representation of steady-state, quasi-three-dimensional head response to pumping within an optimization . An embedding model which describes exactly the nonlinear flow of an unconfined aquifer is presented. In contrast with the embedding models presented in the Utah State University Ground Water Model, it directly achieves the optimal solution without a "cycling." To address the nonlinearity of the flow system, response matrix models couple superposition with the cycling procedure. Their linear influence coefficients are generated using a modified McDonald and Harbaugh model. First, these models are tested for a hypothetical, 625 cell, nonlinear aquifer system and compared in terms of computational accuracy and efficiency. All of the models achieve the same optimal solution. The fully nonlinear embedding model attains the same optimal solution regardless of how far the initial guess is from that solution. Thus, global optimality is probably obtained. A predictive program for comparing a priori the embedding and response matrix models in terms of computational size is also developed. This computes the required memory for running each model, an important factor in computational efficiency. It is based on the number of nonzero elements in the matrix of the optimization scheme. The model most appropriate for a given aquifer and desired management scenarios is dependent upon required simulation accuracy, flow conditions (steady or unsteady) , spatial scale, model computational resources requirement, and the computational capacity of available hardware and software. The linear embedding model coupled with a cycling procedure, as incorporated within a modified version of the USUGWM, is most appropriate for the subject reconnaissance level study of the East Shore Area. Here, the demand for sufficient water of adequate quality is increasing. The underlying aquifer is three-layered, unconfined/confined and is discretized into 4,880 finite-difference cells. To overcome the difficulties of solving many nonsmooth functions describing evapotranspiration, discharge from flowing wells, and drain discharge, a former cycling procedure is improved by optimizing the purely linearized models repeatedly. Using the modified version of the USUGWM, optimal sustained-yield pumping strategies are computed for alternative future scenarios in the East Shore Area.