Three dimensional aquifer deformation and flow modeling

Abstract

Horizontal strain and three-dimensional aquifer deformation resulting from groundwater pumping have been shown to represent important parameters when simulating land subsidence at the basin scale. Horizontal strain is typically ignored in available programs for simulating groundwater flow and subsidence. When the strain field within an aquifer due to an applied stress is not simplified to the one-dimensional Terzaghi assumption, the resulting governing equations yield four dependent variables—three for the displacement field of solids and one for hydraulic head. Typical solution schemes for three-dimensional field-scale problems that include the displacement field and hydraulic head are not practical or possible due to the large amount of storage required by such systems of equations. An efficient one-dimensional row-indexed storage scheme that takes advantage of the sparse matrices produced by this system of equations stores only about two times the number of non-zero matrix elements. The storage scheme is used in conjunction with a preconditioned conjugate-gradient solver that involves the implementation of several row-indexed storage arrays, greatly reducing the required storage and yielding a rapid solution for even large field problems with large grids networks. The modified incomplete-Cholesky factorization scheme is used to produce the preconditioning matrix. The granular displacement model, or GDM, is written as a module for the U.S. Geological Survey's MODFLOW finite-difference flow model to provide fully three-dimensional deformation capability for subsidence prediction. The code is based on Biot's three-dimensional consolidation theory. The displacement field of solids and hydraulic head dependent variables are weakly coupled allowing for the head to be calculated with the original MODFLOW code. The head solution is then used in the solution for the displacement field. The resulting displacements result in a change in storage, which affects the heads. The heads are then recalculated during the same time step until a user specified convergence is reached. Once both the calculated head and displacement field have converged for a given time step, the model moves on to the next time step until the total simulation time is reached. The GDM program has an advantage over one-dimensional subsidence models because the one-dimensional models tend to over-predict subsidence, particularly near pumping wells. In addition one-dimensional models predict more subsidence in far-field areas from the well. In these regions the predictions using the GDM program indicate that storage released from horizontal deformation is large.