Pyrite oxidation (FeS2) causes acidification and mobilization of metals. Mathematical modeling of pyrite oxidation in variably saturated reactive flow systems is challenging because (1) it occurs through a complex interplay of multi-phase flow and transport processes, and (2) aqueous concentrations of key species vary over tens orders of magnitude in different redox conditions. Here, we present a general multi-phase reactive transport model for redox processes. Two alternative implementations were made in the TOUGHREACT and TOUGH2-CHEM simulation codes which use sequential iteration and simultaneous solution, respectively. Both codes are used to simulate a fully and a variably saturated pyrite oxidation problem with simple 1-D flow and reaction conditions. Results from both codes indicate that the effects of oxygen partial pressure reduction due to reactions on the fluid flow is not significant under ambient conditions. However, it must be noted that when fluid flow and chemical reactions are strongly coupled, such as when boiling takes place in geothermal reservoirs, this could be essential. The fully simultaneous approach has a complete process description. The sequential iteration approach is found to be more efficient computationally. The oxygen gas diffusion process plays a dominant role in the chemical evolution for pyrite oxidation in unsaturated conditions. An example in 2-D fractured rock is presented to demonstrate pyrite oxidation under complex flow and geochemical conditions. This example shows that pyrite oxidation exerts strong influence on hydrogeochemical evolution in variably saturated flow systems. The alteration of primary rock minerals and the development of secondary mineral assemblages predicted are consistent with field observations. This example serves as a prototype for oxidative weathering processes with broad significance for geoscientific, engineering, and environmental applications.