We introduce a comprehensive methodology that combines multi-objective optimization, global sensitivity analysis (GSA) and 3D groundwater modeling to analyze subsurface flow dynamics across large-scale domains. In this way, we effectively consider the inherent uncertainty associated with subsurface system characterizations and their interactions with surface waterbodies. We demonstrate the effectiveness of our proposed approach by applying it to the largest groundwater system in Italy.