The natural Egyptian bentonite, collected from south El Hammam area, was modified at three different temperatures 100°C, 200°C and 300°C for 1 h. The raw and modified bentonite samples were characterized by powder X-ray diffraction (XRD), scanning electron microscope (SEM) and BET surface area. The bentonite modified at 100°C exhibited more flaky grains with smooth surface and high surface area as compared to the two other modified types. Response surface methodology in conjunction with central composite rotatable design was used in optimizing and modeling the effect of different parameters such as contact time, initial concentration and dose on the removal of iron ions. Second order quadratic polynomial model was selected to represent the removal process. The mathematical equations of quadratic polynomial model were derived from Design Expert Software (Version 6.0.5). The predicted values from the mathematical equations were highly correlated with the experimental results (R 2 above 0.9) for the required responses in untreated and modified bentonite at 100°C for 1 h. 3D and linear graphs were used to understand the effect of the studied variable parameters and the interaction between them. Under the predicted conditions suggested by the quadratic programming, the modified bentonite at 100°C is more promising and the removal efficiency could be enhanced to 100%. The quadratic polynomial model could be efficiently applied for the modeling of iron removal from aqueous solutions by bentonite.