Galaxy clusters are the largest virialised objects in the Universe and, having formed from the highest and rarest peaks in the primordial density field, their abundance is sensitive to various cosmological parameters as well as to the underlying theory of gravity. They have been used to constrain and forecast the constraints on modified gravity models, with initial results looking promising, especially when considering the observational advancements in the coming years. However, the comparison between observed cluster abundance and the predicted dark matter halo abundance from simulations is complicated by the fact that cluster mass is not directly observable, but has to be inferred from other observations. These cluster observable-mass scaling relations are usually modified in modified gravity models and, if not accounted for, could lead to biased model constraints. In this presentation I will describe a recent effort to deal with this issue and to ensure unbiased model tests, and present some forecasts of the constraining power in f(R) gravity using cluster SZ observations. Working with two leading classes of modified gravity models - f(R) gravity and the DGP braneworld model -- I will also show that, while the law of gravity is modified in complex ways in these models, the statistical properties of dark matter haloes can be very well described by simple models. These latter models are part of the constraint pipeline to be described in this talk.