This paper concerns about a stepwise modelling procedure for groundwater flow simulation in a complex carbonate, folded and faulted, multilayer aquifer, in the Apennine chain in Central Ita-ly, which constitutes a source of good quality water for human consumption. A perennial river acts as the main natural drain for groundwater while sustaining valuable water related ecosys-tems. The spatial distribution of recharge was estimated using the Thornthwaite-Mather method on 60 years of climate data. The system was conceptualized as three main aquifers separated by two locally discontinuous aquitards. Three numerical models were implemented by gradually adding complexity to the model grid: single layer (2D), three layers (Quasi-3D), and five layers (Fully-3D), using an equivalent porous medium approach, in order to find the best solution with a parsimonious model setting. To overcome dry-cell problems in the Fully-3D model, the New-ton-Raphson formulation for MODFLOW-2005 was invoked. Calibration results show that a Fully-3D model was required to match the observed distribution of aquifer outflow to the river baseflow. The numerical model demonstrated the major impact of folded and faulted geological structures on controlling the flow dynamics in terms of flow direction, water heads and spatial distribution of the outflows to the river and springs.