Density-driven groundwater flows are described by nonlinear coupled differential equations. Due to its importance in engineering and earth science, several linearizations and semi-linearization schemes for approximating their solution have been proposed. Among the more efficient are the combinations of Newtonian iterations for the spatially discretized system obtained by either scalar homotopy methods, fictitious time methods, or meshless generalized finite difference method, with several implicit methods for the time integration. However, when these methods are used, some parameters need to be determined, in some cases, even manually. To overcome this problem, this paper presents a novel generalized finite differences scheme combined with an adaptive step-size method, which can be applied for solving the governing equations of interest on non-rectangular structured and unstructured grids. The proposed method is tested on the Henry and the Elder problems to verify the accuracy and the stability of the proposed numerical scheme.