AbstractOptimizing the remediation of light nonaqueous phase liquids (LNAPLs) to achieve an acceptable endpoint status for a site is not trivial. Recently, Sookhak Lari, Johnston, et al. (2018, https://doi.org/10.1016/j.jhazmat.2017.11.006), Sookhak Lari, Rayner, and Davis (2018, https://doi.org/10.1016/j.jenvman.2018.07.041) conducted three‐dimensional multiphase, multicomponent simulations to address LNAPL remediation endpoints for a single recovery well. However, optimized LNAPL remediation for multiple wells is not addressed in the literature. In the first part of this paper, we establish a matrix of 10 simulation scenarios to show the sensitivity of the remedial endpoint (i.e., what is feasibly achieved) to several parameters including viscosity and partitioning attributes of the LNAPL, heterogeneity of the formation, and the location and number of the recovery wells. While this addresses the variability of LNAPL removal from the subsurface and is valuable in its own right, it does not address the optimal removal of LNAPL. We address this in the second part of the paper by linking a genetic algorithm to TMVOC‐MP to allow, for example, the assessment of the optimal number and location of LNAPL recovery wells in a field‐scale problem. Using supercomputing facilities and within 49 genetic algorithm generations, each including 150 members, highly optimized answers to different objective functions were obtained. For the first time, such a multiphase, multicomponent optimization tool promises the possibility of optimizing LNAPL remediation at field scale to achieve practicable endpoint conditions, within computational affordability.