Sea water intrusion (SWI) simulators are essential tools to assist the sustainable management of coastal aquifers. These simulators require the solution of coupled variable-density partial differential equations (PDEs), which reproduce the processes of groundwater flow and dissolved salt transport. The solution of these PDEs is typically addressed numerically with the use of density-dependent flow simulators, which are computationally intensive in most practical applications. To this end, model surrogates are generally developed as substitutes for full-scale aquifer models to trade off accuracy in exchange for computational efficiency. Surrogates represent an attractive option to support groundwater management situations in which fast simulators are required to evaluate large sets of alternative pumping strategies. Reduced-order models, a sub-category of surrogate models, are based on the original model equations and may provide quite accurate results at a small fraction of computational cost. In this study, a variable-density flow reduced-order model based on proper orthogonal decomposition (POD) and utilizing a fully coupled flow and solute-transport model is implemented with a finite-difference (FD) approach for simulating SWI in coastal aquifers. The accuracy and computational efficiency of the FD-POD approach for both homogeneous and-more realistic-heterogeneous systems are investigated using test cases based on the classic Henry's problem (Henry 1964). The findings demonstrate that the combined FD-POD approach is effective in terms of both accuracy and computational gain and can accommodate the output of the most popular variable-density flow models, such as those from USGS's MODFLOW family.