The aim of the present work is to address the closure problem for hemodynamic simulations by developing a flexible and effective model that accurately distributes flow in the downstream vasculature and can stably provide a physiological pressure outflow boundary condition. To achieve this goal, we model blood flow in the sub-pixel vasculature by using a non-linear 1D model in self-similar networks of compliant arteries that mimic the structure and hierarchy of vessels in the meso-vascular regime (radii [Formula: see text]). We introduce a variable vessel length-to-radius ratio for small arteries and arterioles, while also addressing non-Newtonian blood rheology and arterial wall viscoelasticity effects in small arteries and arterioles. This methodology aims to overcome substantial cut-off radius sensitivities, typically arising in structured tree and linearized impedance models. The proposed model is not sensitive to outflow boundary conditions applied at the end points of the fractal network, and thus does not require calibration of resistance/capacitance parameters typically required for outflow conditions. The proposed model convergences to a periodic state in two cardiac cycles even when started from zero-flow initial conditions. The resulting fractal-trees typically consist of thousands to millions of arteries, posing the need for efficient parallel algorithms. To this end, we have scaled up a Discontinuous Galerkin solver that utilizes the MPI/OpenMP hybrid programming paradigm to thousands of computer cores, and can simulate blood flow in networks of millions of arterial segments at the rate of one cycle per 5 min. The proposed model has been extensively tested on a large and complex cranial network with 50 parent, patient-specific arteries and 21 outlets to which fractal trees where attached, resulting to a network of up to 4,392,484 vessels in total, and a detailed network of the arm with 276 parent arteries and 103 outlets (a total of 702,188 vessels after attaching the fractal trees), returning physiological flow and pressure wave predictions without requiring any parameter estimation or calibration procedures. We present a novel methodology to overcome substantial cut-off radius sensitivities.