In addressing the multi-criteria Wind Farm Layout Optimization (WFLO) problem, the literature has been focused on the simple weighted sum approach using single-objective stochastic and evolutionary algorithms, in addition to Pareto formulations using evolutionary algorithms. There is no single solution to a multi-criteria problem with conflicting objectives; therefore, the Pareto approach is useful to provide the developer with a non-dominated set of solutions. However, the evolutionary optimization algorithms tend to be computationally prohibitive, especially when optimizing large-scale wind farms. Additionally, most of WFLO problems are highly constrained, where many unfeasible zones can exist inside the proposed wind farm boundaries, which in turn complicates the optimization process. To remedy these drawbacks, we propose a gradient-based approach to Pareto optimization of the multi-criteria WFLO problems considering land footprint, energy output, electrical infrastructure and environmental impact. Mathematical functions and their derivatives are developed to represent the four objectives, land based constraints, and their gradients. The developed models were validated by devised numerical experiments; and the optimized layouts using the proposed algorithm were compared to those by the Non-Dominated Sorting Genetic Algorithm (NSGA-II). Our results provide some evidence regarding the inability of the NSGA-II to cover the objective space when optimizing wind farms with large power densities. In contrast, our proposed approach succeeds in obtaining high-density layouts efficiently. Furthermore, we demonstrated the superiority of the developed algorithm, in the aspects of coverage, spread, and computational cost. (C) 2017 Elsevier Ltd. All rights reserved.