Flow heterogeneity controls dissolution dynamics in topologically complex rocks

Rock dissolution is a common subsurface geochemical reaction affecting pore space properties, crucial for reservoir stimulation, carbon storage, and geothermal energy. Predictive models for dissolution remain limited due to incomplete understanding of the mechanisms involved. We examine the influence of flow, transport, and reaction regimes on mineral dissolution using time-resolved data from 3D rocks. We find that initial pore structure significantly influences the dissolution pattern, with reaction rates up to two orders of magnitude lower than batch conditions given solute and fluid-solid boundary constraints. Flow unevenness determines the location and rate of dissolution. We propose two models describing expected dissolution patterns and effective reaction rates based on dimensionless metrics for flow, transport, and reaction. Finally, we analyze feedback between evolving flow and pore structure to understand conditions that regulate/reinforce dissolution hotspots. Our findings underscore the major impact of flow arrangement on reaction front propagation and provide a foundation for controlling dissolution hotspots.

Methods
This data set contains (a) the initial steady-state Eulerian flow field at the pore-scale and (b) the temporal evolution of the dissolution reaction rate for various rock samples considered in this study. Briefly, for resolving the flow field, binary 3D rock images at time t=0 are used as the geometric boundaries in which the Navier-Stokes equations are solved numerically with the finite volume method (FlowDict, Math2Market). The boundary conditions imposed a uniform pressure gradient of 1kPa at opposite boundaries of each sample volume. For tracking the evolution of the reaction rate, reactive transport simulations are performed with the Digital Rock Physics module of GeoDict (Math2Market), starting with the initial rock image. Briefly, each dissolution simulation is run in so-called batches. During each batch, (i) the pore geometry is updated, (ii) the flow within that geometry is resolved to achieve an average linear velocity corresponding to the desired Peclet number of the study (iii) N=50,000 acidified particles are tracked in the domain via a flux-weighted injection, and (iv) collisions between particles and the solid surface are recorded. At the end of a batch, the solid voxels where collisions occurred are removed to update the geometry for the next batch. Only one voxel layer is permitted to be dissolved in each batch.  Steps i-iv are repeated for each batch until the sample meets quasi-equilibrium for its reaction rate (i.e., dre/dt < 5%).