Coupling Bayes’ Theorem with a two-dimensional (2D) groundwater solute advection-diffusion transport equation allows an inverse model to be established to identify a set of contamination source parameters including source intensity (M ), release location ( X0 , Y0) and release time (T0), based on monitoring well data. To address the issues of insufficient monitoring wells or weak correlation between monitoring data and model parameters, a monitoring well design optimization approach was developed based on the Bayesian formula and information entropy. To demonstrate how the model works, an exemplar problem with an instantaneous release of a contaminant in a confined groundwater aquifer was employed. The information entropy of the model parameters posterior distribution was used as a criterion to evaluate the monitoring data quantity index. The optimal monitoring well position and monitoring frequency were solved by the two-step Monte Carlo method and differential evolution algorithm given a known well monitoring locations and monitoring events. Based on the optimized monitoring well position and sampling frequency, the contamination source was identified by an improved Metropolis algorithm using the Latin hypercube sampling approach. The case study results show that the following parameters were obtained: 1) the optimal monitoring well position (D) is at (445, 200); and 2) the optimal monitoring frequency (Δt) is 7, providing that the monitoring events is set as 5 times. Employing the optimized monitoring well position and frequency, the mean errors of inverse modeling results in source parameters (M, X0 ,Y0 ,T0 ) were 9.20%, 0.25%, 0.0061%, and 0.33%, respectively. The optimized monitoring well position and sampling frequency can effectively safeguard the inverse modeling results in identifying the contamination source parameters. It was also learnt that the improved Metropolis-Hastings algorithm (a Markov chain Monte Carlo method) can make the inverse modeling result independent of the initial sampling points and achieves an overall optimization, which significantly improved the accuracy and numerical stability of the inverse modeling results.