We present forecasts for constraints on the Hu and Sawicki f(R) modified gravity model using realistic mock data representative of future cluster and weak lensing surveys. We create mock thermal Sunyaev-Zel'dovich effect selected cluster samples for SPT-3G and CMB-S4 and the corresponding weak gravitational lensing data from next-generation weak-lensing (ngWL) surveys like Euclid and Rubin. We employ a state-of-the-art Bayesian likelihood approach that includes all observational effects and systematic uncertainties to obtain constraints on the f(R) gravity parameter log(10)|f(R0)|. In this analysis we vary the cosmological parameters [Omega(m),Omega(nu)h(2),h,A(s),n(s),log(10)|f(R0)|], which allows us to account for possible degeneracies between cosmological parameters and f(R) modified gravity. The analysis accounts for f(R) gravity via its effect on the halo mass function which is enhanced on cluster mass scales compared to the expectations within general relativity (GR). Assuming a fiducial GR model, the upcoming cluster dataset SPT-3GxngWL is expected to obtain an upper limit of log(10)|f(R0)|<-5.95 at 95% credibility, which significantly improves upon the current best bounds. The CMB-S4xngWL dataset is expected to improve this even further to log(10)|f(R0)|<-6.23. Furthermore, f(R) gravity models with log(10)|f(R0)|>=-6, which have larger numbers of clusters, would be distinguishable from GR with both datasets. We also report degeneracies between log(10)|f(R0)| and Omega(m) as well as sigma(8) for log(10)|f(R0)|>-6 and log(10)|f(R0)|>-5 respectively. Our forecasts indicate that future cluster abundance studies of f(R) gravity will enable substantially improved constraints that are competitive with other cosmological probes.