Monitoring groundwater (GW) level variations, or anomalies in multiple wells, over long periods of time is essential to understanding changes in regional groundwater resource availability. However, it is challenging to predict these GW anomalies over the long term in agricultural areas due to complicated boundary conditions, heterogeneous hydrogeological characteristics, and groundwater extraction, as well as nonlinear interactions among these factors. To overcome this challenge, we developed an advanced modeling framework based on a recurrent neural network of long short-term memory (LSTM) as an alternative to complex and computationally expensive physical models. GW anomalies were forecast two months in advance (t + 2) based on the evaluation of drivers that influence GW dynamics in densely irrigated agricultural regions. An application and evaluation of this new approach was conducted in the Wisconsin Central Sands (WCS) region in the U.S., one of the most productive agricultural regions. The modeling framework was developed for the period 1958–2020 by utilizing easily accessible dynamic and static variables to represent hydrometeorological and geological characteristics. GW anomaly observations were acquired from 26 piezometers (wells) installed in the sandy aquifer in WCS over 10–60 years. The subset of GW observations from ∼ 10–60 years, not used in model training, can forecast GW anomalies two months out with a coefficient of determination R2 ∼ 0.8. Additionally, MAE was less than 0.34 m/month across the study region for both temporal and spatial modeling frameworks. Groundwater anomalies showed high spatiotemporal variability, and their responses are influenced differently by boundary conditions, catchment geology, climate, and topography across locations. Sites with higher autocorrelation with previous two-months GW anomalies reduced bias by increasing R2. Land use change and irrigation pumping have interactive effects on GW anomalies forecasting. The novelty of this study is identifying the regional drivers of GW fluxes. This case-specific information and location-related simplification, modification, and assumption of LSTM is a unique contribution to the existing literature. Our framework can be used as an alternative method of simulating water availability and groundwater level changes in areas where subsurface properties are unknown or difficult to determine.