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.　©　2024　Editorial　Office　of　Hydrogeology　and　Engineering.　All　rights　reserved.