Abstract: Variable-density groundwater flow (VDGF) is jointly driven by hydraulic and density gradient, leading to strong nonlinearity, large computational burden of numerical models, and therefore huge computational cost of Monte Carlo simulation for uncertainty analysis. This study developed the reduced-order model (ROM) for VDGF and built the Gaussian process (GP) for simulating the numerical error of the ROM. The coupled model can obtain solutions of head and salinity across the study domain while GP simulates observation information at limited locations. Moreover, the coupled model can provide higher solution accuracies of head and salinity at the observation locations than the ROM. A two-dimensional (cross-section) VDGF test case was considered, where hydraulic conductivity was taken as a spatially random field. MC simulations were performed using three models, including the full-system model, the ROM, and the coupled model, with corresponding MC strategies denoted as FSMC, ROMC, and GP-ROMC, respectively. The results show that ROMC can be an alternative to FSMC for conducting uncertainty quantification. The relationship between head (or salinity) and the dimensional of ROM can be characterized using power functions with determinate coefficients larger than 0.99. GP-ROMC has higher solution accuracy than ROMC, which indicates that GP is capable for simulating the numerical error of ROM. The results in this study are significant for performing simulation, uncertainty quantification, risk assessment, and parameter estimate in the context of groundwater.