Summary Subsurface development involves well-placement decisions considering the highly uncertain understanding of the reservoir in the subsurface. The simultaneous optimization of a large number of well locations is a challenging problem. Traditional gradient-based methods can be adapted for well location optimization (WLO) when these problems are converted into real-valued representations and equipped with protocols to handle noisy objective functions. However, their application to large-scale scenarios often remains impractical. This impracticality arises because computing gradients of the objective function can be prohibitively expensive in realistic settings, particularly without using the adjoint method. In this paper, we explore the application of a novel quasi-Newton trust-region (TR) method that employs the stochastic simplex approximate gradient (StoSAG). We have implemented the Broyden-Fletcher-Goldfarb-Shanno (BFGS) quasi-Newton updating algorithm together with line-search (LS) and TR optimization strategies. The StoSAG-based optimization methods have been applied to a realistic synthetic reservoir featuring 26 wells considering two distinct cases: Each includes 20 realizations of porosity and permeability. The first case exhibits mild heterogeneity, while the second exhibits significant heterogeneity with a large correlation length. We have conducted a series of runs to evaluate the performance of these algorithms in addition to comparisons to the finite-difference (FD) and particle-swarm-optimization (PSO) algorithms. We introduce a novel approach to enhance the accuracy of StoSAG gradients by proposing modified StoSAG formulations. These formulations are tailored to exploit the structure of the objective function and to capture the relationships between its components and the individual optimization parameters. This approach involves using a correction matrix W informed by problem-specific knowledge. The entries of W vary from 0 to 1 and are proportional to the interference effects the neighboring wells have on each other. We determine those entries (or weights) based on the radii of investigation around the wells and the distance between the well pairs. Results indicate that the steepest-descent (SD) algorithm coupled with StoSAG has superior performance to PSO and FD. Although the objective function is prone to numerical noise and not continuously differentiable with respect to well locations, StoSAG overcomes this challenge because it acts as a smooth approximation. Comparative tests further confirm that the TR-BFGS is more effective than the LS-BFGS. Moreover, we show that using the proposed gradient correction procedure results in a significant acceleration in convergence, indicating an enhancement in the StoSAG gradient approximation quality. This enhancement allows the TR-BFGS algorithm to achieve considerably higher performance than SD, illustrating that the accuracy of the BFGS approximation benefits from improved gradient quality.