This work presents a computational formulation of the Bayesian maximum entropy (BME) approach to solve a stochastic partial differential equation (PDE) representing the advection-reaction process across space and time. The solution approach provided by BME has some important features that distinguish it from most standard stochastic PDE techniques. In addition to the physical law, the BME solution can assimilate other sources of general and site-specific knowledge, including multiple-point nonlinear space/time statistics, hard measurements, and various forms of uncertain (soft) information. There is no need to explicitly solve the moment equations of the advection-reaction law since BME allows the information contained in them to consolidate within the general knowledge base at the structural (prior) stage of the analysis. No restrictions are posed on the shape of the underlying probability distributions or the space/time pattern of the contaminant process. Solutions of nonlinear systems of equations are obtained in four space/time dimensions and efficient computational schemes are introduced to cope with complexity. The BME solution at the prior stage is in excellent agreement with the exact analytical solution obtained in a controlled environment for comparison purposes. The prior solution is further improved at the integration (posterior) BME stage by assimilating uncertain information at the data points as well as at the solution grid nodes themselves, thus leading to the final solution of the advection-reaction law in the form of the probability distribution of possible concentration values at each space/time grid node. This is the most complete way of describing a stochastic solution and provides considerable flexibility concerning the choice of the concentration realization that is more representative of the physical situation. Numerical experiments demonstrated a high solution accuracy of the computational BME approach. The BME approach can benefit from the use of parallel processing (the relevant systems of equations can be processed simultaneously at each grid node and multiple integrals calculations can be accelerated significantly, etc.).