  We present forecasts for constraints on the Hu \& 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_{\rm m}, \Omega_\nu h^2, h^2, 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-3G$\times$ngWL 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-S4$\times$ngWL 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}| \geq -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_{\mathrm{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.
