Discontinuous Galerkin (DG) method is presented for numerical modeling of melt migration in a chemically reactive and viscously deforming upwelling mantle column. DG methods for both advection and elliptic equations provide a robust and efficient solution to the problems of melt migration in the asthenospheric upper mantle. Assembling and solving the elliptic equation is the major bottleneck in these computations. To address this issue, adaptive mesh refinement and local timestepping methods have been proposed to significantly improve the computational wall time. The robustness of DG methods is demonstrated through two benchmark problems by modeling detailed structure of high-porosity dissolution channels and compaction-dissolution waves.