An HDG-type discretization of 1D embedded blood flow in complex networks with Riemann-coupled junctions

We present a hybridizable discontinuous Galerkin (HDG) method for one-dimensional blood flow in arterial networks with Windkessel boundary conditions. The scheme uses a Harten--Lax--van Leer (HLL) numerical flux and implicit time integration. After spatial discretization, the complete network---including cell states, hybrid traces, junction conditions, and dynamic RCR terminal pressures---is written as a residual index-one differential--algebraic equation (DAE) and integrated with variable-step, variable-order backward differentiation formulas. The implicit formulation removes the CFL stability restriction of explicit schemes and allows variable-order, variable-step solvers. The resulting single global residual and Jacobian for the network and its terminal dynamics also provide a natural interface for monolithic coupling with external nonlinear models. We establish local solvability and a semi-discrete energy balance, and show numerically optimal convergence rates for both solution variables in the L² norm and the broken H¹ seminorm. The solver is validated against established benchmarks, including an aortic bifurcation and a $37$-artery network, and is further applied to a cerebral arterial network. Mass conservation is observed in a synthetically generated anatomical liver network of 8\,k vessels, and strong-scaling studies are performed on synthetic networks of 50\,k and 100\,k vessels.