In this study, we present a machine learning framework utilizing Temporal Fusion Transformers (TFTs) to predict ionospheric variability, which is difficult to predict due to the nonlinear coupling between solar, geomagnetic, and thermospheric drivers. We predict the total electron number density (TEC) derived from GNSS observations, using a variety of inputs including solar radiation, geomagnetic indices, and GNSS-based vertical TEC. Experimental results from 2010 to 2025 demonstrate a low root mean square error (RMSE) of 3.33 TECU for predictions up to 24 hours in advance, with solar EUV radiation providing the strongest predictive signal. Reproducibility is enhanced through the open-source toolkit \texttt{ionopy}.