The embedding optimization modeling approach is adapted to aid sustainable groundwater quantity and quality management of complex nonlinear multilayer aquifers. Implicit block-centered finite difference approximations of the quasi three-dimensional unsteady flow equation and Galerkin finite element approximations of the two-dimensional advection-dispersion transport equation are embedded directly as constraints in the model. Also used are nonlinear constraints describing river-aquifer interflow, evapotranspiration, and vertical flow reduction due to unconfinement. These circumvent use of large numbers of integer variables. The use of both linear and nonlinear formulations in a cyclical manner reduces execution time and improves confidence in solution optimality. The methodology is demonstrated for Salt Lake valley where groundwater quantity and quality management are needed, the proportion of pumping cells and cells needing head constraint is large, and many flows are described by discrete nonlinear or piecewise linear functions.