  A noteworthy aspect in blood flow modeling is the definition of the
mechanical interaction between the fluid flow and the biological structure that
contains it, namely the vessel wall. It has been demonstrated that the addition
of a viscous contribution to the mechanical characterization of vessels brings
positive results when compared to in-vivo measurements. In this context, the
implementation of boundary conditions able to keep memory of the viscoelastic
contribution of vessel walls assumes an important role, especially when dealing
with large circulatory systems. In this work, viscoelasticity is taken into
account in entire networks via the Standard Linear Solid Model. The
implementation of the viscoelastic contribution at boundaries, namely inlet,
outlet and junction, is carried out considering the natively hyperbolic nature
of the mathematical model. A non-linear system is established based on the
definition of a specific Riemann Problem at junctions. Basic junction test
cases are analyzed, such as a trivial 2-vessels junction, for both a generic
artery and a generic vein, and a simple 3-vessels junction, considering an
aortic bifurcation scenario. The chosen IMEX-SSP2 Finite Volume scheme is
demonstrated to be second-order accurate in the whole domain and well-balanced,
even when including junctions. Two benchmark arterial networks are then
implemented, differing in number of vessels and in viscoelastic parameters.
Comparison of the results obtained in the two networks underlines the high
sensitivity of the model to the chosen viscoelastic parameters. In these
numerical tests, the network-wise conservation of the contribution provided by
the viscoelastic characterization of vessel walls is assessed. Finally,
numerical results of pressure waveforms in the common carotid artery and in the
femoral artery of the second network simulated are compared with in-vivo data.
