  The parameters in the governing system of partial differential equations of
multicompartmental poroelastic models typically vary over several orders of
magnitude making its stable discretization and efficient solution a challenging
task. In this paper, inspired by the approach recently presented by Hong and
Kraus~[Parameter-robust stability of classical three-field formulation of
Biot's consolidation model, ETNA (to appear)] for the Biot model, we prove the
uniform stability, and design stable disretizations and parameter-robust
preconditioners for flux-based formulations of multiple-network poroelastic
systems. Novel parameter-matrix-dependent norms that provide the key for
establishing uniform inf-sup stability of the continuous problem are
introduced. As a result, the stability estimates presented here are uniform not
only with respect to the Lam\'e parameter $\lambda$, but also with respect to
all the other model parameters such as permeability coefficients $K_i$, storage
coefficients $c_{p_i}$, network transfer coefficients $\beta_{ij},
i,j=1,\cdots,n$, the scale of the networks $n$ and the time step size $\tau$.
  Moreover, strongly mass conservative discretizations that meet the required
conditions for parameter-robust stability are suggested and corresponding
optimal error estimates proved. The transfer of the canonical (norm-equivalent)
operator preconditioners from the continuous to the discrete level lays the
foundation for optimal and fully robust iterative solution methods. The
theoretical results are confirmed in numerical experiments that are motivated
by practical applications.
