Upscaling Strategies of Porosity-Permeability Correlations in Reacting Environments from Pore-Scale Simulations

In geochemically reacting environments, the mineral dissolution and precipitation alters the structural and transport properties of the media of interest. The chemical and structural heterogeneities of the porous media affect the temporal evolution of the permeability with respect to porosity. Such correlations follow a nonlinear trend, which is difficult to estimate a priori and without knowledge of the microstructure itself, especially under the presence of strong chemical gradients. Macroscopic fieldscale codes require such an input, and in the absence of exact descriptions, simplified correlations are used. After highlighting the diversity of microstructural evolution paths, due to dissolution, we discuss possible upscaling strategies.


Introduction
Precipitation and dissolution reactions in porous media dominate and control a large number of geochemical processes and industrial applications. The precipitation and dissolution of minerals from aqueous solutions alters the pore space and its connectivity. This has a strong effect on the mass convection and diffusion through the porous medium. When a mineral precipitates/dissolves at the reactive porous surface, the overall porosity decreases/increases, leading to a subsequent decrease/increase in permeability and effective diffusivity. At the same time, the connectivity of the pores can also change in a way to block or to facilitate the mass diffusion processes.
Reactive transport modelling at the field scale is usually based on a macroscopic finite element or a finite volume discretization scheme [1]. In such descriptions, the computational domain is divided into small elements/volumes, the so-called voxels. The voxels are typically several orders of magnitude larger than the typical pore diameter, and as a consequence, all chemical and transport properties within such volumes are homogenized and smoothed out. The pore space and its transport properties are therefore represented using macroscopic parameters, such as the porosity, the tortuosity, the diffusivity, and the permeability. In such a description, the small scale geometrical characteristics and the heterogeneities of the porous materials are neglected. Such an assumption allows making accurate numerical predictions in the case of relatively mild chemical gradients, as well as in the case where chemical reactions can be approximated by equilibrium values. However, when strong chemical gradients are present with simultaneous dissolution and precipitation reactions, the pure macroscopic reactive transport codes fail to make accurate predictions of the evolution of the system. Dissolution and precipitation reactions change the pore space and the resulting material properties, in a nonlinear way, and therefore have a strong feedback in the transport properties of the porous medium. The lack of explicit pore structure description and of appropriate material-specific correlations is responsible, for example, for the numerical artefact, relevant to the dependency of the resulting clogging times on the spatial grid discretization [1][2][3][4][5]. An improvement to the numerical predictions can be achieved (a) by coupling pore-level solvers to the macroscopic ones in a multiscale simulator [6][7][8] or (b) by providing the necessary microscopic feedback in terms of appropriate correlations or tabulated values, which can be defined from pore-level simulations (upscaling of results). Pore-level methods allow the simulation of the advection-diffusion-reaction processes at the pore space, where surface charges and reactive surface areas can be explicitly considered. Representative structures can be generated via computer models or can be obtained via X-ray or other microtomography techniques. When combined with appropriate kinetic and thermodynamic solvers that act at the submicrometer level, it is possible to reproduce accurately the experimental observations.
Depending on the level of abstraction, different porelevel methodologies exist. The more detailed ones solve the relevant flow equations or some approximation depending on the flow regime and the flow physics that are involved, in realistic geometries. Such methods are the lattice Boltzmann method [9][10][11][12][13], the methods based on particle hydrodynamics [14,15], as well as the standard finite volume methods when applied in complex geometries with moving boundaries [16][17][18]. Lattice Boltzmann models can resolve transport processes in realistic complex geometries, involving complex interactions between species and phases, but are more computational intensive compared to pore network models. A significant advantage of the lattice Boltzmann methodology is the minimum effort to discretize the realistic computational domain, as well as the efficient continuous solid structure update per time step. Such an example is the evolution of a system when simultaneous dissolution and precipitation processes are competing and drive the evolution of the system [11,19]. Efficient parallelization though allows running simulations with many billion degrees of freedom in relatively small computer clusters, especially when GPGPUs are used [20,21]. We note that the numerical extraction of microscopic properties using realistic geometries has the potential to provide useful input to the effective medium theories, which are used to upscale porous medium flows [22] In this paper, we construct pore geometries with target porosity and initial permeability following a methodology similar to [23,24]. The target permeability is selected in a way to represent limestone rock samples commonly found in hydrocarbon reservoirs or in geothermal fields. The implemented chemical reaction is representing the calcite dissolution under the presence of acid, a common process used in the field stimulation, in order to enhance the permeability of the formation. The evolution of selected geometries is examined using the lattice Boltzmann framework, and permeability-porosity correlations are numerically extracted. Upscaling strategies that allow passing information to the macroscopic solvers and therefore bridge pore level and macroscopic scales, are discussed based on the output of the simulations.

Reactive Transport Modelling
2.1. Pore-Level Modelling. For the simulation presented in the next sections, the lattice Boltzmann methodology is implemented. This method is a special discretization of the Boltzmann equation. The elementary variables are the so-called populations or velocity probability distribution functions f i . At every distribution function corresponds to a discrete velocity vector [25][26][27]. Different discretization schemes lead to different numbers of discrete velocities, which results in several lattice models [28,29]. For twodimensional simulations, the standard D2Q9 square lattice with 9 discrete velocities is usually implemented due to its simplicity and robustness in complex geometry domains (see Figure 1).
For the modelling of the advection-diffusion and precipitation processes, a multicomponent LB model is used as described in [19]. For the sake of completeness, we briefly present it also here. The model is composed of a basis fluid medium that recovers the Navier-Stokes equations at the macroscopic limit, plus a passive scalar-coupled population set that simulates the diffusion of ions. The isothermalguided equilibrium nine-velocity model (D2Q9 lattice) of Prasianakis et al. [30] is selected as the basis model. The discrete velocities of populations f i for i = 0-8 are c i = 0, 0 for i = 0, c i = ±1, 0 and (0, ±1) for i = 1-4, and c i = ±1, ±1 for i = 5-8 [29].
The following population moments correspond to the density of the solution ρ and the momentum j a in the direction a = x, y: The guided equilibrium populations f eq i are given in a closed form, where T 0 = 1/3: The Boltzmann BGK equation is solved: where τ is the relaxation parameter, and μ = τρT 0 is the resulting macroscopic dynamic viscosity. BGK stands for the Bhatnagar-Gross-Krook collision model as depicted in the right hand side of the aforementioned equation and describes the relaxation of populations f i to their equilibrium state f eq i with a single relaxation time τ. Implementation of the BGK model for porous medium flows needs specific care, since under specific circumstances, high slip velocities might arise at the solid-fluid interface. This could affect the numerical measurements of permeability and the evolution of the geometries due to reactions. A detailed study on this issue has been presented in [31]. Here, we operate the model in conditions as described in the aforementioned reference. For the advection-diffusion-reaction equations of the reactive species, a D2Q9 model is also implemented.
The equilibrium populations ξ i eq for the reactive species ξ i are given: where C i is the concentration of the considered ions, and u α is obtained from the basis model. The relevant population moment that corresponds to the concentration is Computer-Generated Stochastic Porous Media. Characteristic material microstructures can be resolved and digitalized, using a variety of experimental techniques that span from X-ray tomographic techniques to the combination of FIB-SEM microscopy. At the same time, there is a lot of effort in the algorithmic stochastic reconstruction of porous geometries (clays, membranes, etc.) that respect certain structural and statistical properties. In Figure 2, two such examples are shown. Figure 2(a) is obtained using the methodology of Tyagi et al. [32]. This pore map is constructed such that it matches the targeted rock's mineralogical components and its macroscopic properties such as porosity, grain, and pore size distributions. This methodology allows generating anisotropic structures, which are composed of different grain types. The map in Figure 2(a) represents a clay material, where grains and interlayers within the grains are present. The pore map in Figure 2(b) was created by following the methodology of [23,24]. A Voronoi tessellation of randomly distributed points on a plane is used as a template. Subsequently, the edges of the tessellation are used as guides to open the pore space, resulting in fully connected porosity. Here, we increase the degree of heterogeneity by varying the channel size and by distributing variable size spherical pores at the junctions of the pore channels. We construct our pore maps based on the latter methodology due to its simplicity and in order to work with a generic porous medium.

Description of the Reactive System and Boundary
Conditions. The evolution of a porous geometry, due to mineral reactions, depends strongly on the aqueous mixture species concentrations and the flow regime. The quantification of these conditions is aided from dimensional analysis. The nondimensional numbers that are used to characterize geochemical reactive flows at low Reynolds number (Darcy regime) are the Péclet number: Pe = UL/D 0 and the Damköhler number: Da = kL/D 0 C 0 , where D 0 is the mass diffusivity, U is the convective flow velocity, C 0 in this work is the acid concentration at inlet, and L is the characteristic linear dimension of the system of interest. Whenever mentioned in this text, the Pe and Da values correspond to the initial conditions of the flow setup.
The solid domain of the generated microstructures represents pure calcite. In the connected void space, a water solution is allowed to advect, diffuse, and react at the solid interface.
The reactant enters the domain from the left boundary with a uniform flow rate. The top and bottom boundary conditions are treated as periodic, while the right boundary as a zero gradient boundary.
For simplification, we here consider a single-step heterogeneous reaction of the dissolution of calcite applicable for low pH conditions: H + + CaCO 3 ➔ Ca 2+ + HCO − 3 . The reaction constant is k 1 = 8 9 × 10 −1 mol/m 2 s [33,34], and the reaction rate is a first-order reaction R = k 1 α H+ , where α H+ is the activity of H + in the solution. The diffusion coefficient is D 0 = 10 −9 m 2 /s, and the kinematic viscosity of the fluid ν = 10 −6 m 2 /s.

Results and Discussion
3.1. Permeability of Computer-Generated Structures. Several pore map realizations have been generated using the same set of rules (distribution of channel width and distribution of spherical pore sizes), differing only in the number of original random seed points. The two-dimensional domain size was selected to be square and is discretized by a 1500 × 1500 computational grid, whose total size corresponds to a 100 μm × 100 μm domain. This allows constructing random pore maps with variable target porosities that range from ε = 0.1 to 0.9. The permeability (K (m 2 )) of the generated structures was numerically measured by flow simulations similar to [31]. The results are plotted in the log-linear plot of Figure 3. The permeability correlates positively with the porosity. In the framework of macroscopic reactive transport codes, the microstructural information is implicitly considered using permeability-porosity correlations. In reactive environments, for example, when dissolution or precipitation processes take place, the porosity increases or decreases at a pace dictated by the flow and the chemical kinetics. In the case where the generated geometries are representing a real material, a first approximation, due to the lack of mechanistic understanding, would be to update the porosity according to the thermodynamics or kinetics. Subsequently, the permeability of the respective voxel would be updated following the permeability-porosity trend of the pristine structures (black circles in Figure 3). As it will be shown in the next section, such an approximation would Figure 1: The two-dimensional 9-velocity lattice (D2Q9).
3 Geofluids not be accurate for reactive transport calculations in the presence of strong chemical gradients.

Evolution of Permeability in Reactive Environments.
The use of the lattice Boltzmann framework to study the evolution of pore structures in reactive environments provides process understanding of the underlying mechanisms ([10, 12, and references within, 9,19,35]). The effect of different transport and chemical conditions has been studied by several authors [10,36]. Depending on the chemical conditions and the flow properties, the evolving geometries due to dissolution follow distinct paths. The change in pore topology and connectivity will have an effect on permeability. Different regimes have been identified and have been categorized in phase diagrams based on the characteristic nondimensional numbers that describe the reactive environment. For high Pe and Da numbers, the cross sections of the preferential flow paths increase fast, due to the strong supply of reactants (Pe) and the fast reaction rates (Da). Soon, the so-called wormholes start to appear resulting in the fast increase of the permeability versus the porosity of the reacting structure. At the other end, permeability increases at the slowest pace, when Pe is moderate and Da is high. The reason for that is the fast depletion of reactants at the incoming boundaries leading to face dissolution [12]. This categorization is verified also here.
We selected two different geometries stemming from the analysis of the previous section, with initial pristine porosity of ε = 0 25 (see Figure 4, top row) and ε = 0 39 (see Figure 4, bottom row). The evolution of each of the two pore  Figure 3: Permeability-porosity trends for several scenarios. Black-filled circles connected with a black line correspond to the permeability versus the porosity of the computer-generated unreacted samples. Blue (ε = 0 25) and red (ε = 0 39) curves correspond to the distinct permeability-porosity temporal evolution paths under different reactive transport conditions (circles: high Pe-high Da, crosses: moderate Pe-high Da). 4 Geofluids realizations was examined under two different flow conditions. In the first case, HCl acid 0.001 M (pH = 3) was injected from the left boundary resulting in Pe = 4 and Da = 87 × 10 6 . This case evolved as face dissolution. In the second case, HCL acid 0.01 M (pH = 2) was injected at higher rates resulting in higher Pe = 400 and Da = 87 × 10 5 . These conditions lead to the formation of wormholes. The evolution of the latter cases is shown in Figure 4. After specific time intervals, the evolved structures were extracted, and their permeability was measured in a similar manner as done in the previous section. The resulting permeabilities that represent the evolved pore structures are plotted in Figure 3. Blue lines and symbols represent the evolution of the ε = 0 25 structure, and red lines and symbols represent the evolution of the ε = 0 39 one.

Extracting Correlations from Pore-Level Simulations and
Upscaling Strategies. The pore-level simulations provide the basic understanding of the underlying mechanisms, which dictate the structural evolution of porous media in reactive conditions. At the same time, it is possible to provide the necessary input to the macroscopic algorithms via upscaling of the results. This can be simply done by replacing, for example, the permeability-porosity correlations that are modelled in a macroscopic code, with more precise and case-specific correlations. To that end, there are mainly three different approaches. First, the result can be transmitted in the form of power law or Kozeny-Carman type of function as K = f ε . Second, tabulated values can be provided instead of a power law, such that during the macroscopic simulations, specific values can be calculated, after interpolating between successive points. Third, the macroscopic code could call on demand the pore-level solver, to deliver the prediction of evolution in a fully coupled multiscale manner. The permeability porosity correlations that describe the evolution of the structures in the presence of velocity and chemical gradients have been extracted after fitting known type of relationships. As a first remark, it should be mentioned that it was not possible to fit the results using as a model of Kozeny-Carman-type function. This was especially true for the case of high Pe-Da (wormholes). The most appropriate and simple type of model was found to be the K = a × ε b , where a and b are the fitting parameters. The results are shown in Figure 5 for the pore map ε = 0 39 (see Videos S1 and S2 in the Supplementary Materials, which depicts the temporal evolution of the porous medium and the respective mass transport).
The scaling with respect to the porosity was found to be of the order b = 0 9 for the moderate Pe case. The high Pe-Da results could not be fitted using only one function of that type due to the existence of two distinct evolution regimes. Evolution starts with a relatively weak scaling, where b = 2 4 (for 0.39< ε <0.45) and corresponds to the wormhole formation, and continues with a strong scaling regime where b = 13 7 (for 0.45< ε <0.5) and corresponds to the point where the flow is accelerated due to the growing of the major wormholes. For the upscaling of these results, specific care has to be taken that the simulation At the moment, it is seems that such a fitting cannot be fully automated in the classical sense. The existence of (a) very small numbers among the fitted parameters and (b) different scaling regimes, will always be a source of error when transferring such functionals to macroscopic codes. Nevertheless, if done with care, this kind of microscopic feedback has the potential to substantially improve the modelling perspectives, provided that digitalized characteristic microstructures are available. By passing, we note that the two geometries (ε = 0 25 and ε = 0 39) seem to evolve in a similar way under similar conditions (see Figure 3). The comparison between the blue and the red set of curves in Figure 3 goes beyond the scope of this paper. For the extraction of safe conclusions, it would be needed to test several geometries in reactive environments, and proceed with statistical analysis of the results.
The second approach would be to tabulate the a priori executed pore-level simulation results and provide these tables as input to a macroscopic code. Exact values for specific cases can be determined via interpolation schemes. This would not require a fitting process and shall be efficient as soon as a minimum number of information (size and accuracy of tables) has to be passed to the next level of description. We note that for complex chemically reactive systems, this option would require the precalculation of a large number of different scenarios. The subsequent storage of the resulting information with very high accuracy could become impractical.
The third option would be to fully couple the pore-level codes with macroscopic codes in a multiscale manner. Executing several pore-level simulations for every voxel, in order to extract the microscopic physics, would defeat the purpose of pore-level simulations. At the same time, it would result in an extremely slow simulation. Nevertheless, novel much promising machine learning techniques started recently to appear as, for example, in [37]. In that paper, the authors demonstrate a technique to accelerate the solution of chemical equilibrium equations, by using previous equilibrium calculations of similar input conditions. Such algorithms, when applied to the problem described in this paper, can identify voxels with similar input and save the most representative cases, for use in future calculations. This would accelerate a fully coupled multiscale code.

Conclusions
The evolution of permeability with respect to porosity in reactive environments is of great importance in geochemical reactive transport modelling and subsequent predictions. Within the porous structure, the mineral dissolution and precipitation processes alter the pore topology in a strongly nonlinear way. In this paper, we have constructed calcite microstructures in a stochastic way and have traced their temporal evolution in the presence of an acid. Permeability-porosity correlations have been extracted for each of the different evolution paths. Using this paradigm, we discussed possible upscale strategies and the numerical bridging between pore level understanding and macroscopic modelling. We conclude that such approaches are viable under the condition of fully automating the communication between the solvers.

Data Availability
Data can be requested from the authors, and all information necessary to reproduce results are mentioned within the article.  Figure 5: Fitted permeability porosity correlations, describing the evolution of permeability of the pore structure ε = 0 39, for two different reactive and flow conditions. Open red circles correspond to the evolution under high Pe-high Da conditions and have been fitted using two correlation functions (black lines). Red solid line and crosses depict the low Pe-high Da case. Black circles represent the permeability of the pristine structures. 6 Geofluids