<p>This paper presents a comprehensive numerical framework for simulating convection-dominated transport phenomena in heterogeneous porous media, combining a weak Galerkin finite element method for spatial discretization on general polygonal meshes with a high-order diagonally implicit Runge–Kutta scheme for temporal integration. The method effectively addresses key challenges in porous media simulations, including handling of sharp gradients in convection-dominated flows, accurate representation of complex geometries and material interfaces, and robust treatment of nonlinear reaction terms. Theoretical analysis establishes unconditional stability and optimal convergence rates of <InlineEquation ID="IEq1"> <EquationSource Format="TEX">\(\mathcal {O}(\tau ^p + h^k)\)</EquationSource> <EquationSource Format="MATHML"><math> <mrow> <mi mathvariant="script">O</mi> <mo stretchy="false">(</mo> <msup> <mi>τ</mi> <mi>p</mi> </msup> <mo>+</mo> <msup> <mi>h</mi> <mi>k</mi> </msup> <mo stretchy="false">)</mo> </mrow> </math></EquationSource> </InlineEquation>, where <i>p</i> and <i>k</i> represent the temporal and spatial orders of accuracy, respectively. The analysis is extended to mixed boundary conditions commonly encountered in practical applications. Numerical experiments demonstrate the method’s superiority over conventional approaches, particularly in maintaining solution quality for high Péclet number flows and complex heterogeneous domains. The scheme’s computational efficiency is validated through CPU-time comparisons against established methods for equivalent accuracy targets. Additional investigations confirm the method’s robustness for strongly nonlinear reactions and its extensibility to three-dimensional problems. The conservation properties and ability to handle general meshes make it particularly suitable for practical porous media applications including contaminant transport and enhanced oil recovery, as demonstrated through a realistic hydrogeological simulation.</p>