The overall objective of this thesis is to better understand the nature, magnitude, and spatial variability of groundwater flow through a complex carbonate aquifer and identify how this variability relates to the local flow system. To accomplish this, innovative field methods were applied at a research site in Guelph, Ontario with nine closely-spaced boreholes spanning the full thickness of an important local aquifer. Field measurements using fibre optic Active Distributed Temperature Sensing (A-DTS) in sealed boreholes identified the location and magnitude of natural gradient flow within the aquifer sequence, and to what extent flow is affected by hydraulic cross-connection and pumping. Results show quantifiable changes in head and flow rates at discrete fractures to both hydraulic conditions. Variation in flow direction at different depths was measured using a new multi-borehole A-DTS method that determines natural gradient flow direction by measuring the variation of flow rate around a pumping well. Results show a large variation in flow direction throughout the sequence and are generally consistent with horizontal hydraulic gradient directions at similar depths. Flow in the shallow rock appears to be towards the Eramosa River, 400 m to the north. Slightly deeper, flow swings to the south possibly due to flow towards a rubble layer in a buried bedrock valley, which is at a similar elevation 350 m to the south. The strongest evidence of active flow corresponds to fractures, though there is evidence of flow within a 2 m cavern in the Gasport Formation. Four distinct hydrogeological units were identified on the site and the interfaces between them correspond to laterally extensive bedding parallel fractures. Aquitard properties were observed in the Ancaster Member and not in the Vinemount Member, the latter of which is traditionally associated with being a regional aquitard. This study presents detailed field examples demonstrating the variability in flow magnitude and direction, as well as the possible influence of surface water and bedrock valleys. These details can be used to produce more robust hydrogeologic models that better represent flow paths from recharge to discharge zones and inform wellhead protection and contaminant transport and fate.