Redox cycles, geochemistry, and pH are recognized as key drivers of subsurface biogeochemical cycling in wetland ecosystems but are typically not included in land surface models. These omissions may introduce errors when simulating carbon cycling and greenhouse gas emissions in systems where redox interactions and pH fluctuations are important, such as wetland-rich tundra landscapes where saturated soil conditions combined with high soil iron concentrations can influence biogeochemistry, or coastal regions where sulfate concentrations associated with saltwater influence can drive biogeochemical contrasts across salinity gradients. Here, we coupled the Energy Exascale Earth System Model (E3SM) Land Model (ELM) with geochemical reaction network simulator PFLOTRAN, allowing geochemical processes and redox interactions to be integrated with land surface model simulations. We implemented a reaction network including aerobic decomposition, fermentation, iron oxide reductive dissolution and dissolved iron oxidation, sulfate reduction, sulfide oxidation, methanogenesis, methane oxidation, and pH dynamics. We used the model framework to simulate biogeochemical cycling and methane production across redox gradients in arctic tundra permafrost polygons and thermokarst features, and in coastal wetlands across gradients of salinity. Model simulations were parameterized using laboratory incubations and literature values and evaluated using measured porewater concentrations and surface gas emissions from wetland field sites across coastal and arctic regions of the United States. These results demonstrate how directly simulating biogeochemical reaction networks can improve land surface model simulations of subsurface biogeochemistry and carbon cycling in wetland ecosystems, and highlight the value of porewater biogeochemical data for evaluating process-based wetland biogeochemical models.