Reactive Transport Model of Gypsum Karstiﬁcation 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 karstiﬁcation 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 karstiﬁcation of: (1) the dynamic update of ﬂow and transport parameters due to porosity changes; (2) the spatial distribution of minerals in the rock matrix; (3) the time evolution of water inﬂows through the boundaries of the model; (4) the functions relating permeability, k , to porosity, φ . The average porosity of the matrix after 1000 years of simulation increases from 0.045 to 0.29 when ﬂow, transport, and chemical parameters and the water inﬂows 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 inﬂow 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 modiﬁed Fair-atch k- φ relationships. The gypsum dissolution front propagates into the matrix faster when the power law with n = 2 and 3 and the Verma–Pruess k- φ relationships are used.


Introduction
Gypsum is one of the most common evaporitic rocks. Its high solubility compared to other sedimentary rocks can lead to the formation of karst features [1]. Its dissolution generally causes the development of fissures and enlarges cavities, providing preferential flow paths in which groundwater moves more rapidly than in the rest of the medium. Karst features formed in evaporitic rocks are similar to those found in carbonate rocks. However, the main difference is in the time it takes for karst features to form in both sedimentary rocks. While their development in natural carbonate aquifers could take thousands of years, fractures in gypsum and other evaporitic rocks may evolve on a human timescale [2]. Therefore, detailed knowledge about the evolution of karstic systems containing gypsum is key to many human activities, including freshwater supply, agriculture, waste disposal, and construction projects, as well as in the avoidance of possible catastrophes caused by subsidence [3][4][5][6][7][8].
Evaporitic formations underlie around 25% of the Earth's continental surface, being formed by the evaporation of saline waters causing the crystallization of dissolved salts [9].
Traditional groundwater flow codes are based on Darcy's law, which is valid when the flow is laminar. However, when water flow through fractured media is turbulent, Darcy s law is no longer valid. Scanlon et al. [59] evaluated the ability of equivalent porous media models in fractured media and reported that these models could be used to simulate regional groundwater flow in highly karstified aquifers. By contrast, Karay and Hajnal [60] and Fang and Zhu [61] reported that the incorrect use of the laminar or turbulent flow equation in the fractures can lead to significant errors in numerical simulation, although the importance of these potential errors depends largely on the aperture of the fractures.
The numerical modeling of fractured porous media in carbonate aquifers has been performed extensively [27,[62][63][64][65][66]. However, the numerical simulation of evaporitic karst aquifers has been less common. Birk [67] developed a process-based modeling tool for the characterization of highly complex karst flow systems and applied it to gypsum aquifers. Birk et al. [68] simulated the development of gypsum maze caves under artesian conditions. Campana and Fidelibus [69] applied a combined reactive-transport/density-dependent flow model to evaluate gypsum dissolution in a coastal kart aquifer. Guo et al. [70] described the formation of a gypsum cavity induced by dissolution through a large-scale diffuse interface model.
Here, we present a reactive transport model to study gypsum karstification in physically and chemically heterogeneous fractured evaporitic media and evaluate the long-term (1000 years) evolution of gypsum karstification. The models are based on the hydrogeological and hydrogeochemical conditions prevailing near Villar de Cañas in central Spain, where a site was investigated for its suitability for hosting an interim surface facility for high-level radioactive waste [71,72]. Numerical models account for the changes in porosity caused by mineral dissolution/precipitation and the associated effects on the flow, transport, and chemical parameters of the fractured medium. Several scenarios are analyzed by considering several spatial patterns of mineral distribution and several distributions of boundary water inflows. In addition, gypsum karstification is evaluated for six different permeability-porosity relationships.

Conceptual Model
The conceptual model was established based on the hydrogeological and hydrogeochemical conditions near the village of Villar de Cañas (Cuenca, Spain) within the Loranca Basin [71][72][73]. The basin contains mainly continental Tertiary deposits folded along with the underlying Mesozoic rocks [74]. The main hydrogeological units in the Villar de Cañas area (from bottom to top) include ( Figure 1): (1) the lowermost hydrogeological formations, which include the Lower Balanzas Lutite Formation (LBI) and the Lower Tertiary Formation (UI); (2) the evaporitic Balanzas Gypsum Formation (YB); (3) the Upper Balanzas Lutite Formation (LBS). These formations show large spatial heterogeneity. The UI + LBI formation contains lutites/shales, gypsiferous marls, and sandstones. The YB formation is made up of massive gypsum with intercalations of lutites and limestones, while the LBS formation is composed mainly of gypsiferous lutites [75]. Most of the formations in the Villar de Cañas area have low to very low hydraulic conductivities, with minimum values of about 10 −6 m/d in massive gypsum. However, the hydraulic conductivity can reach high values (>10 m/d) in areas affected by karstic dissolution in the YB formation [76]. The main water inflow into the groundwater system comes from the infiltration of rainfall in the outcrops of the UI+LBI formation east of the Villar de Cañas area. The groundwater flow shows a regional component from this area to the west, discharging into the Záncara River and its alluvial, as well as into a creek located in the outcrop of the YB formation [75]. Geological map of the Villar de Cañas area showing the location of boreholes S1 and S2 and the Balanzas Gypsum formation (YB), which are affected by gypsum dissolution due to the inflow of groundwater that is undersaturated with respect to the gypsum from the lowermost LBI (Lower Balanzas Lutites) and UI (Lower Tertiary) formations.
Calcium, magnesium, and sulphate are the dominant dissolved ions in the groundwaters of the study area. The pH ranges from 6.8 to 7.3, while the electrical conductivity ranges from 2000 to 10 5 µS/cm. The chemical composition of the groundwater in the Villar de Cañas area is the result of the interaction of the groundwater with several reactive minerals present in the subsurface. Gypsum is the most abundant mineral in the system, with which water reacts immediately. The groundwater from the UI + LBI formation is undersaturated with respect to gypsum and, in its discharge towards the creek in the YB formation, this mineral dissolves, causing the development of discontinuities and small fractures. The dissolution of gypsum and carbonates increases the concentration of dissolved calcium and sulfate until water reaches equilibrium with gypsum. However, the high concentrations of sulfate found in the groundwater reveal that there are other sources of sulfate in the groundwater, such as the dissolution of other sulfate minerals found in small amounts at the site, such as bloedite Na 2 Mg(SO 4 ) 2 ·4H 2 O, thenardite, Na 2 SO 4 , and glauberite, Na 2 Ca(SO 4 ) 2 [76][77][78][79]. In these low-permeability media, the groundwater flows preferentially through the discontinuities created through mineral dissolution. A reactive transport model was used to evaluate the gypsum karstification along the contact of the YB and UI+LBI formations. A discontinuity/fracture was included at the bottom of the 2 × 2 m 2 two-dimensional model domain, while the rest of the domain corresponded to the rock matrix ( Figure 2). The thickness of the model was uniform and equal to 0.1 m. Although the thickness did not play any role in the model's results because the model was two dimensional, it was taken into account to calculate the water inflow through the boundary. Here, the water inflow into the model domain was equal to 742.78 L/year based on the groundwater flow models used by Águila [80] in the Villar de Cañas area. Initially, the water inflow into the model occurred only through the discontinuity from east to west (right to left in Figure 2). The water flowing into the model domain was undersaturated with respect to gypsum and caused its dissolution. The geochemical model accounted for the following reactions: (1) aqueous complexation, (2) acid/base, and (3) mineral dissolution/precipitation. The reactive transport model was used under isothermal conditions at 25 • C. The chemical system was defined in terms of the concentrations of the following primary species: H 2 O, H + , Ca 2+ , Cl − , HCO 3 − , K + , Mg 2+ , Na + , SO 4 2− , and SiO 2 (aq). The model considered 7 minerals and 35 aqueous species identified from speciation runs performed with the EQ3/6 code [81]. Chemical reactions and their equilibrium constants at 25 • C for aqueous species and mineral dissolution/precipitation are listed in Table 1. Isothermal reactive transport simulations are usually performed at 25 • C because the available thermodynamic data in existing thermodynamic databases correspond to 25 • C. The temperatures of groundwater in the Villar de Cañas area range from 14 • C at the water table to 20 • C at 200 m depth. Although the equilibrium constants change with temperature, their values at 25 • C are assumed to be representative of the conditions of the study area. All of the reactions except for the dissolution/precipitation of magnesite are assumed to be at chemical equilibrium. The kinetic dissolution/precipitation of magnesite is modeled with the following kinetic rate law: where r is the dissolution/precipitation rate (mol/m 2 /s); s is a function that is equal to 1 for mineral dissolution and −1 for precipitation; Ω is the saturation index, which is equal to ratio of the ion activity product and the equilibrium constant; e −Ea/RT is a thermodynamic factor that takes into account the apparent activation energy of the reaction Ea (KJ/mol), the gas constant R (KJ/mol·K), and the absolute temperature T (K); k is the kinetic rate constant (mol/m 2 /s) at 25 • C; θ and η are empirical parameters. The kinetic parameters of magnesite were obtained from Palandri and Kharaka [82], while the specific surface area was assumed to be equal to 0.6987 dm 2 /L.

Numerical Model
A 2-D discrete fracture matrix (DFM) groundwater flow and reactive transport model under saturated conditions was used to study the dissolution of gypsum and evolution of a discontinuity/fracture. The mathematical model was created with the finite element code of CORE 2D V5 [21,84,85]. This is a code for transient saturated and unsaturated water flow, heat transport, and multicomponent reactive solute transport under both local chemical equilibrium and kinetic conditions in heterogeneous and anisotropic media. CORE 2D V5 has been extensively verified against analytical solutions and other reactive transport codes [30,86,87] and is widely used to model groundwater flow and solute transport in real-world aquifers [80,88,89] and laboratory or in situ experiments [90,91], as well as to study the long-term geochemical evolution of radioactive waste repositories [92][93][94]. Moreover, the code is able to consider the porosity feedback effect caused by mineral dissolution/precipitation reactions and update the flow and transport parameters at every time step, thus enabling the investigation of the spatial and temporal evolution of fractures. CORE 2D V5 uses the Kozeny-Carman equation to calculate the changes in permeability (k) due to porosity (φ) variations according to [23]: where k 0 and φ 0 are the initial or reference permeability and porosity, respectively. The model domain has a surface of 2 × 2 m 2 and is discretized with a 2-D finite element mesh of triangular elements ( Figure 2). The mesh has 1845 nodes and 3520 elements and includes two material zones: (1) a matrix composed mainly of gypsum and (2) a 0.05-mwide discontinuity at the bottom boundary of the model. Several scenarios are considered by varying the spatial distribution of the minerals in the matrix. The hydrodynamic and transport parameter values used in the matrix [95,96] and the discontinuity of the numerical model are listed in Table 2. The model parameters were derived from the flow, heat transfer, solute transport, and groundwater age models developed by Águila [80] in the Villar de Cañas area. The hydraulic conductivity (K) of the fracture is more than six orders of magnitude larger than that of the matrix. The porosity in the matrix is equal to 0.045, while that of the fracture is assumed to be equal to 0.95. Numerical simulations were performed for a total time of 1000 years. Table 2. Hydrodynamic and transport parameters of the matrix and the discontinuity (fracture) of the 2-D reactive transport model (adapted from [80,95,96] Figure S1 of the electronic Supplementary Materials (ESM) shows the finite element mesh with the flow boundary conditions in the numerical model. The water outflow occurs through the left boundary of the model and is simulated with a Dirichlet condition by assuming a prescribed hydraulic head (H = 0 m). The water inflow (Q) to the model occurs through the right boundary and is simulated with a Neumann condition in which Q at each node depends on its porosity. CORE 2D V5 was updated to allow for the dynamic calculation of the water inflow as a function of the porosity, Q(φ), and to enable Q to be distributed along the right boundary of the model as φ increases due to the dissolution of gypsum and other minerals present in the system. The modified version of the code ensures that the total water inflow remains constant over the 1000 years of the simulation.
The modified code only allows for a water inflow through the nodes of the right boundary, in which the computed porosity exceeds a threshold value that is assumed to be equal to 0.05. The initial porosity in the nodes of the matrix is equal to 0.045, so initially, the water inflow to the model occurs only through the discontinuity. At each time step, the code also calculates the sum of the computed porosities in the nodes located at the right boundary of the model with φ larger than 0.05 to create a weight that allows the distribution of the water inflow through the right boundary. The water inflow rate through each node of the right boundary, Q i , is calculated according to: where Q i and φ i are the water inflow rate and porosity for node i, respectively; Q t is the total water inflow into the model domain; n is the number of nodes along the right boundary of the model domain with porosities larger than 0.05. Two initial porewaters and a boundary water with different chemical compositions were defined in the reactive transport model. The chemical compositions of these porewaters are listed in Table 3 and were obtained from chemical speciation by using chemical analyses of water samples collected in two boreholes [77][78][79]. The chemical analyses of water from the S2 borehole, drilled in the YB formation (see Figure 1), were used to define the chemical composition of the initial porewater in the matrix. Instead, the chemical composition of the initial water in the discontinuity and the inflow water through the right boundary of the model were defined from the chemical analyses of water samples collected in the S1 borehole, which was drilled in the UI + LBI formation. Table 3. Chemical compositions of the boundary water and the initial porewaters used in the reactive transport model (values in mol/L). Four different scenarios (1, 2, 3, and 4) are analyzed to study gypsum karstification for several mineral patterns and distributions of boundary water inflows. In addition, gypsum karstification is evaluated in scenario 1 for six different permeability-porosity relationships.

Scenario 2
The mineral composition of the matrix in scenario 2 is similar to that of scenario 1 ( Figure 2). However, the water inflow is not distributed along the right boundary in scenario 2, so the water inflow into the model occurs only through the discontinuity during the entire simulation. Figure S2 of the electronic Supplementary Materials shows the flow boundary conditions of scenario 2. The model's results in this scenario are computed by taking into account the dynamic update of flow and transport parameters due to changes in porosity and by assuming that porosity and the rest of the model parameters remain constant in time. In the last procedure, the porosity variations are calculated externally from mineral volume fraction variations, which is a common modeling practice to greatly reduce computation times and convergence issues [20]. The results obtained in both cases are compared to evaluate the significance of the dynamic update of the flow, transport, and chemical parameters in reactive transport models with mineral dissolution/precipitation.

Scenario 3
The flow boundary conditions imposed in the model used in scenario 3 are those defined in Figure S1 (similar to scenario 1). However, a stratified model with four mineral zones in the matrix is considered in scenario 3. Figure 3 shows the finite element mesh and the geometry of the four mineral zones defined in this scenario. Three additional zones with different mineral compositions have been introduced into the matrix, forming horizontal bands parallel to the fracture that cover the entire length of the model. The mineral volume fractions are listed in Table 4. The volume fractions of the highly soluble sulfate minerals (bloedite, thenardite, and glauberite) increase when moving away from the discontinuity. Bloedite, thenardite, and glauberite are not present in the matrix (zone 1).  The model used in scenario 4 is similar to that in scenario 3 (stratified model and flow boundary conditions ( Figure S1)), except for the location and length of mineral zones 2, 3, and 4. These zones are horizontal bands that are 1 m long and located in the center of the matrix (see Figure 3). Unlike scenario 3, mineral zones 2, 3, and 4 do not reach the model boundaries. The mineralogical compositions of the material zones are listed in Table 4. Figure 4 shows the spatial distribution of the computed porosity, φ, hydraulic conductivity, K, and flow velocity after 100, 300, and 1000 years of simulation. The water inflow into the model domain is undersaturated with respect to the gypsum, causes the progressive dissolution of this mineral in the matrix, and causes the increase in φ and K near the right boundary of the model. Initially, the water inflow occurs only through the fracture at the bottom of the model, so the gypsum dissolution in the matrix is limited to the surroundings of the fracture. However, after about five years of simulation, the first node of the matrix at the right boundary exceeds the threshold porosity of 0.05, enabling the inflow of water into the model through nodes of the matrix as the fracture increases in size. After 100 years of simulation, the computed porosity in all of the nodes of the finite element mesh at the right boundary exceeds the threshold porosity of 0.05, so the total water inflow into the model is distributed (proportionally to the porosity of each node) along the entire right boundary of the model. Nevertheless, there are no significant variations in porosity after 100 years, since the average porosity of the matrix remains very low (0.071) and the maximum porosity is smaller than 0.18. As the simulation progresses, the porosity near the right boundary of the model increases. After 300 years of simulation, the porosity in some nodes of the matrix reaches values larger than 0.86, and the average porosity of the matrix increases to 0.143. The computed porosity exceeds 0.85 in almost 21% of the matrix at the end of the simulation, leading to a significant increase in the size of the original fracture. The average porosity at the nodes considered initially in the matrix is greater than 0.29 after 1000 years of simulation. CORE 2D V5 updates the hydraulic conductivity and pore diffusion coefficient at every time step according to the porosity variations by using the Kozeny-Carman equation (Equation (2)) and Archie's law [97], respectively. Therefore, the spatial distributions of both parameters are closely related to the spatial distribution of the porosity (see Figure 4 and Figure S3 of the electronic Supplementary Materials). The average of the hydraulic conductivity in the matrix increased by five orders of magnitude during the 1000 years of simulation, thus affecting the spatial distribution of the flow velocities in the model. Initially, the flow velocity through the matrix is 3000 times smaller than the velocity in the fracture. However, as water inflow occurs through nodes of the matrix along the right boundary of the model, the flow velocity in the matrix increases and, after 100 years, it is greater than 10 −7 m/s in a large part of the matrix (see Figure 4). The groundwater in the matrix moves down towards the fracture from right to left as the porosity and hydraulic conductivity increase near the right boundary of the model. The increase in K near the right boundary of the model facilitates the groundwater to flow towards the fracture and causes a decrease in flow velocity in the rest of the matrix. The area of the matrix with low water velocities (<1.6 × 10 −8 m/s) increases with time, taking up 76% of the matrix surface after 1000 years of simulation.

Scenario 1
The increase in porosity with time in the model matrix is mainly caused by the dissolution of gypsum. Figure 5 presents the spatial distribution of the computed mineral volume fraction and porosity along a vertical line (x = 1.7 m) after 100 and 1000 years of simulation. Gypsum shows a dissolution front that is displaced through the matrix with time, causing the porosity to increase. The rest of the minerals considered in the model have a very small effect on the porosity. Chalcedony dissolves slightly and calcite precipitates in areas where gypsum has completely dissolved. Magnesite barely reacts during the simulation, while bloedite, thenardite, and glauberite dissolve completely during the first times of simulation because they have high solubility and reactivity. The dissolution of these highly soluble sulfate minerals, which are found in very small amounts in the Villar de Cañas area [98], causes the large concentrations of dissolved sulfate measured in the groundwaters.   Figure 6 shows the spatial distribution of the computed concentration of dissolved SO 4 2− at 0.01, 9, 100, 500, and 1000 years of the simulation in scenario 1. The computed concentration of dissolved SO 4 2− increases throughout the matrix at initial times (>1 mol/L) due to the dissolution of bloedite, thenardite, and glauberite (see Figure 6). Afterwards, the sulfate concentrations in the matrix decrease as water with a lower sulfate concentration enters into the model (t = 9 years). However, the concentration of dissolved SO 4 2− in the matrix is larger than that of the boundary water because of the dissolution. After 100 years of simulation, the concentration of dissolved SO 4 2− is equal to 0.018 mol/L throughout the matrix. The sulfate concentrations remain constant until gypsum begins to be exhausted at the right boundary of the rock block. Then, the computed concentration of dissolved SO 4 2− decreases until reaching values similar to those of the boundary water (0.0067 mol/L). The decrease in flow velocity over time at the right boundary of the rock block causes an increase in sulfate concentrations in this area, as attested by the spatial distribution of the flow velocity and the concentration of dissolved SO 4 2− at 1000 years in Figures 4 and 6, respectively. The dissolution of gypsum in the matrix also causes significant variations in the concentrations of dissolved Ca 2+ . Figure 7 illustrates the time evolution of the computed concentrations of dissolved Ca 2+ and SO 4 2− and the volume fraction of gypsum at the mesh of three nodes. They are located at the same distance from the discontinuity (y = 0.325 m), but at different distances from the right boundary of the model (x = 2 m, x = 1.75 m, and x = 1.65 m). Initially, the computed concentrations of dissolved SO 4 2− increase due to the dissolution of bloedite, thenardite, and glauberite. Then, the dissolution of gypsum causes an increase of the concentration of dissolved Ca 2+ . Afterwards, the calcium and sulfate concentrations remain nearly constant, while gypsum is dissolved. The computed concentrations of dissolved Ca 2+ and SO 4 2− decrease to values similar to those in the boundary water when gypsum has completely dissolved. Gypsum is exhausted at the nodes at x = 2, 1.75, and 1.65 m after 400, 600, and 900 years, respectively. Therefore, the gypsum dissolution rate decreases over time. The 25 cm band of gypsum closest to the right boundary of the model dissolves at a rate of 1.25 mm/year, while the dissolution rate in the neighboring 10 cm band is 0.33 mm/year. The average gypsum dissolution rate in the entire model domain during the 1000 years of simulation is 0.375 mm/year. Figure 8 shows the spatial distribution of the porosity after 100, 500, and 1000 years of simulation in scenario 2 computed with and without the porosity feedback effect. When the water inflow into the model occurs through the fracture and the flow and transport parameters are updated by taking changes in porosity into account, gypsum dissolves near the fracture. Therefore, the fracture width increases with time by up to three times after 1000 years of simulation. The porosity is greater than 0.85 in 2.3% of matrix after 500 years of simulation, while this value increases to 5.9% of the matrix after 1000 years. Moreover, the average porosity in the matrix is equal to 0.13 at the end of the simulation. On the other hand, gypsum hardly dissolves around the fracture when the porosity feedback effect is not taken into account during the simulation. The average porosity in the matrix, calculated externally from changes in mineral volume fractions, increases only up to 0.05 after 1000 years of simulation. It is clear, therefore, that the type of flow boundary condition imposed on the model and the porosity feedback effect are especially relevant in the dissolution of gypsum and the evolution of the fracture width in the long term.

Scenario 3
This scenario addresses gypsum karstification by considering several spatial patterns of mineral distribution in the matrix that may introduce additional discontinuities. The dissolution of bloedite, thenardite, and glauberite during the initial stages of the simulation causes a rapid increase in porosity in zones 2, 3, and 4 ( Figure 3). The porosity in these areas exceeds the threshold value (φ > 0.05), so the code enables, from early simulation times, the inflow of unsaturated water with respect to gypsum directly into the matrix through the nodes at the right boundary of the model in zones 2, 3, and 4. The increase in porosity and, consequently, the water boundary inflow into the model domain are larger in the zones of the model where the amounts of bloedite, thenardite, and glauberite are larger (zone 4 in Figure 3 and Table 4). The dissolution of these minerals along with the onset of the gypsum dissolution front in the matrix leads to an increase in the concentrations of dissolved Ca 2+ and SO 4 2− , subsequently causing the precipitation of gypsum in zones 2, 3, and 4 of the model. Gypsum precipitation in zones 2-4 occurs during the first two years of simulation. This precipitation occurs slightly earlier in zone 4, where the porosity and permeability are larger than those in zones 2 and 3 due to the greater initial volume fractions of bloedite, thenardite, and glauberite in zone 4. Figure 9 shows the spatial distribution computed for the cumulative gypsum, bloedite, thenardite, and glauberite dissolution/precipitation along a vertical line at x = 1 m after 2 years of simulation. Bloedite, thenardite, and glauberite dissolve while gypsum precipitates in zones 2, 3, and 4. In addition, gypsum dissolves near the fracture at the bottom of the model domain.  Figure 10 presents the spatial distribution of the computed porosity after 100, 300, and 1000 years of simulation in scenario 3. The porosity increases more rapidly in the zones where the highly soluble sulfate minerals are present at the start of the simulation, leading to the development of additional discontinuities/fractures in the matrix. The average porosity in the matrix increases to 0.1 and 0.2 after 100 and 300 years of simulation, respectively. At the end of the simulation, the computed porosity exceeds 0.85 in 24.7% of the matrix, and its mean value in the matrix increases to 0.343 (more than 7.5 times greater than the initial porosity). Even though bloedite, thenardite, and glauberite are exhausted at the beginning of the simulation, the distribution of these minerals in the model domain has a significant impact on gypsum dissolution and long-term porosity changes. Figure 11 shows the spatial distribution of the computed gypsum volume fraction and the porosity along a vertical line in the middle of the model (x = 1 m) at 0, 100, and 1000 years. The slight increase in the volume fractions of highly soluble sulfate minerals in zones 2 (3% vol.), 3 (6% vol.), and 4 (9% vol.) has a significant impact on the porosity reached in these zones after 1000 years. The maximum porosity at the end of the simulation in zone 4 at 1 m from the right boundary of the model is around 0.45, and it is 2.5 and 4 times larger than those of zones 2 and 3, respectively. In addition, the porosity of the matrix near zone 4 also increases significantly, which can lead to a long-term connection between the fractures, thus favoring gypsum dissolution.

Scenario 4
This scenario addresses the effect on gypsum karstification of the presence in the matrix of areas with different mineral compositions that are not in contact with the water inflow and outflow boundaries of the model (Figure 3). For this, the numerical model assumes that the highly soluble sulfate minerals can dissolve in the early stages of the simulation, causing the porosity in zones 2, 3, and 4 to increase slightly with respect to that in the matrix. However, the rate of water inflow into the model is not affected by the dissolution of bloedite, thenardite, and glauberite. Figure 12 shows the spatial distribution of the computed porosity after 100, 300, and 1000 years of simulation. After 100 years, the computed porosity increases in zones 2, 3, and 4 of the model due to the dissolution of the highly soluble sulfate minerals. Nevertheless, there are no changes in porosity yet at the right boundary of the model, and the average porosity of the matrix increases only to 0.053. In contrast, the porosity begins to increase near the right boundary of the model after 300 years of simulation. In addition, the increase in porosity to the right of zones 2, 3, and 4 (near the inflow boundary of the model) is greater than in the left boundary. For instance, the computed porosity to the right of zone 3 increases to 0.24, while φ to the left is equal to 0.14 after 300 years of simulation. The average porosity in the matrix increases to 0.1 at this time. At the end of the simulation, the computed porosity is larger than 0.85 in 22% of the matrix, and the average porosity in the matrix increases to 0.295. The gypsum dissolution front reaches zones 2 and 3 after 1000 years. The computed porosity increases more rapidly in these zones than in the rest of the matrix.  Figures 13 and 14 illustrate the spatial distribution of the concentration of dissolved SO 4 2− and Mg 2+ at selected times of the simulation in scenario 4, respectively. The initial concentration of dissolved SO 4 2− in the matrix is about 0.025 mol/L, but the dissolution of bloedite, thenardite, and glauberite in zones 2, 3, and 4 causes a significant increase in the computed sulfate concentrations in these zones after 1 year of simulation. Dissolved sulfate diffuses through the matrix, reaching values of about 2 mol/L in most of the matrix after 10 years of simulation. The water inflow into the matrix increases due to the dissolution of gypsum, causing water with lower sulfate concentration to ingress into the matrix (t = 50 years). After 100 years of simulation, the concentration of dissolved SO 4 2− is approximately uniform throughout the matrix. However, the computed concentration of dissolved SO 4 2− decreases to the concentration of the boundary water in areas where gypsum has been exhausted (t = 500 and 1000 years).
The spatial distribution of the computed concentration of dissolved Mg 2+ is similar to that of dissolved sulfate. Initially, the concentration of dissolved Mg 2+ in the matrix is about 0.0055 mol/L. However, the dissolution of bloedite at early times causes the increase in the computed magnesium concentrations in zones 2, 3, and 4. Then, dissolved magnesium diffuses through the matrix (t = 1 and 10 years). The dissolution of gypsum leads to the increase in the water inflow into the matrix and the decrease in the concentration of dissolved Mg 2+ (t = 50 years). The concentration of dissolved Mg 2+ in the matrix is constant and similar to the concentration in the boundary water from 100 to 1000 years.

Analysis of Sensitivity to Permeability-Porosity Relationships
Sensitivity runs were performed to investigate the relevance of the permeabilityporosity (k-φ) relationships in gypsum karstification. The following three additional k-φ relationships were implemented in CORE 2D V5: (1) the modified Fair-Hatch equation, (2) the Verma-Pruess model, and (3) the power-law relationships. The simulation in scenario 1 was performed for several k-φ relationships.
The modified Fair-Hatch equation is commonly used to describe the dissolution of minerals. It was originally derived by Bear [99] and later modified by Chadam et al. [31] as follows: where φ f is the final porosity after complete dissolution of the soluble mineral and E 1 is a constant. Here, E 1 was assumed to be equal to 0.7 [100]. The Verma-Pruess equation was derived by Verma and Pruess [32] from a pore-body and throat model as: where φ c is the critical porosity at which the permeability becomes zero, and n is an empirical parameter. Verma and Pruess [32] reported that φ c should be around 0.8-0.9 of the original porosity. Here, we take φ c = 0.8φ 0 and n = 2 [34]. Power-law relations are widely used to model permeability changes in experiments. Power-law relationships have the following general form: Unlike the Kozeny-Carman equation, the exponent n is usually calibrated from experimental data, which makes the power-law relationships slightly more adaptive. However, they are often considered less accurate due to oversimplification of the complex processes that change pore structure [35]. A wide range of n values have been reported in experimental and modeling studies [101][102][103][104]. Here, we performed sensitivity runs for n = 2, 3, and 5. Figure 15 shows the computed porosity and hydraulic conductivity along the horizontal line at y = 1 m after 1000 years of simulation in scenario 1 by using the six k-φ relationships. Additionally, the spatial distributions of the computed porosity, hydraulic conductivity, and flow velocity after 100, 300, and 1000 years of simulation using the modified Fair-Hatch ( Figure S4), Verma-Pruess ( Figure S5), and power-law (from Figures S6-S8) permeability-porosity relationships are presented in the electronic Supplementary Materials. The gypsum dissolution front is affected by the type of permeability-porosity relationship. At the end of the simulation, the size of the matrix where gypsum has been exhausted is greater for the Verma-Pruess and power-law (with n = 2) relationships. The gypsum dissolution front moves more slowly with the Kozeny-Carman, modified Fair-Hatch, and power-law (with n = 5) relationships. The average porosity in the matrix after 1000 years of simulation is equal to 0.373 for a power-law (with n = 2) relationship, while it is less than 0.3 for the Kozeny-Carman, modified Fair-Hatch, and power-law (with n = 5) relationships. The differences in the model predictions increase with time. The average porosities in the matrix after 100 years in the sensitivity runs range from 0.068 to 0.073. The hydraulic conductivities at the end of the simulation are largest for the permeabilityporosity relationships with the largest exponents n (power law with n = 5, Kozeny-Carman, and modified Fair-Hatch relationship). The changes in k affect groundwater flow, gypsum dissolution, and, consequently, the changes in matrix porosity. Figure 16 shows the spatial distribution of the computed flow velocity after 500 years of simulation in scenario 1 using the six permeability-porosity relationships. The spatial distributions of flow velocity in the models using the power law with n = 5, Kozeny-Carman, and modified Fair-Hatch relationships are similar at this time. The higher values of hydraulic conductivity achieved in these three models cause water to flow preferentially from the right boundary of the model towards the fracture at the bottom of the model, while the computed flow velocity remains at low values (less than 2 × 10 −8 m/s) in a large part of the matrix (blue areas in Figure 16). Consequently, the gypsum dissolution front and the increases in porosity progress more slowly in these models with higher values of n. The area of the matrix with flow velocities greater than 2 × 10 −8 m/s after 500 years of simulation is larger in the models using the Verma-Pruess (Figure 16c) and power law with n = 3 ( Figure 16e) and n = 2 (Figure 16f) relationships than in the other three models.

Discussion
The investigation of karstification processes in evaporitic formations has not received as much attention as karstification in carbonate rocks due to their poorer water quality in evaporitic media [15,105,106]. However, understanding the evolution of dissolution and karstification processes in these formations (much faster than in carbonate aquifers) is essential to avoid the collapse of structures and catastrophes [107,108]. For this, numerical modeling plays an important role due to the excessive research times required for the study of karst aquifers in the laboratory and in the field [109]. The current investigation demonstrated the importance of reactive transport modeling for understanding and simulating gypsum karstification through a discrete fractured matrix. The use of numerical models capable of dynamically updating the flow and transport parameters as minerals dissolve proved fundamental for a realistic representation of the enlargement of discontinuities and fractures. Table 5 compares the average porosities in the matrix at selected times and the percentages of the matrix volume with porosities larger than 0.85 at the end of the simulation in the four scenarios analyzed by using different modeling hypotheses. The simulations carried out in scenario 2 without considering the porosity feedback effect (PFE) show very small variations in porosity over the 1000 years of simulation. On the other hand, the average porosity of the matrix at the end of the simulation for the same scenario is almost three times greater than its initial value when the PFE is taken into account. Furthermore, the use of numerical models that account for the PFE allows the dynamic updating of the boundary conditions of the model when the boundaries are subject to changes in porosity due to mineral dissolution/precipitation. The average porosity computed in the matrix after 1000 years in the model of scenario 1 with the update of the inflow boundary is 0.29. In this case, 20.96% (vol.) of the matrix has a porosity greater than 0.85 at the end of the simulation. Thus, the dynamic update of the water inflow boundary according to porosity changes (scenario 1) causes an increase in the average porosity of the matrix of 121% with respect to the model considering the PFE but without the update of the inflow boundary (scenario 2). It is well known that the relevance of the dynamic update of flow and transport parameters increases with time, and therefore, the PFE is especially relevant for long-term simulations. While the average porosity of the matrix increases by 45% after 100 years when the PFE is taken into account, the average porosity increases by 469% when both the PFE and the update of the inflow boundary condition are taken into account. Nevertheless, as the gypsum dissolves and the fracture enlarges, water finds other, new preferential flow paths, causing a decrease in the gypsum dissolution rate over time. Thus, although the dissolution rate of gypsum is higher than that of other types of rocks in carbonate aquifers, the dominant evolution time of gypsum dissolution may be greater than the human lifetime, as was also previously reported by Campana and Fidelibus [69]. The average gypsum dissolution rate calculated with the numerical model with the PFE is equal to 0.375 mm/year. This rate is in agreement with the gypsum denudation rates reported by Calaforra [110] and Guerrero and Gutiérrez [111] for the Sorbas karst in southeast Spain, which range from 0.28 to 0.42 mm/y. The differences found in the current investigation agree with the findings of Berkowitz [109], who highlighted the importance of modeling guidelines in fractured geological media in correctly interpreting field measurements and laboratory experiments. Therefore, the modeling of fracture propagation is one of the challenges to be addressed in the future, as pointed out by Berre et al. [37]. Table 5. Average porosity in the matrix after 100, 300, and 1000 years and the percentage of the matrix with porosity greater than 0.85 at the end of the simulations of scenarios 1 to 4 of the 2-D reactive transport model. PFE = porosity feedback effect.

Scenario
Permeability The model's results show that the dissolution of highly soluble sulfate minerals, such as bloedite, thenardite, and glauberite, may lead to the creation of additional discontinuities or fractures, thus significantly modifying the groundwater flow in low-permeability media and accelerating dissolution processes in areas near the new fractures. The size of the matrix where the porosity is greater than 0.85 in the stratified model and where the mineral zones with bloedite, thenardite, and glauberite covered extend from the left to the right boundaries (scenario 3) is 3.7% larger than that of the non-stratified model of scenario 1. The increase in porosity in the stratified model in which the mineral zones do not extend from the left to the right boundaries (scenario 4) is only 1%. This result agrees with the findings of Min et al. [19] and Zhao et al. [26], who reported that heterogeneities in mineral distribution can greatly affect the evolution of karst systems. In addition, the reactive transport models confirm that the dissolution of these highly soluble sulfate minerals is responsible for the excess of dissolved sulfate, sodium, and magnesium detected in the hydrogeochemical conceptual model of the Villar de Cañas area [77][78][79]. The models also show that, once these minerals dissolve, high concentrations of dissolved sulfate, sodium, and magnesium can remain in the system for long periods of time due to the low permeability of the medium. This result confirms that dissolution processes in low-permeability evaporitic formations can harm groundwater quality for long periods of time [10].
The results of the sensitivity runs carried out with six different permeability-porosity relationships demonstrated the importance of selecting suitable relationships when simulating karstification processes. The models in which permeability-porosity relationships with large exponents, n (Kozeny-Carman, modified Fair-Hatch, and power law with n = 5), lead to the largest values of hydraulic conductivities, K. The increase in K facilitates the water flow through the fractures and discontinuities near the water inflow boundary and causes the gypsum dissolution front to advance less compared to the models with permeabilityporosity relationships with lower values of n. Thus, the average porosity in the matrix after 100 years of simulation for a power-law relationship with n = 2 is 28.6% larger than that calculated with the Kozeny-Carman relationship (see Table 5). This agrees with the results reported by Lai et al. [34], who reported a similar morphological development of the dissolution front when using the modified Fair-Hatch and Kozeny-Carman relationships, but a relatively different one when using the Verma-Pruess relationship.
CIEMAT performed short-term laboratory leaching tests on gypsum rock samples from the Villar de Cañas site [112]. Preliminary experimental results show patterns similar to those of the numerical model presented here. Longer laboratory tests should be designed in the future to test model predictions in a manner similar to what was done by Sadeghiamirshahidi and Vitton [113].
The work presented here could extended by comparing our model's results with those calculated by using other reactive transport codes with different approaches to simulating karstification processes. Moreover, further research is needed to evaluate the relevance of common simplifications when simulating karst aquifers, such as the laminar flow assumption.

Conclusions
A 2-D groundwater flow and reactive transport model of gypsum karstification in physically and chemically heterogeneous systems has been presented. The model is inspired by the hydrogeochemical conditions of Villar de Cañas in Cuenca (Spain) and consists of a low-permeability matrix composed mainly of gypsum and a discontinuity/fracture that acts as a preferential water pathway. Various scenarios were simulated to numerically assess the long-term evolution of the karst systems and to investigate the relevance of: (1) the dynamic update of flow and transport parameters due to porosity changes provoked by mineral dissolution/precipitation; (2) the spatial distribution of minerals in the rock matrix; (3) the time evolution of the boundary water inflows; (4) the permeability-porosity relationships.
The results of the model simulation allow us to draw the following conclusions:

•
Reactive transport codes and models capable of dynamically updating the flow, transport, and chemical parameters in models with mineral dissolution/precipitation reactions are essential for a realistic representation of karstification processes. The model's results show that the average porosity of the matrix increases from 0.045 to 0.29 after 1000 years of simulation when flow, transport, and chemical parameters and the boundary water inflows are dynamically updated when porosity increases due to gypsum dissolution.

•
The dynamic update of model boundary conditions as porosity and permeability increase due to mineral dissolution leads to the development of karstic formations on shorter time scales. The porosity hardly changes after 1000 years of simulation when the porosity feedback effect is not considered in the numerical model, while the average porosity increases from 0.045 to 0.13 when boundary water inflows through rock discontinuities are taken into account. • Accounting for small quantities of highly soluble sulfate minerals, such as bloedite, thenardite, or glauberite in the numerical model is especially relevant for simulating gypsum karstification. The dissolution of these minerals plays an important role in the development of additional fractures/discontinuities in the rock matrix.

•
The spatial distribution of the highly soluble sulfate minerals in the rock matrix plays a significant role in the long-term evolution of gypsum karstic formations.

•
The dissolution of low-permeability evaporitic formations can harm groundwater quality for long periods of time.

•
Model simulations of karstification processes show significant differences when different permeability-porosity relationships are used. The increase in hydraulic conductivity is largest for the power law with n = 5, Kozeny-Carman, and the modified Fair-Hatch relationships. However, the advance of the gypsum dissolution front is largest for the power law with n = 2 and n = 3 and the Verma-Pruess relationships.
The results of the numerical model also confirm that the dissolution of highly soluble sulfate minerals, such as bloedite, thenardite, or glauberite is responsible for an excess of dissolved sulfate, sodium, and magnesium detected in the groundwater at the Villar de Cañas site and can remain in the system for long periods of time due to the low permeability of the medium. In addition, the present research work demonstrates that the use of numerical models that take the porosity feedback effect into account is essential to reproduce the larger values of porosity and permeability observed along the contacts between geological formations at the Villar de Cañas site.

Supplementary Materials:
The following material is available online at https://www.mdpi.com/ article/10.3390/en15030761/s1, Figure S1: Flow boundary conditions imposed in the numerical model. The water inflow through the right boundary of the matrix was defined as a function of the porosity (Q(φ)), Figure S2: Flow boundary conditions imposed in scenario 2, Figure S3: Spatial distribution of the computed diffusion coefficient after 100, 300, and 1000 years of simulation in scenario 1, Figure S4: Spatial distribution of the computed porosity (top), hydraulic conductivity (middle), and flow velocity (bottom) after 100, 300, and 1000 years of simulation in scenario 1 using the modified Fair-Hatch permeability-porosity relationship. The size of the arrows is not proportional to the magnitude of the flow velocity vector, Figure S5: Spatial distribution of the computed porosity (top), hydraulic conductivity (middle), and flow velocity (bottom) after 100, 300, and 1000 years of simulation in scenario 1 using the Verma-Pruess permeability-porosity relationship. The size of the arrows is not proportional to the magnitude of the flow velocity vector, Figure S6: Spatial distribution of the computed porosity (top), hydraulic conductivity (middle), and flow velocity (bottom) after 100, 300, and 1000 years of simulation in scenario 1 using the power-law permeability-porosity relationship with n equal to 2. The size of the arrows is not proportional to the magnitude of the flow velocity vector, Figure S7: Spatial distribution of the computed porosity (top), hydraulic conductivity (middle), and flow velocity (bottom) after 100, 300, and 1000 years of simulation in scenario 1 using the powerlaw permeability-porosity relationship with n equal to 3. The size of the arrows is not proportional to the magnitude of the flow velocity vector, Figure S8: Spatial distribution of the computed porosity (top), hydraulic conductivity (middle), and flow velocity (bottom) after 100, 300, and 1000 years of simulation in scenario 1 using the power-law permeability-porosity relationship with n equal to 5. The size of the arrows is not proportional to the magnitude of the flow velocity vector.

Funding:
The research leading to these results was funded by ENRESA through a Research Contract with Fundación Universidad de A Coruña. The research was also partly funded by the Spanish Ministry of Economy and Competitiveness (Project PID2019-109544RB-I00) and the Galician Regional Government (Grant number ED431C2021/54 from "Consolidación e estruturación de unidades de investigación competitivas"). The first author had a contract from the FPI Program of the Spanish Ministry of Economy and Competitiveness.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.
Data Availability Statement: Data supporting the reported results can be found at https://osf.io/ ep8gh/ accessed date: 30 December 2021.