This study extends the Algebraic Dynamic Multilevel (ADM) method for simulating contaminant transport and retention in vadose zones. Building upon a fully implicit scheme that couples variably saturated flow and contaminant transport, the developed ADM framework effectively predicts contaminant plume migration across both unsaturated and saturated media under heterogeneous conditions. During the simulation, ADM dynamically adjusts grid resolution based on the spatial gradients of primary variables, applying fine-scale grids in regions with steep gradients and coarsening the mesh where fields remain smooth. These dynamic adjustments are achieved through prolongation and restriction operators that transfer solutions across multilevel grid systems. As both water content and contaminant concentration evolve spatiotemporally, dual coarsening criteria are introduced to simultaneously capture flow and transport dynamics. Results show that the developed model reproduces the contaminant migration obtained from the fully resolved solution using substantially fewer grids. Moreover, it offers the flexibility to trade off numerical accuracy against computational cost by selecting an appropriate coarsening criterion.

Underground hydrogen (H₂) storage in saline aquifers is a viable solution for large-scale H₂ storage. Due to its remarkably low viscosity and density, the flow of H₂ within saline aquifers exhibits strong instability, which needs to be thoroughly investigated to ensure safe operations at the storage site. For the first time, we develop a lattice Boltzmann model tailored for pore-scale simulations of the H₂-brine system under typical subsurface storage conditions. The model captures the significant contrast of fluid properties between H₂ and brine, and it offers the flexibility to adjust the contact angle to suit varying wetting conditions. We show that the snap-off is enhanced in a system with a high capillary number and a small contact angle. These conditions lead to a low recovery factor, which is unfavorable for H₂ production from the aquifer. Moreover, the relative permeability curves, computed from the simulation results, exhibit distinct behaviors for H₂ and brine. In the case of the wetting phase, the relative permeability can be quantified using the quadratic expression, whereas for the non-wetting phase, the relative permeability exhibits a nearly linear behavior, and saturation alone appears insufficient to characterize the relative permeability at large saturations of non-wetting phase. This implies that different formula for liquid and gas phases may be employed for continuum-scale simulations.