We consider a linearized dynamical system modeling the flow rate of water along the rivers and hillslopes of an arbitrary watershed. The system is perturbed by a random rainfall in the form of a compound Poisson process. The model describes the evolution, at daily time scales, of an interconnected network of linear reservoirs and takes into account the differences in flow celerity between hillslopes and streams as well as their spatial variation. The resulting stochastic process is a piece-wise deterministic Markov process of the Orstein–Uhlembeck type. We provide an explicit formula for the Laplace transform of the invariant density of streamflow in terms of the geophysical parameters of the river network and the statistical properties of the precipitation field. As an application, we include novel formulas for the invariant moments of the streamflow at the watershed’s outlet, as well as the asymptotic behavior of extreme discharge events, and conditions for the statistical scaling of streamflow with respect to Horton order.