  This paper investigates the impact of the effects of bulk viscosity on the
accelerating expansion and large-scale structure formation of a universe in
which the underlying gravitational interaction is described by $f(Q)$ gravity.
Various toy models of the $f(Q)$ gravity theory, including power-law
($f_{1\rm{CDM}}$), exponential ($f_{2\rm{CDM}}$), and logarithmic
($f_{3\rm{CDM}}$) are considered. To test the cosmological viability of these
models, we use 57 Hubble parameter data points (OHD), 1048 supernovae distance
modulus (SNIa), their combined analysis (OHD+SNIa), 14 growth rate ($f$-data),
and 30 redshift-space distortions ($f\sigma_8$) datasets. We compute the
best-fit values $\Omega_m, H_0$ and exponents $\{n, p, \gamma\}$ including the
bulk viscosity coefficient $\zeta_0$, through a detailed statistical analysis.
Moreover, we study linear cosmological perturbations and compute the density
contrast \(\delta(z)\), growth factor $D(z)$, growth rate $(f(z)$, and
redshift-space distortion $f\sigma_8(z)$. Based on the Akaike Information
Criterion (AIC) and Bayesian / Schwartz Information Criterion (BIC), a
statistical comparison of the $f(Q)$ gravity models with \lcdm is made. From
our statistical analysis of cosmic measurements, we found an underestimation of
all models on the OHD data; therefore, statistical viability led us to weigh it
in favour of SNIa data which resulted in the exponential ($f_{2\rm{CDM}}$)
model without bulk viscosity being the most plausible alternative model, while
on average the $f_{1\rm{CDM}}$ model performed statistically the weakest. We
also found that adding bulk viscosity did not improve the fit of the models to
all observational data.
