We combined electrical resistivity tomography (ERT) on land and in a stream with zone-based hydraulic conductivities (from multi-level slug testing) to investigate the local geological heterogeneity of the deposits in a wetland–stream system. The detailed geology was incorporated into a numerical steady-state groundwater model that was calibrated against average head observations. The model results were tested against groundwater fluxes determined from streambed temperature measurements. Discharge varied up to one order of magnitude across the stream and the model was successful in capturing this variability. Water quality analyses from multi-level sampling underneath the streambed and in the wetland showed a stratification in groundwater composition with an aerobic shallow zone with oxygen and nitrate (top ∼3 m) overlying a reduced, anoxic zone. While NO3- concentrations up to 58 mg L−1 were found in the top of the aquifer and immediately underneath the streambed no NO3- was detected deeper within the aquifer. An inverse relationship between NO3- and SO42- suggests that pyrite oxidation takes place in the deeper parts of the aquifer. Simulated flow path lines showed very different origins for deeper groundwater samples. No nitrate reduction is believed to occur in the shallow zone, where oxygen is present, and the residence time is on the order of 1 year. Nitrate reduction can, however, occur in the deeper parts, which are oxygen-free, and where the residence time is on the order of 7 years. A simulation with a homogeneous model did not match the observations nearly as well as a heterogeneous model based on ERT and a spatially distributed hydraulic conductivity. Furthermore, the origin of the sampled groundwater could not have been predicted from groundwater hydraulic head and the groundwater chemistry alone. The presented approach of integrating such methods in groundwater–surface water exchange studies, proved efficient to obtain information of the controlling factors.