Pump-and-treat (P&T) is a widely-adopted solution for the containment of solute plumes in contaminated aquifers. A cost-effective design of P&T systems requires optimizing (minimizing) the overall pumping rates (Q). This optimization is a stochastic process, as Q is a random variable linked to the randomness of the aquifer hydraulic conductivity (K). Previously presented stochastic approaches to minimize Q adopted two-dimensional (2D) Gaussian random spatial fields (r.s.f.) of log-transformed K. Recent studies based on geological entropy have demonstrated the limited ability of Gaussian r.s.f. to reproduce extreme K patterns, which mostly control transport in heterogeneous aquifers, when compared to non-Gaussian r.s.f. Moreover, 2D models generate different flow and transport connectivity than three-dimensional (3D) models. On these premises, this work aimed at extending previous works on P&T optimization in heterogeneous aquifers through Monte-Carlo groundwater simulations of 2D and 3D Gaussian and non-Gaussian r.s.f. The results indicated that the mean (Q<overbar></mml:mover>n) and variance (sigma Qn2</mml:msubsup>) of the optimal Q distribution depend strictly on the chosen model dimensionality and r.s.f. generator. In particular, 2D models and models embedding indicator-based (i.e. non-Gaussian) r.s.f. tended to generate higher Q<overbar></mml:mover>n and sigma Qn2</mml:msubsup> than 3D models with increasing number of model layers (K-L) and Gaussian models. This behavior can be explained considering the spatial ordering of K clusters in the simulated aquifers, which is measured through metrics derived from the concept of geological entropy. It was found that 2D models and models embedding non-Gaussian r.s.f. displayed more spatially-persistent ordered K structures than 3D models and Gaussian models, resulting in higher Q<mml:mo stretchy="false"><overbar></mml:mover>n and sigma Qn2. This is attributed to the relative amount of heterogeneity sampled by the solute source and the increased likelihood of more ordered K clusters to generate preferential flow and solute transport channeling than more disordered and chaotic systems, which enhance solute mixing. Combining P&T with physical barriers (i.e. cut-off walls) was helpful to reduce both <mml:msub><mml:mover accent="true">Q<mml:mo stretchy="false"><overbar></mml:mover>n and sigma Qn2 in all tested scenarios, corroborating previous findings. However, the relative efficacy of a specific physical barrier geometry to reduce <mml:msub><mml:mover accent="true">Q<mml:mo stretchy="false"><overbar></mml:mover>n and sigma Qn2 also depends on the chosen model dimensionality and r.s.f. generator.
