This paper proposes a hybrid mixed-integer linear programming-linear programming-linear programming (MILP-LP-LP) methodology for simultaneous optimal design and operation of unconfined groundwater utilization systems. The proposed model is an extension to the earlier LP-LP model introduced by the authors for the optimal design and operation of confined groundwater utilization systems. The proposed model can be used to minimize the total cost, including the well drilling, pump installation cost and energy cost, of utilizing a two-dimensional unconfined aquifer under both steady-state and transient flow conditions. The solution of the problem is defined by the well numbers and location, well drilling depth and the corresponding pumping rates, satisfying a downstream demand, lower/upper bound on the pumping rates, and lower/upper bound on the water level drawdown in the wells. A discretized version of the differential equation governing the aquifer is first embedded into the model formulation as additional constraints. The resulting mixed-integer non-linear programming problem is then decomposed into three sub-problems with different sets of decision variables, namely, hydraulic head, well numbers and locations and their pumping rates, and well depths. Starting with a set of random values for all decision variables, the three sub-problems are solved assuming fixed values for the other variables. This process is iterated until convergence is achieved.