Nitrification-denitrification processes in the hyporheic zone control the dynamics of dissolved inorganic nitrogen (DIN) species and can lead to production of nitrous oxide, which contributes to the greenhouse effect. We consider DIN dynamics in an advection-dominated regime, wherein transport and reactions occur along streamlines crossing hyporheic sediments. Our focus is on the impact of uncertainty in both stream-water quality and rate constants of the subsurface reactions on predictions of DIN concentrations. We derive equations for a joint probability density function (PDF), and corresponding marginal PDFs and cumulative distribution functions (CDFs), of the species concentrations. Their derivation requires a novel closure, which depends on the mean and (co)variance of the species concentrations. We use streamline coordinates to reduce the dimensionality of the PDF/CDF equations and, hence, the computational effort of solving them. For the sake of completeness, we also present similar equations in Cartesian coordinates. By providing a complete probabilistic description of species concentrations at the bedform scale, our PDF/CDF equations allow one to evaluate the impact of random spatiotemporal variability in the inputs on the DIN dynamics. They yield physically-based prior distributions for data assimilation and can be deployed to guide measurement campaigns by identifying regions with largest predictive uncertainty.