Discharge in many rivers is often fed by outflow from a shallow groundwater reservoir. It is becoming clear that the outflow from this aquifer is not linearly proportional to storage as is commonly assumed in many algorithms. Numerical analysis of flow recession curves from about 100 river gauging stations instead reveals a nonlinear relationship between baseflow, Q, and storage, S, for which the equation S = aQ^b was adopted. Values of the exponent b are found by calibration to be between 0 and 1 but with a high concentration around 0.5, which is in accordance with the findings of other studies and theoretical approaches yielding b = 0.5 for unconfined aquifers and relating the coefficient a to catchment properties, primarily area and shape of basin, pore volume and transmissivity. This non-linear reservoir function is proposed as a more realistic alternative to the linear reservoir function. The relatively fast response of groundwater flow to rainfall is mainly a result of the increase of hydraulic head of the groundwater reservoir accelerating the exfiltration of 'old', pre-event water into the river bed. As fissure and pore volumes communicate hydraulically, it appears physically reasonable to model the system by one non-linear reservoir for catchments, or parts of them, instead of applying independent parallel linear reservoirs. The non-linear reservoir algorithms are supported by an analytical derivation. They are extended for the automatic separation of baseflow from a time-series of daily discharge in rivers and the computation of storage and effective recharge of groundwater in river basins by inverse nonlinear reservoir routing. The time-series obtained allow the identification and quantification of long-term changes to the water balance. Relationships between computed groundwater storage and observed groundwater level can also be established.