Concentration-discharge relationship is crucial in river hydrology, as it reflects water quality dynamics across both low- and high-flow regimes. However, its mathematical description is still challenging owing to the underlying complex physics and chemistry. This study proposes an infinite-dimensional stochastic differential equation model that effectively describes the concentration-discharge relationship while staying analytically tractable, along with the computational aspects of the model. The proposed model is based on the superposition of the square-root processes (or Cox-Ingersoll-Ross processes) and its variants, through which both the long-term moments and autocovariance of river discharge and the fluctuation of water quality index can be derived in closed forms. Particularly, the model captures both long (power decay) and short (exponential decay) memories of the fluctuation in a unified manner, while quantifying the hysteresis in the concentration-discharge relationship through mutual covariances with time lags. Based on a verified numerical method, the model is computationally applied to weekly data on total nitrogen (TN. long memory with moderate fluctuation), total phosphorus (TP. short memory with large fluctuation), and total organic carbon (TOC. short memory with moderate fluctuation) from a rural catchment to validate its applicability to real-world datasets. Based on the identified model and its mutual covariance, our findings indicate that, on average, the peak concentrations of these water quality indices appear approximately 1 day after discharge. Finally, the study discusses the effects of model uncertainty on mutual covariance.