Reactive Transport Model of Gypsum Karstification in Physically and Chemically Heterogeneous Fractured Media

Gypsum dissolution leads to the development of karstic features within much shorter timescales than in other sedimentary rocks, potentially leading to rapid deterioration of groundwater quality and increasing the risk of catastrophes caused by subsidence. Here, we present a 2-D reactive transport model to evaluate gypsum karstification in physically and chemically heterogeneous systems. The model considers a low-permeability rock matrix composed mainly of gypsum and a discontinuity (fracture), which acts as a preferential water pathway. Several scenarios are analyzed and simulated to investigate the relevance for gypsum karstification of: (1) the dynamic update of flow and transport parameters due to porosity changes; (2) the spatial distribution of minerals in the rock matrix; (3) the time evolution of water inflows through the boundaries of the model; (4) the functions relating permeability, k, to porosity, phi. The average porosity of the matrix after 1000 years of simulation increases from 0.045 to 0.29 when flow, transport, and chemical parameters and the water inflows through the boundary are dynamically updated according to the porosity changes. On the contrary, the porosity of the matrix hardly changes when the porosity feedback effect is not considered, while its average increases to 0.13 if the water inflow occurs through the discontinuity. Moreover, the dissolution of small amounts of highly soluble sulfate minerals plays a major role in the development of additional fractures. The increase in hydraulic conductivity is largest for the power law with an exponent of n = 5, as well as the Kozeny-Carman and the modified Fair-atch k-phi relationships. The gypsum dissolution front propagates into the matrix faster when the power law with n = 2 and 3 and the Verma-Pruess k-phi relationships are used.