Semi-analytical method for matrix diffusion in heterogeneous and fractured systems with parent-daughter reactions

A semi-analytical/numerical method for modeling matrix diffusion in heterogeneous and fractured groundwater systems is developed. This is a significant extension of the Falta and Wang (2017) method that only applied to diffusion in an aquitard of infinite thickness. The current solution allows for the low permeability matrix to be embedded within a numerical gridblock, having finite average thickness, a specified volume fraction and a specified interfacial area with the high permeability domain. The new formulation also allows for coupled parent-daughter decay reactions with multiple species that each have independent retardation factors, decay rates, and yield coefficients in both the high and low permeability parts of the system. The method uses a fitting function to approximate the transient concentration profile in the low permeability part of each gridblock so that the matrix diffusion flux into the high permeability part of the gridblock can be computed as a concentration dependent source-sink term. This approach is efficient because the only unknowns in each gridblock are the concentrations in the high permeability domain, so there is practically no increase in computational effort compared to a conventional transport simulation. The method is shown to compare favorably with an analytical solution for matrix diffusion in fractured media with parallel fractures, with an analytical solution for matrix diffusion with parent-daughter decay reactions, with laboratory experiments of matrix diffusion in a layered system, with a laboratory experiment involving lens shaped inclusions, and with fine grid numerical simulations of transport in highly heterogeneous systems.