Inferring CO2 saturation from seismic data is important when seismic methods are applied at CO2 sequestration sites for verification and accounting purposes, such as verifying the total injected CO2 volume, comparing with model predictions for concordance evaluation, tracking the migration of CO2 plume, and detecting possible leakage from the storage reservoir. In this work, we infer CO2 saturation levels at three depths from simulated surface seismic, downhole pressure and total dissolved solids (TDS) data using machine learning (ML) methods. The simulated monitoring data are based on 6000 numerical multi-phase flow simulations of hypothetical wellbore CO2 and brine leakage from a legacy well into shallow aquifers at a model CO2 storage site. We conduct rock physics modeling to estimate changes in seismic velocity due to the simulated CO2 and brine leakage at each time step in the flow simulation outputs, resulting in 120,000 forward seismic velocity models. 2D finite-difference acoustic wave modeling is performed for each velocity model to generate synthetic shot gathers, along a sparse 2D seismic line with only 5 shots and 40 receivers. We extract 6 time-lapse seismic attribute anomalies from each trace in the time window relevant to each geologic layer, and use the seismic features, together with downhole pore pressure, TDS features to train the machine learning algorithms. The impact of seismic noise on the performance of the trained machine learning models has also been investigated. Inferred CO2 saturations from the trained classifiers are in good agreement with observations. Direct pressure and TDS measurements from downhole monitoring can increase the accuracy of the inferred CO2 saturation class from the forward modeled 2D surface seismic data. Our ML workflow represents a promising way to combine measurements from multiple monitoring techniques, together with seismic monitoring to achieve more accurate seismic quantitative interpretation.