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 log10|fR0|. In this analysis we vary the cosmological parameters [Ωm,Ωνh2,h,As,ns,log10|fR0|], 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−3G×ngWL is expected to obtain an upper limit of log10|fR0|<−5.95 at 95% credibility, which significantly improves upon the current best bounds. The CMB−S4×ngWL dataset is expected to improve this even further to log10|fR0|<−6.23. Furthermore, f(R) gravity models with log10|fR0|≥−6, which have larger numbers of clusters, would be distinguishable from GR with both datasets. We also report degeneracies between log10|fR0| and Ωm as well as σ8 for log10|fR0|>−6 and log10|fR0|>−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. Published by the American Physical Society 2024