In the current study, the hieratical zeolitic imidazole framework 8 (H/ZIF-8) and bimetallic H/ZIF-8, H/ZIF-8@La and H/ZIF-8@Cu, were prepared with a simple and green method using water as a solvent. Process parameters (temperature and time of crystallization), as well as compositional parameters (amount of 2-Methylimidazole and 2-(methylamino)ethanol), were studied. H/ZIF-8 samples were characterized with XRD, EDS, FE-SEM, TGA/DTA, FTIR, N-2 isotherms, and their performance was evaluated for sulfate uptake from aqueous media. An experimental design was utilized in this study to optimize the independent variables using central composite design (CCD) under the response surface methodology (RSM) method. A significant agreement between the models and experimental data was verified by analysis of variance (ANOVA). The adsorption equilibrium models of Langmuir, Jovanovic, Freundlich, and Temkin isotherms were evaluated and the results described that the Langmuir model was the best with the experimental data. Different kinetic models were estimated and found that pseudo-second-order kinetic data were well-fitted for removal reaction. The determination of different thermodynamic parameters reflected that the sulfate uptake was spontaneous and feasible and that three H/ZIF-8 samples had an endothermic nature. The adsorption on H/ZIF-8 samples was not significantly influenced by the competing anions of nitrate, chloride and fluoride but phosphate displayed slightly greater negative effects. The used H/ZIF-8 samples could be regenerated and reused in eight consecutive cycles with a proper desorption agent. The results of sulfate removal from a real sample revealed that using H/ZIF-8 and two bimetallic H/ZIF-8 samples for sulfate removal from polluted waters is a promising alternative for sulfate recovery.