A mathematical model and a corresponding coupled numerical model that incorporate the concepts of three dimensional poroelasticity based on Biot's consolidation theory are developed for simulating the displacement field of solids within unconsolidated aquifers in response to induced changes in water pressure. The granular displacement model (or GDM) is incorporated into MODFLOW, a U.S. Geological Survey's finite-difference program, and uses an efficient row-indexed storage mode and biconjugate gradient solver, so that it is tractable at the field scale. The application to three dimensions is an improvement over two-dimensional axisymmetric models incorporating Biot's equations and over the commonly used one-dimensional subsidence models based on Terzaghi's method of effective stress. Most subsidence due to ground-water withdrawal occurs in the inelastic range of specific storage within clay interbeds or confining units. However, most horizontal deformation occurs within coarser aquifer units. The GDM program uses both elastic and inelastic storage and Poisson's ratio as key parameters. Conversion from elastic to inelastic specific storage occurs when the previous maximum volume strain for a particular cell is exceeded. Model outputs are compared with a currently available one-dimensional subsidence model (IBS1) developed for MODFLOW and an axisymmetric finite-element Biot program. The results indicate that under traction free conditions subsidence is a three-dimensional problem and one-dimensional subsidence models tend to focus excessive amounts of vertical deformation near the pumped well. The magnitude of vertical deformation in one-dimensional subsidence models is exacerbated as the grid size becomes smaller in the vicinity of the pumping well. This is due to increased calculated drawdown in the vicinity of the well for more finely-spaced grids.