Groundwater level time series are increasingly simulated with both physically motivated transfer function-noise (TFN) models and purely data-driven machine learning (ML) models, but the two families are rarely compared under an identical, forcing-only protocol. This study compared a TFN model (Pastas, gamma impulse response, linear recharge, AR(1) noise) with random forest (RF), extreme gradient boosting (XGBoost), support vector regression (SVR) and a long short-term memory (LSTM) network for simulating monthly mean hydraulic heads in three observation wells in the Netherlands (B58C0698 near Heibloem; B32C0609 and B32C0639 near De Bilt), with records of 24–37 years. All models were driven only by precipitation and reference evaporation, without lagged heads, and were trained on the first 70–80% of each record and tested on the remaining 79–95 months. Heads showed a regular winter–spring maximum and late-summer minimum and were most strongly correlated with the 6-month cumulative precipitation surplus (r = 0.74–0.93). No significant trend was found in annual heads or annual precipitation surplus. In the test periods, all models performed well at B58C0698 (NSE = 0.90–0.93) and B32C0639 (NSE = 0.79–0.87), with differences between models that were not statistically significant (Diebold–Mariano p > 0.27). At B32C0609, whose test period coincided with the start of daily logger measurements, all models underestimated heads by 8–12 cm and NSE fell to 0.36–0.58. The tree ensembles were the least affected but not significantly better than TFN (p ≥ 0.05). Averaged over the three wells, SVR (NSE = 0.76) and RF (0.76) ranked slightly above XGBoost (0.75), LSTM (0.71) and TFN (0.70). Permutation importance attributed 54–91% of the ML skill to the 6-month precipitation surplus, and the TFN step responses reached 95% of their final value after 323–441 days. With forcing-only inputs and a few hundred monthly observations, a parsimonious TFN model is as accurate as the ML models and easier to interpret, whereas the ML models add value mainly where the response is not captured by a linear recharge assumption.