History matching is a critical process used for calibrating simulation models and assessing subsurface uncertainties. This common technique aims to align the reservoir models with the observed data. However, achieving this goal is often challenging due to the non uniqueness of the solution, underlying subsurface uncertainties, and usually the high computational cost of simulations. The traditional approach is often based on trial and error, which is exhaustive and labor-intensive. Some analytical and numerical proxies combined with Monte Carlo simulations are utilized to reduce the computational time. However, these approaches suffer from low accuracy and may not fully capture subsurface uncertainties. This study proposes a new robust method utilizing Bayesian Markov Chain Monte Carlo (MCMC) to perform assisted history matching under uncertainties. We propose a novel three-step workflow that includes 1) multi-resolution low-fidelity models to guarantee high-quality matching; 2) Long-Short Term Memory (LSTM) network as a low-fidelity model to reproduce continuous time-response based on the simulation model, combined with Bayesian optimization to obtain the optimum low fidelity model; 3) Bayesian MCMC runs to obtain the Bayesian inversion of the uncertainty parameters. We perform sensitivity analysis on the LSTM’s architecture, hyperparameters, training set, number of chains, and chain length to obtain the optimum setup for Bayesian-LSTM history matching. We also compare the performance of predicting the recovery factor using different surrogate methods, including polynomial chaos expansions (PCE), kriging, and support vector machines for regression (SVR). We demonstrate the proposed method using a water flooding problem for the upper Tarbert formation of the tenth SPE comparative model. This study case represents a highly heterogeneous nearshore environment. Results showed that the Bayesian-optimized LSTM has successfully captured the physics in the high-fidelity model. The Bayesian-LSTM MCMC produces an accurate prediction with narrow ranges of uncertainties. The posterior prediction through the high-fidelity model ensures the robustness and accuracy of the workflow. This approach provides an efficient and practical history-matching method for reservoir simulation and subsurface flow modeling with significant uncertainties.