The spatial distribution of physical and chemical heterogeneities greatly impacts solute transport in groundwater and is critical in many subsurface applications. In this work we perform flow-through experiments in multidimensional setups and we combine them with forward and inverse reactive transport modeling to explore the capability of imaging pyrite inclusions in the subsurface. We studied the oxidative dissolution of pyrite in multidimensional setups, including 1-D columns and a 2-D flow-through chamber. Pyrite dissolution was recently studied in batch setups where the consumption of oxygen, the change in pH and the release of dissolved iron and sulfur were measured to constrain the dissolution kinetics of the reactive minerals. In this work, we embedded pyrite inclusions in a sandy porous matrix. The reactive mineral and the inert sandy matrix had the same grain size in order to obtain physically homogeneous but chemically heterogeneous media. Spatially distributed optode measurements of dissolved oxygen allowed tracking the propagation of an oxic groundwater front in the initially anoxic 1-D columns and 2-D flow-through chamber. The non-invasive optode sensor measurements of O2 were performed at high-resolution (i.e., spacing of 2.5 mm) along the axis of the 1-D column and at different cross sections in the 2-D flow-through setup. These data, as well as the analysis of dissolved species originating from pyrite oxidation and measured at the outlet of the flow-through systems, were compared with the outcomes of a reactive transport model describing the main physical and geochemical processes in the setups. Finally, the Principal Component Geostatistical Approach (PCGA) was applied to image the spatial distribution of the pyrite inclusions. The results show the capability of the proposed methodology to accurately identify both the spatial locations and the concentration of the reactive mineral zones embedded in the 1-D and 2-D systems.