This paper compares the results produced by MODFLOW, a constant-density model, to results produced by SEAWAT, a variable-density model, to investigate the feasibility of using MODFLOW in a saline environment below an estuary known as the Indian River Lagoon. The comparison was conducted over sixteen numerical simulation cases at different conditions of estuarine salinity CL, hydraulic conductivity anisotropy ratio Kr, and water table elevations on the freshwater boundaries in a two-dimensional vertical domain. The use of MODFLOW at the study site under the calibrated Kr distribution ranging from 1000-20,000 was found to accurately match the field-measured and SEAWATsimulated results with a remarkable increase in accuracy at higher groundwater elevations. The study determined a critical value of Kr of 1000 above which, MODFLOW simulations of the variable-density problem produced results that agreed well with those produced by SEAWAT. However, MODFLOW starts to produce significant errors with Kr below the critical value and hence, it should not be used for simulating variable-density environments when Kr<1000. The amount of submarine groundwater discharge (SGD) predicted by either model, and also MODFLOW accuracy in predicting the SGD are directly proportional to the head difference between the groundwater divide elevation and the lagoon water surface, but to a lower extent, are inversely proportional to CL.