  Modified gravity (MG) theories predict, in general, that the ratio of
gravitational wave (GW) to electromagnetic (EM) luminosity distances, $\Xi$,
differs from its general relativity (GR) value of unity at cosmological scales,
thus providing another perturbative probe to MG. In this paper, we introduce
new phenomenological parametrizations for both the
Friedmann-Lema\^itre-Robertson-Walker (FLRW) background evolution of $f(R)$
models, via the dark energy equation of state parameter, $w_{\text{DE}}$, and
for $\Xi$ in this class of theories. We simulate a mock dataset for the
Einstein Telescope (ET) of 1000 GW signals from binary neutron star (BNS)
mergers and redshift information from their EM counterpart, exploring the
consequent constraints on the relevant gravitational, cosmological and
phenomenological parameters. As a model of particular interest, we take
$\gamma$-gravity theory and investigate whether it could be distinguished from
GR's $\Lambda$CDM model. We then combine our results with actual data from type
Ia supernovae (SNIa) and combined baryon acoustic oscillations (BAO) and cosmic
microwave background (CMB) observations. Additionally, we also investigate the
potential bounds to $f_{R0}$ for any viable $f(R)$ whose background evolution
is indistinguishable from the standard model of cosmology above a certain
redshift, showing that, for a $\Lambda$CDM fiducial model, ET data would
provide $|f_{R0}|<10^{-2}$ at a 95\% level. We conclude altogether that probing
the redshift evolution of the GW luminosity distance from detections of the ET
in its first running decade will not substantially help constraining $f(R)$
theories of gravity.
