A new MODFLOW package (Nonlinear Flow Process; NLFP) simulating nonlinear flow following the Forchheimer equation was developed and implemented in MODLFOW-2005. The method is based on an iterative modification of the conductance calculated and used by MODFLOW to obtain an effective Forchheimer conductance. The package is compatible with the different layer types, boundary conditions, and solvers as well as the wetting capability of MODFLOW. The correct implementation is demonstrated using four different benchmark scenarios for which analytical solutions are available. A scenario considering transient flow in a more realistic setting and a larger model domain with a higher number of cells demonstrates that NLFP performs well under more complex conditions, although it converges moderately slower than the standard MODFLOW depending on the nonlinearity of flow. Thus, this new tool opens a field of opportunities to groundwater flow simulation with MODFLOW, especially for core sample simulation or vuggy karstified aquifers as well as for nonlinear flow in the vicinity of pumping wells.

The theory of gravity-driven regional groundwater flow was first proposed in 1962/3 based on the Laplace equation. Hydraulic-head patterns were calculated for a two dimensional trapezoidal and homogeneous flow domain with flow lines drawn by hand. The flow region was intended to represent one flank of a stream basin with a periodically undulating water table. At the dawn of numerical modeling the results generated international interest. Numerical models began to be produced with progressively increasing complexity of basin geometry, types and distributions of permeability and time dependent flow. One of the most important results of the first analyses was the birth of the flow-system concept. In a flow system groundwater moves from relatively highly elevated recharge areas, through medium high mid-line regions to relatively low lying discharge areas where it may resurface. Because flow systems are associated with topographic elements of different scale, they are self-organized in hierarchically nested geometric patterns.
The understanding of the systematized structure of basinal groundwater flow soon resulted in the recognition that flow systems act like subsurface conveyor belts. They mobilize and remove matter and heat from the recharge area, pick up more or/and emplace some of it en route, and deposit them in the discharge region. In short: flowing groundwater is a general geologic agent. The original „Theory of regional groundwater flow” became thus expanded into a bimodal umbrella theory with two component theories: i) „The hydraulics of basin-scale groundwater flow” and ii) „The geologic agency of regional groundwater flow”. More than half a century after its conception the theory is extensively analyzed and continues to be applied to a growing number of groundwater related disciplines