We adapt a Bayesian Hierarchical Network Model (BHNM) for modeling and ensemble forecasting of daily river stages at multiple gauges on a rainfed river network. The stage at a gauge on any day is modeled as a Probability Density Function (PDF), with parameters varying temporally and spatially across the gauges. The PDF parameters vary as a function of covariates that include stage, streamflow, and precipitation from previous times at upstream gauges and catchment areas between gauges. This leverages the network structure of the river to capture spatial correlation along with the river basin hydrologic processes encapsulated in them. With suitable priors on the model parameters, likelihood functions, and, using Markov Chain Monte Carlo simulation approach predictive posterior distributions of all the space-time model parameters are obtained and, consequently, that of river stages across the network for any day. The best model, which includes the PDF type and sub-set of covariates, is obtained via objective criteria, Deviance Information Criteria (DIC). The model is demonstrated by applying it to daily stages at four gauges on the monsoon rainfed Narmada River in western India during the peak monsoon period (July-August) of 1978 to 2014 period. Model validation in cross-validation mode shows skillful and reliable forecasts of river stages, including higher stages that correspond to floods relative to a null-model of linear regression. Since river stages are often used in flood disaster management and preparedness efforts, this model and the skillful outputs enhance the potential for effective real time flood warning and mitigation.