We　develop　a　reduced-order　modeling　strategy　aimed　at　providing　numerical　Monte　Carlo　simulations　of　groundwater　flow　in　randomly　heterogeneous　transmissivity　fields.　We　rely　on　moment　equations　for　groundwater　flow　and　conduct　space　reductions　for　both　transmissivity,　T,　and　hydraulic　head,　h.　A　truncated　singular　value　decomposition　(SVD)　solver　is　employed　to　cope　with　the　ill-conditioned　stiffness　matrix　caused　by　(negative　and　thus)　unphysical　values　of　T　that　might　arise　due　to　possible　low　accuracy　stemming　from　the　order　of　model　reduction.　The　performance　of　the　approach　is　assessed　through　the　analysis　of　various　synthetic　reference　scenarios.　These　encompass　diverse　degrees　of　heterogeneity　of　the　transmissivity　field　and　various　values　of　reduced-order　dimensions,　n　and　m,　associated　with　h　and　T,　respectively.　Transmissivity　is　conceptualized　as　a　composite　(spatial)　random　field　where　there　is　uncertainty　in　the　locations　of　regions　associated　with　diverse　geomaterials　as　well　as　in　the　heterogeneity　of　transmissivity　therein.　Our　results　are　also　compared　against　their　counterparts　that　one　could　obtain　upon　performing　a　model　reduction　solely　on　the　basis　of　hydraulic　heads.　Our　findings　show　that:　(i)　resting　on　the　truncated　SVD　solver　is　beneficial　for　coping　with　ill-conditioned　stiffness　matrices;　(ii)　the　two　model　reduction　strategies　provide　comparable　solution　accuracy　for　m　≥　5n,　while　(iii)　the　computational　cost　associated　with　the　reduced-order　model　based　on　space　reduction　for　both　T　and　h　is　always　significantly　smaller　than　that　associated　with　space　reduction　based　solely　on　h.　©　2024　Elsevier　B.V.