Abstract The mathematical model and its numerical solution presented here were constructed for a hydraulic system composed of two deep aquifers that are partially separated by an aquiclude in the Negev Desert, Israel. Groundwater in these aquifers is confined, except for small areas, and it is not homogeneous because of large variations in salinity and temperature. The model is designed to predict aquifer response to pumpage in terms of water pressure and density. Increasing the rate of pumpage is likely to enlarge the unconfined area. The conceptual model has been transformed to a three-dimensional mathematical model. Two factors made it possible to reduce the mathematical model to two dimensions in each aquifer by averaging equations along the vertical: little variation in state variables along the vertical plane and the presence of a hydrostatic distribution of pressure. We used differential mass and volume balance equations for which the averages of water density and pressure along the vertical at any given point were the dependent variables. The piezometric head at any geographical point, however, could be used as a means for distinguishing between confined and water-table conditions in the hydraulic system at the initial state and at various pumpage states. The transformation of these equations into finite differences results in a large system of nonlinear equations. The sparse structure of this system suggests the use of a succession over relaxation (S.O.R.) scheme. However, the presence of highly nonlinear terms in the difference equations introduces considerable difficulty in applying this method. We therefore designed a special linearization technique to overcome this difficulty. Other features that complicate the solution process are the dynamic confined-unconfined nature of the aquifers and the irregular shape of the aquiclude. We tested and verified the validity of the model by a series of simple test cases and then used the model to study aquifer response to various pumpage alternatives.