Mixing-limited reactions are central to a wide range of processes in natural and engineered porous media. Recent advances have shown that concentration gradients sustained by flow at the pore-scale influence macroscopic reaction rates over a large range of reactive transport regimes. Yet, resolving concentration gradients driven by fluid mixing at the pore-scale is challenging with current simulation methods. Here, we introduce a computational methodology to resolve concentration gradients at the pore scale in mixing-limited reactions. We consider a steady-state reactive transport problem characterized by reactive fluids flowing in parallel in a porous material. Given a mesh representation of the pore space and a steady velocity field, we solve the steady advection-diffusion equation for conservative scalar transport using a stabilized finite-element method combined with mesh refinement adapted to local scalar gradients. Based on this solution and assuming instantaneous reaction kinetics in the fluid, we infer the distribution of species involved in an irreversible bi-molecular reaction. We validate the method by comparing our results for uniform flow with analytical solutions and then apply it to simulate mixing-limited reactions in a three-dimensional random bead pack and Berea sandstone sample. Chaotic flow within the pore space leads to sustained concentration gradients, which are captured by our numerical framework. The results underscore the ability of the methodology to simulate transverse mixing and mixing-limited reactions in complex porous media and to provide bottom-up numerical data to improve the prediction of effective reaction rates at larger scales.