The mass discharge emanating from dense non‐aqueous phase liquid (DNAPL) source zones (SZs) is often used as a key metric for risk assessment. To predict the temporal evolution of mass discharge, upscaled models have been developed to approximate the relationship between the depletion of SZ and the mass discharge. A significant challenge stems from the choice of the SZ parameterization, so that a limited number of domain‐averaged SZ metrics can suffice as an input and accurately predict the complex mass‐discharge behavior. Moreover, existing deterministic upscaled models cannot quantify prediction uncertainty stemming from modeling parameterization. To address these challenges, we propose a method based on a Bayesian Neural Network (BNN) which learns the non‐linear relationship between SZ metrics and mass discharge from multiphase‐modeling training data. The proposed BNN‐based upscaled model allows uncertainty quantification since it treats trainable parameters as distributions, and does not require a manual parameterization of the SZ a‐priori. Instead, the BNN model chooses three physically meaningful SZ quantities related to mass discharge as input features. Then, we use the expected gradients method to identify the feature importance for mass‐discharge prediction. We evaluated the proposed model on laboratory‐scale DNAPL dissolution experiments. The results show that the BNN model accurately reproduces the multistage mass‐discharge profiles with fewer parameters than existing upscaled models. Feature importance analysis shows that all chosen features are important and sufficient to reproduce complex mass discharge. This model provides accurate mass‐discharge predictions and uncertainty estimation, therefore holds a great potential for probabilistic risk assessments and decision‐making.