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.