The design of an optimal pumpage or recharge strategy for contaminated groundwater remediation is considered. The formulation of the control problem combines a two-dimensional Galerkin finite element simulation model of steady state groundwater flow and transient convective-dispersive transport with non-linearly constrained optimization algorithms. The ‘discrete gradients’ of the constraints are computed by the ‘discrete adjoint state method’ or by the solution of ‘discrete sensitivity equations’, according to the number of control parameters compared to the number of gradients to be computed. This technique is shown to be very efficient when used in conjunction with the ‘constrained variable metric method’ of Han, modified by Schittkowsky. Numerical experiments indicate that problems with up to fifty injection or extraction wells may be solved in a computer time equivalent to one or two hundred simulations of the state equations.