Agricultural nitrate pollutants infiltrate into the subsurface and contaminate groundwater. The redox environment in the subsurface is important for the natural removal of nitrate by denitrification. Detailed knowledge of the redox conditions is needed in order to make better-targeted nitrogen regulations for farmers. However, unveiling three-dimensional (3D) redox architectures is challenging because one only observes redox conditions in boreholes. Therefore, this work proposes a combination of towed transient electromagnetic resistivity (tTEM) geophysical surveys and redox boreholes to model 3D redox architecture stochastically. The tTEM survey reveals the geological structure in high resolution. However, the tTEM survey and redox boreholes are often non-colocated. To address this issue, geostatistical simulations are performed to generate multiple resistivity data colocated with redox boreholes. Then, a statistical learning method, namely multinomial logistic regression, is leveraged to predict multiple 3D redox architectures given the uncertain surrounding resistivity structures. In this way, the aggregated prediction of multiple redox architectures has a higher prediction accuracy than a redox prediction model with interpolated resistivity. The trained statistical model can also identify significant resistivity structures for redox predictions. An inverse problem has also been formulated to better match the redox borehole data. In summary, the proposed workflow models 3D resistivity and redox architecture jointly, aggregates to a highly accurate redox architecture, and provides important resistivity structures for domain experts. The highly accurate redox architecture supports a better agricultural regulation decision.