The high-resolution 3D groundwater flow and transport simulation problem requires massive discrete linear systems to be solved, leading to significant computational time and memory requirements. The domain decomposition method is a promising technique that facilitates the parallelization of problems with minimal communication overhead by dividing the computation domain into multiple subdomains. However, directly utilizing a domain decomposition scheme to solve massive linear systems becomes impractical due to the bottleneck in algebraic operations required to coordinate the results of subdomains. In this paper, we propose a two-level domain decomposition method, named dual-domain decomposition, to efficiently solve the massive discrete linear systems in high-resolution 3D groundwater simulations. The first level of domain decomposition partitions the linear system problem into independent linear sub-problems across multiple subdomains, enabling parallel solutions with significantly reduced complexity. The second level introduces a domain decomposition preconditioner to solve the linear system, known as the Schur system, used to coordinate results from subdomains across their boundaries. This additional level of decomposition parallelizes the preconditioning of the Schur system, addressing the bottleneck of the Schur system solution while improving its convergence rates. The dual-domain decomposition method facilitates the partition and distribution of the computation to be solved into independent finely grained computational subdomains, substantially reducing both computational and memory complexity. We demonstrate the scalability of our proposed method through its application to a high-resolution 3D simulation of chromium contaminant transport in groundwater. Our results indicate that our method outperforms both the vanilla domain decomposition method and the algebraic multigrid preconditioned method in terms of runtime, achieving up to 8.617× and 5.515× speedups, respectively, in solving massive problems with approximately 108 million degrees of freedom. Therefore, we recommend its effectiveness and reliability for high-resolution 3D simulations of groundwater flow and transport.