Magnetic resonance imaging (MRI) was used to visualize the NAPL source zone architecture before and after surfactant-enhanced NAPL dissolution in three-dimensional (3D) heterogeneously packed flowcells characterized by different longitudinal correlation lengths: 2.1 cm (aquifer 1) and 1.1 cm (aquifer 2). Surfactant flowpaths were determined by imaging the breakthrough of a paramagnetic tracer (MnCl(2)) analyzed by the method of moments. In both experimental aquifers, preferential flow occurred in high permeability materials with low NAPL saturations, and NAPL was preferentially removed from the top of the aquifers with low saturation. Alternate flushing with water and two surfactant pulses (5-6 pore volumes each) resulted in approximately 63% of NAPL mass removal from both aquifers. However, overall reduction in mass flux (Mass Flux 1) exiting the flowcell was lower in aquifer 2 (68%) than in aquifer 1 (81%), and local effluent concentrations were found to increase by as high as 120 times at local sampling ports from aquifer 2 after surfactant flushing. 3D MRI images of NAPL revealed that NAPL migrated downward and created additional NAPL source zones in previously uncontaminated areas at the bottom of the aquifers. The additional NAPL source zones were created in the direction transverse to flow in aquifer 2, which explains the higher mass flux relative to aquifer 1. Analysis using a total trapping number indicates that mobilization of NAPL trapped in the two coarsest sand fractions is possible when saturation is below 0.5 and 0.4, respectively. Results from this study highlight the potential impacts of porous media heterogeneity and NAPL source zone architecture on advanced in-situ flushing technologies.