Groundwater management practices, such as sustainable drinking water extraction or contamination protection, can have significant socioeconomic impacts. A complete uncertainty analysis should ideally be performed to anticipate all possible outcomes and assess any risk. Uncertainties arise as a result of our limited understanding of the physical processes involved, as well as a scarcity of measurement data, whether directly or indirectly related to the physical parameters of interest. In this paper, we use a small number of tracing experiments (predictor) to predict the wellhead protection area (WHPA, target), the shape and extent of which are influenced by the distribution of hydraulic conductivity (K). Our first goal is to make stochastic predictions of the WHPA within the Bayesian Evidential Learning (BEL) framework, which uses machine learning to find a direct relationship between predictor and target. This relationship is learned using a small number of training models (400) drawn from the prior distribution of K. Forward modelling is used to obtain the 400 pairs of simulated predictors and targets. Newly collected field data can then be used directly to predict the approximate posterior distribution of the corresponding WHPA, obviating the need for the traditional step of data inversion. The number and location of data sources (injection wells) influence the posterior WHPA distribution's uncertainty range. Our second goal is to extend BEL to determine the optimal design of data source locations that minimises the WHPA's posterior uncertainty. Because the BEL model, once trained, allows the computation of the posterior uncertainty corresponding to any new input data, experimental design can be done explicitly, without averaging or approximating. We estimate the WHPA's posterior uncertainty range using the Modified Hausdorff Distance (MHD) and Structural Similarity (SSIM) index metrics. Because the breakthrough curves store information on a large area of the K field surrounding the pumping well, increasing the number of injection wells reduces the derived posterior WHPA uncertainty. Our method can also estimate which injection wells are more informative than others, as demonstrated by a k-fold cross-validation procedure. Overall, the application of BEL to experimental design allows for identifying data sources that maximise the information content of any measurement data while keeping budget constraints and computational costs to a minimum.