Gravitational wave (GW) standard sirens are well-established probes with which one can measure cosmological parameters, and are complementary to other probes like the cosmic microwave background or supernovae standard candles. I will focus on dark GW sirens, specifically binary black holes (BBHs) for which there is only GW data. Relying on the assumption of a source mass model for the BBH distribution, we consider four models that are representative of the BBH population observed so far. In addition to inferring cosmological and mass model parameters, we use dark sirens to test modified gravity theories. These theories often predict a different GW propagation, leading to a changed GW luminosity distance. We perform a joint estimate of the population parameters governing mass, redshift, cosmology, and GW propagation. We use data from the third LIGO-Virgo-KAGRA observation run (O3) and find — for the four mass models and three signal-to-noise cutoffs — that GR is consistently the preferred model to describe all observed BBH GW signals to date. Furthermore, all modified gravity parameters have posteriors that are compatible with the values predicted by GR at the $90\%$ confidence interval. We then focus on future observation runs O4 and O5. We show that there are strong correlations between cosmological, astrophysical and modified gravity parameters. If GR is the correct theory of gravity, and assuming narrow priors on the cosmological parameters, we recover the modified gravity parameter $\Xi_0 = 1.47^{+0.92}_{-0.57} $ with O4, and $\Xi_0 = 1.08 ^{+0.27}_{-0.16}$ with O4 and O5 (General relativity (GR) corresponds to $\Xi_0= 1$.). If, however, Nature follows a modified gravity model with $\Xi_0=1.8$ and $n=1.91$, we exclude GR at the 1.7 sigma level with O4 and at the 2.3 sigma level with O4 and O5 combined.