Uncertainty quantification for mineral precipitation and dissolution in fractured porous media

In this work we present an uncertainty quantification analysis to determine the influence and importance of some physical parameters in a reactive transport model in fractured porous media. An accurate description of flow and transport in the fractures is key to obtain reliable simulations, however, fractures geometry and physical characteristics pose several challenges from both the modeling and implementation side. We adopt a mixed-dimensional approximation, where fractures and their intersections are represented as objects of lower dimension. To simplify the presentation, we consider only two chemical species: one solute, transported by water, and one precipitate attached to the solid skeleton. A global sensitivity analysis to uncertain input data is performed exploiting the Polynomial Chaos expansion along with spectral projection methods on sparse grids.


Introduction
The Paris agreement, adopted by 196 parties in 2015, aims at limiting global warming to below 2 • C, preferably to 1.5 • C, compared to pre-industrial levels.The reduction of greenhouse gas emissions in the atmosphere is crucial to achieve such long term goal, and requires the transition towards renewable energies, and the subsequent need for effective energy storage; moreover, a safe long term sequestration of CO 2 is considered a promising strategy to reduce emissions into the atmosphere.Many of the aforementioned strategies entail a massive use of the subsurface for fluids injection, storage and production.If, on one hand, it is necessary to guarantee the mechanical integrity of the subsurface to avoid unwanted fracturing and induced seismicity, it is also important to evaluate the effect of chemical reactions on the hydraulic properties of the porous media.Indeed, the injection of water at different temperature and with different solutes concentration with respect to the formation can cause dissolution and precipitation of minerals, with an impact on porosity and permeability.In some cases, as in the case of in situ carbon mineralization, these reactions can be exploited to our advantage to obtain a better storage of CO 2 by 1) chemically trapping carbon by binding it into existing mineral and, at the same time 2) thanks to the change of specific volume associated with the transformation of the minerals, creating a more effective cap-rock, [23].Even if this process occurs over very long time scales in natural conditions, its controlled exploitation is an interesting technological challenge.
A realistic mathematical model of these phenomena encompasses a model for flow in porous media (we focus on single phase flow, assuming small concentrations of gases, if present), coupled with the transport of mobile species and the modeling of both kinetic and equilibrium reactions.Moreover, since reaction rates are usually influenced by temperature, and in view of the possible application of the model to low temperature geothermal plants, the coupled model is completed by the heat equation, [14].
Since fractures are ubiquitous in porous media and have a major impact on flow and transport (both of solutes and heat) this work is focused on the modeling of fractured porous media, where fractures are modeled as lower dimensional objects, which can be much more permeable than the surrounding medium, or nearly impermeable if, for instance, precipitation occurs reducing their aperture, [15].
Many physical, geometrical and geochemical parameters involved in the model are affected by uncertainty, therefore a purely deterministic evaluation of the model is useless without a suitable analysis of the solution variance.In this work we propose an uncertainty quantification workflow based on Polynomial Chaos [30,18] and sparse grids [6] for the sampling of the parameters space: this choice allows us to obtain accuracy with a manageable number of evaluations of the model, which, being coupled and time dependent, has a non-negligible computational cost.The goal is to compute the Sobol indices associated with some input parameters to quantify the impact of the uncertainty of such quantities on a variable of interest, typically the medium porosity.The space distribution of the Sobol indices, or the partial variances, can also give interesting insight into the problem, which, in spite of being reduced to the minimum possible complexity, already exhibits a non-trivial, fully coupled behavior.
The paper is structured as follows.In Section 2 we present the model equation in a homogeneous porous medium; then, in Section 3, we extend the equations to the hybrid dimensional case to account for the coupling with fractures.Section 4 briefly discusses the proposed numerical approximation schemes and some implementation details.Section 5 describes the uncertainty quantification workflow and Section 6 is devoted to the presentation of a a complete set of test cases.Finally, conclusions are drawn in Section 7.

Problem description and mathematical model
We consider the coupled problem of single-phase flow and reactive transport in porous media, accounting for the porosity changes linked to mineral reac-tions.Single phase flow can be assumed, even in the presence of CO 2 , at low gas concentrations, for instance at the edge of the plume, away from injection wells.Typically, a large number of species is needed for a realistic modeling of reactive transport, including solid mineral species and aqueous complexes, and the chemical reactions include equilibrium and kinetic ones.In the following we will discuss the general case and then specialize the model to the simple case of the dissolution/precipitation of a single mineral species.

Single phase flow
Under the assumption of single-phase flow, the fluid pressure p and Darcy velocity q can be computed as the solution of a Darcy problem complemented by suitable boundary conditions prescribing a given pressure or normal flux on the boundary.Here k is the permeability, µ the viscosity, ρ w is the fluid density and g is the gravity acceleration.Note that in (1) we have assumed constant density, but we are accounting for porosity changes through a source term.Indeed, porosity depends on the changes of the mineral volume fraction as detailed in Section 2.3.The intrinsic permeability, which we assume to be isotropic in the bulk porous medium, can be modeled as a function of porosity.We consider the following law where k 0 is the reference value at the initial porosity φ 0 .Note that this dependence will introduce a non-linear coupling among the model equations.

Advection-diffusion-reaction equations for mobile species
Let u i denote the molar concentration of the i-th mobile species, for instance a dissolved mineral or gas.The governing equation reads where d is the diffusion tensor (which may account for mechanical dispersion, and is considered to be the same for all species), q is the Darcy velocity, N r is the number of reactions and ν ir is the stoichiometric coefficient of species i in reaction r.
In particular, we consider a set of N r reactions in the form where X i is a generic species, solid or mobile, ν ir is zero if X i is not involved in the reaction r, negative if X i is a reactant in the forward reaction, and positive if it is a product.The reaction rate R r should be interpreted as the net rate, i.e.
where R + r is the forward reaction rate and R − r is the backward one.Their expressions are often empirical and problem dependent; the one considered in this work is illustrated in Section 2.5.Finally, equation ( 2) should be complemented by initial conditions u i (x, 0) = u i,0 (x), and boundary conditions on the concentration or normal flux.

Mineral species
In the case the species is immobile the evolution equation for its concentration, denoted as w i , reduces to an ordinary differential equation since we do not have the transport and diffusion terms, i.e.
where an initial condition is supplemented as w i (x, 0) = w i,0 (x).However, it may be more convenient to consider a different measure of concentration for solid species, in particular we want to compute the solid volume fractions φ i as the ratio of the volume of mineral i for a unit volume of rock.The volume fractions can be obtained from the concentrations as where η i is the molar volume of the mineral.If we let φ I be the volume fraction of inert minerals (i.e.species that are not affected by chemical reactions) we have that where N s is the number of solid, immobile species.

Heat equation
Since reaction rates usually depend on temperature, we include in our model the heat equation for thermal conduction (based on Fourier's law) and convection in the porous media.We assume thermal equilibrium between rock and water, so that only one equation can be considered instead of two.The temperature field is indicated as θ and its evolution is described by Here c is the effective thermal capacity defined as the porosity-weighted average of the water c w and solid c s specific thermal capacities, where ρ w and ρ s are the densities of the water and solid phase respectively.The effective thermal conductivity Λ is computed in a similar way, as where Λ s and Λ s are the water and solid thermal conductivity.Finally, j models a source or sink of heat in the system.Equation ( 5) is completed by initial conditions, θ(x, 0) = θ 0 (x), and boundary conditions on the temperature or heat flux.

Simplified dissolution-precipitation model
In this work we consider a single, simple kinetic reaction in the form where U and V are the positive and negative ion respectively, and W is the precipitate they can form, [27].The net reaction rate results from the difference between precipitation (forward reaction) and dissolution (backward reaction).Under the hypotheses of electrical equilibrium, we can assume that the concentration of U and V are equal and denoted by u, which from now on is the molar concentration of the mobile species, whereas w is the molar concentration of the precipitate.The reaction rates are modeled as moreover, at equilibrium we have that r d = r p , thus, where u e is the equilibrium solute concentration.The net reaction rate is then Note that precipitation proceeds at a constant rate until w = 0, and is then set to zero. γ

The complete model
The complete model, in the case of our interest, computes (q, p, θ, u, w, φ) by solving the following system of non-linear and coupled equations given by complemented by constitutive laws and suitable initial and boundary conditions.

Hybrid dimensional model for fractured porous media
In the following we introduce the extension of the coupled model (6) to the case of fractured porous media, following [15].For the sake of simplicity we will present the model in the case of a single fracture, geometrically reduced to its centerline, see Figure 1.This model reduction strategy is often adopted in the simulation of fractured porous media to reduce the computational cost by avoiding excessive mesh refinement; moreover, in our case it is particularly convenient since the aperture can change in time due to reactions.For more references on this approach see [24,19,28,26,5,1,3,14] and references therein.
Let Ω be the fractured domain.Following [13], we define γ as a non selfintersecting C 2 curve (if n = 2) or surface (if n = 3).In an equi-dimensional setting a fracture can be defined as the following set of points Thus, we have a subdomain of Ω, separated by the surrounding porous medium Ω by the interfaces γ + and γ − with associated normal vectors n + and n − .We replace the fracture Γ with its centerline γ and we assign a unique normal n γ = n γ,+ to the fracture, see Figure 1.We assume that the fracture is open, with porosity φ γ = 1.

Reduced variables
In the following we denote with the subscript γ the variables in the fracture, and with Ω the variables in the porous medium.Note that, after geometrical reduction, we denote with Ω the domain Ω \ γ.Reduced vectors variables in the fracture are defined as the integral of the tangential components of the corresponding equi-dimensional variables, T (x)q(x, s)ds with T := I − N and N := n γ ⊗ n γ the tangential and normal projection matrices, respectively.The reduced scalar variables are instead the integral average, for each section of the fracture, of the equi-dimensional counterparts w(x, s)ds.
Moreover, we denote by ∇ τ , ∇ τ • the tangential gradient and divergence defined on the tangential space of the fracture.
Following [24] we assume that the permeability k and diffusivity d, can be decomposed in normal and tangential components as

Reduced Darcy flow model
The reduced model for the Darcy flow, which describes the evolution of the reduced Darcy velocity q γ and pressure p γ in the fracture is obtained, following [24], by the following steps: • integration of the mass balance equation in each section of the fracture, to obtain a conservation equation for q γ ; • integration of the tangential component of Darcy law in each section of the fracture to obtain a relationship between q γ and p γ ; • integration, with a suitable approximation, of the normal component of Darcy law in each section of the fracture to obtain two coupling conditions between the fracture and the surrounding medium.
The coupled problem in Ω and γ then reads: 9) where )ds is the reduced source or sink term.Following lubrication theory, the fracture tangential permeability k γ can be expressed as a function of the aperture, as described in more detail in Subsection 3.5.
The conditions on γ ± model the fact that the flux exchange between the fracture and the surrounding porous media is related to the pressure jump via κ γ .Note that, since κ γ , as k γ , can be modeled as a function of γ , if the aperture goes to zero the flux exchange vanishes.

Reduced heat model
The reduced model that describes the evolution of temperature θ is obtained, similarly to the Darcy problem, by integrating the conservation equation in each section of the fracture; the coupling conditions however, which stem from a suitable approximation of the total normal heat flux, should take into account the different nature of the advective and diffusive fluxes in the coupling.The coupled model in Ω and γ reads: where the conservation equation in the fracture accounts for heat flux exchanged with the fracture on both sides, through the terms ψ ± , defined as where θ± is selected as the upwind value, i.e.

Reduced solute and precipitate model
The reduced model that describes the evolution of the solute u γ is similar to the reduced heat equation.By integrating the solute conservation equation in the fracture we obtain the reduced equation, such that Note that the balance equation in the fracture accounts for exchanges with the porous medium, in particular we have that where once again ũ± is the upwind value, i.e.
For the precipitate in the fracture w γ , being the original model an ordinary differential equation valid for each point of the domain, the reduced model in the fracture becomes simply Note that ( 14) is not directly coupled with the corresponding equation in the porous matrix since both describe local phenomena.

Permeability and aperture model
As already mentioned we assume that both components of the permeability k in the fracture follow a law which relates them to the aperture, more precisely where k γ,0 and κ γ,0 are reference coefficients along and across the fracture, respectively, and γ,0 > 0 is the initial aperture.As the porosity changes with mineral precipitation, we consider a similar law to describe the evolution of the fracture aperture γ .We have where η γ represents the molar volume of the mineral as it precipitates on the fracture walls.

Numerical approximation
The numerical schemes for the solution of the deterministic problem (6) are implemented in the PorePy library [20] which provides support for multidimensional coupling, allowing for an easy implementation of the problem in fractured media.

Time integration and splitting strategy
Problem ( 6) is fully coupled in a non-linear way.For the sake of computational efficiency, in this work its solution is based on a non-iterative splitting strategy, with the underlying assumption that the changes to the flow parameters due to chemical reactions are relatively slow.In particular, at each time step we follow the scheme proposed in [15] and: 1. we first solve the Darcy problem to obtain the advective fields q Ω , q γ .The flow problem is discretized in time by approximating the time derivative ∂ t φ by finite differences as where φ * = 2φ n − φ n−1 is the extrapolated value; 2. with q Ω , q γ we solve the heat equation, discretized in time with the Implicit Euler method; 3. then, given the temperature field we solve the advection-diffusion-reaction problem which is in turn split into (a) the advective step, discretized in time with the Implicit Euler scheme, which gives and intermediate solute concentration u * ; (b) the reaction step to compute the final solute and precipitate concentrations (note that the precipitate is not affected by transport).This step is integrated explicitly in time, with the addition of an event location procedure to avoid negative precipitate concentrations.
4. Finally, we update porosity and permeability for the next step.

Space discretization
Space discretization is based on a standard, conforming approach where fractures are honoured by the computational grid and each element of the fractures grid is a face of the porous media grid.However, this assumption could be relaxed allowing for different grid resolutions with the use of mortar variables.Finally, since equations are in mixed-dimensions, all the numerical schemes are applied in different dimensions, i.e. in 2D and 1D.Since the Darcy flux is involved in the advective terms of the transport and heat equations it is of fundamental importance that local mass conservation is fulfilled.Therefore, we approximate the Darcy problem in its mixed formulation and employ a suitable pair of discrete spaces for the pressure and the Darcy flux.In particular we employ the lowest order Raviart-Thomas element pair RT 0 , P 0 .
For the numerical solution of the heat equation and the advection-diffusionreaction equation we apply, in hybrid dimensions, the Finite Volume method.Consistently with the continuous model we consider an upwind approximation of the advective term, whereas the diffusive term is approximated with the two point flux approximation (TPFA), see [10], [11], [9].Since we are considering constrained triangulations to honour the fractures the grid may in principle not be orthogonal.However, we assume that the distortion is small enough to obtain a reliable approximation even with a simple TPFA scheme.

Sensitivity analysis workflow
In this section we present the algorithm employed to approximate stochastic quantities by means of Polynomial Chaos (PC) expansions [30,18].This technique will allow us to compute the sensitivity Sobol indices and to obtain a surrogate model of the problem for a quick evaluation of the quantities of interest.PC expansion have been used to treat a large variety of problems, including elliptic models (see, e.g., [2,25]), fluid mechanics problems [21], and flow and transport in porous media (see, e.g., [17,4]).
The sampling of the uncertain parameters space is performed with pseudospectral projection on sparse grids [6], thus obtaining an accurate estimate with a limited number of evaluation of the deterministic model.This is particularly important since the problem is time dependent and, moreover, the presence of chemical reactions can introduce a fast time scale, constraining the time step amplitude with an increase in the computational cost for each evaluation.

Polynomial Chaos expansion
Let N be the number of parameters ξ = (ξ i ) 1≤i≤N and Ξ the space of possible realizations.For the sake of simplicity the parameters are rescaled so that Ξ = [0, 1] N .Moreover, given a probability measure ρ : Ξ → R + , the inner product of two second-order random variables X(ξ) and Y (ξ) is defined as The Polynomial Chaos (PC) expansion of a variable X(ξ) reads where {X k = X, φ k : k ∈ N N } are the spectral modes of X, the basis functions {φ k (ξ) : k ∈ N N } are multi-variate polynomials chosen to be orthogonal with respect to the product •, • , and the multi-index k = (k 1 , . . ., k N ) denotes the polynomial degree with respect to the parameters ξ i .The PC approximation is then obtained by truncating the expansion in (17) to a finite set K ⊂ N N , which determines the quality of the approximation: The statistical moments of the variables of interest are easily obtained from the PC approximations; e.g., mean, variance, and covariance are given by

Spectral projection method
The coefficients of the PC expansion in (18) can be computed in different ways (see, e.g., [22,7]).This work is based on non-intrusive pseudo-spectral projection, which in our opinion, provides the best trade-off between complexity and precision.The numerical quadrature schemes are constructed as sparse tensorization of a one dimensional formula [16].Then, given M quadrature points and the respective weights {w (q) } 1≤q≤M , the modes (X k ) k∈K are computed as The complexity of the method is governed by the number M of evaluations of the deterministic problem, while the accuracy depends on the PC basis {φ k } k∈K .In order to maximize the accuracy with respect to the computational effort, we adopt a sparse method (cf. Figure 2) hinging on the application of Smolyak's formula [29] directly on the projection operator, rather than on the integration operator.Specifically, the set K is defined as the largest possible one such that the discrete projection is exact for any function spanned by {φ k } k∈K , i.e.: In this work, we have opted for an isotropic sparse tensorization of nested Clenshaw-Curtis quadrature rules, where the level l ∈ N of the sparse grid is the only parameter controlling the quality of the approximation.As l increases, both the number of nodes M and the multi-index set K increase.

Sensitivity analysis
The sensitivity analysis consists in the evaluation of the different contributions of the input parameters on the variance of the solution.This is achieved by the computation of the Sobol indices, defined as the ratio between the partial variance corresponding to the input parameter under investigation ξ i , and the total variance of the quantity of interest X S 1,i := where E(X(ξ)|ξ i ) denotes the conditional expected value of X given ξ i .The indices (S 1,i ) 1≤i≤N in ( 20) are known as principal or first-order Sobol indices and measure the individual contribution of the coefficient ξ i to the variance.Higher-order Sobol indices measure the effect of the concurrent variation of more variables.For instance, second-order Sobol indices read Finally, the total Sobol index is obtained as the sum of the indices involving parameter ξ i S Ti := 1 − Var(E(X(ξ)|ξ \i )) Var(X(ξ)) , where the vector ξ \i = (ξ j =i ) contains all uncertain variables except ξ i .The computation of the Sobol indices results directly from the PC approximation (18) of the variable of interest.Indeed, the partial variances can be explicitly expressed as functions of the spectral modes similarly to the statistical quantities in (19).We refer the reader to [8] for all the details.

Numerical examples
In this section, we present two test cases to validate the proposed approach.In both cases we show the reference numerical solution and discuss the uncertainty quantification related to three parameters affected by uncertainty.
Due to the complexity of the problem, all data in the two cases are the same except for the number of fractures in the network.In the first case, presented in Subsection 6.1, a single fracture touching one boundary is considered, while in Subsection 6.2, ten intersecting fractures are considered.The data for the Darcy problem, the heat equation, and the precipitattiondissolution process are given in Table 1.Furthermore, we assume that the following three parameters are affected by uncertainty and uniformly distributed with mean and variance given by where θ inflow denotes the inflow temperature on the bottom boundary.Being the construction of the sparse grids dependent only on the chosen level and the number of uncertain parameters, the number of simulations needed to construct the PC expansion are: 31 runs for level 2, 111 runs for level 3, 351 for level 4, 1023 for and level 5. Once constructed, the evaluation of the PC expansion takes a small fraction of the time used by the full order model to evaluate the solutions for different times and for different values of the uncertain parameters.Hence, we will use the PC expansion as a surrogate model and evaluate its performances and accuracy properties.In both cases, we will discuss the convergence properties of the PC expansion by increasing the level of the sparse grid considered, the analysis of the Sobol indices, conditioned variances and covariances for selected solutions and finally the probability distribution functions.The simulations are developed with the library PorePy, a simulation tool for fractured and deformable porous media written in Python, see [20].

Single fracture network
Let us consider a domain Ω = (0, 1) 2 with a single immersed fracture defined by the following vertices: (0.1, 0) and (0.9, 0.8).The fracture thus touches the bottom boundary of Ω as depicted in Figure 3. Data and uncertain parameters are reported in the beginning of Section 6.We point out that the fracture is permeable at the beginning of the simulation and due to the solute precipitation its aperture diminishes in time.As a result the effective fracture permeability decreases and until the fracture behaves as a barrier and not any more as a preferential path.
In the following parts, we detail some aspects related to the uncertainty quantification analysis.In Subsection 6.  in Subsection 6.1.2we discuss the variances and covariances of the solutions and in Subsection 6.1.3we introduce the computed probability distribution functions of some of the components of the solutions along the fracture.
The reference solution, corresponding to the average input parameters, is reported in Figure 4, where it is possible to notice the variation of the pressure distribution over time due to the sealing of the fractures and the transport of the solute when the fractures are still highly permeable.

Convergence
In this part, we discuss the convergence and accuracy properties of the surrogate model built with the PC expansion.In fact, the latter can be used to make fast simulations without the need of running the full order model.Since we are dealing with a time dependent problem, we analyse the PC expansion for two different simulation times: after few time steps (t = t 1 = 0.1T ) and at the end of the simulation (t = T ).
Figure 5 presents, for both times, the error decay of the computed solutions  by increasing the level of the sparse grid.For a smooth relation between uncertain parameters and the solutions, we expect exponential decay of the error with respect to the level, which is the behaviour observed in the figure.Additionally, for t = t 1 the error computed is much smaller compared to the end of the simulation, showing a temporal dependence on the quality of the PC expansion.
In Figure 6, we compare the porosity φ computed with the full order model with the one constructed by the PC expansion and the corresponding relative error.The two solutions are in good agreement for both times and the error is rather low.As before, the latter is smaller at the beginning of the simulation and tends to increase at the end, in particular near the inflow boundary at the bottom of Ω.Moreover, at the beginning of the simulation the fracture is highly permeable and, consequently, we observe also a region in the proximity of γ where the error is higher due to stronger geochemical effects.
Finally, Figure 7 compares some of the variables in the fractures computed with the full order model or reconstructed with the PC expansion.Also in this case, the quality of the latter is high and in good agreement with the reference solution.Moreover, in Figure 7 the green dashed lines represent the solutions obtained for each run to construct the PC expansion and can be useful to visualize the variance of the solution.

Analysis of variance and correlations
An important factor is the impact on observed variables of the uncertain input data.In Figure 8 we report the Sobol indices for some variables in the fracture for t = t 1 and t = T .We notice that the three variables are influenced in a similar manner by the uncertain data, and the activation energy is the most important factor, followed by the temperature at the inflow boundary.We notice that while the high temperature front penetrates in the domain and in the fracture, the importance of the activation energy over the temperature inflow   tends to diminish and the latter becomes more important.Note also that the Sobol indices of γ w γ at final time are different from the other two considered variables since η γ dominates the induced uncertainty of the precipitate.
In Figure 9, we compare the variances of the porosity in the domain Ω induced by the uncertain data.Depending on the situation, it might be more convenient to consider the Sobol indices instead of the variances in particular when they span different order of magnitudes.The importance of the activation energy over time is rather interesting.In the beginning the fracture is highly permeable and most of the flow is concentrated around the fracture.However, due to the deposition of new material the fracture becomes very low permeable; thus the water flow tends to avoid it and concentrates more in the right part of the domain, where the solute is transported and becomes precipitate altering the porosity.This effect, less evident, can be seen also in the spatial distribution of the Sobol index of η γ .For the temperature, we still notice the permeability change effect coupled with the inflow of higher temperature that speeds up the precipitation process and, as a result, the porosity decay.
Another important aspect is the interdependency of the output variables, which can be expressed by their covariances.Figure 10 presents this relation between the porosity in the media and several other variables at the two considered times.In the pressure-porosity covariance we observe again the effect of the variation of the fracture behavior in time, fromp permeable to imperme-  able.For the porosity and φ Ω w Ω we observe a negative correlation, namely if more precipitate is deposited in the porous media, less void space is left, and the porosity diminishes.Finally, the temperature front can be seen in the plot of the covariance between φ Ω and θ Ω , i.e. higher values of the latter tend to facilitate the deposition of solute with concentration higher than the equilibrium.This increases the value of the precipitate and consequently lowers the value of the porosity; conversely, in the top part of the domain the effect is the opposite due to the fact that thermal capacity and conductivity depend on the porosity in a complex way.

Probability density functions
We consider now the probability density functions (PDFs) of some variables induced by the uncertain data, which are uniformly distributed.Figure 11 shows, for level 2, the distribution of γ u γ at two points along the fracture, and for both times.We notice that at the beginning of the simulation the PDFs are more spread showing a high variability of the considered variable.However, at the end of the simulation the uncertainty tends to become much smaller and the value is more concentrated.Another important aspect is that the PC expansion outcomes might not fulfill physical bounds, in this case we can get negative values of γ u γ which are not correct.The situation improves by considering a higher level of the sparse grid, indeed as represented in Figure 12 this phenomena is not present any more and the PDFs constructed with the PC expansion are in good agreement with the one computed by the full order model.

Multiple fractures network
We consider now a test case with a network composed of multiple-fractures.The geometry is given by the Benchmark 3 of [12], where fractures at t = 0 are now considered all highly permeable with material properties and problem data equal to the previous test case.A graphical representation of the computational domain is given in Figure 14 on the left.Fractures γ 3 and γ 8 will be considered later for a specific analysis.Some quantities from the reference numerical solution are reported in Figure 13, where it is possible to notice the variation of the pressure distribution over time due to the sealing of the fractures and the transport of the solute when the fractures are still highly permeable.

Convergence
We discuss now the convergence properties of the PC expansion.On the right of Figure 14, we plot the error decay for increasing sparse grid level for multiple variables.Also in this case, the exponential decay expected is confirmed for all the variables.
In Figure 15 we compare the porosity computed by the differential model with the one computed by the PC expansion.On the right we also represent the relative error.The two solutions are in good agreement with a maximum error of 13% confined at the bottom of the domain, the error is much smaller in the other Figure 16 compares some of the variables along the two fractures γ 3 and γ 8 computed by the full order model and by the PC expansion.Also in the fractures, we observe a high quality for the solutions computed with the PC expansion even for γ 8 that has two intersections with other fractures.The jump across the intersection is properly captured, confirming also in this case that the surrogate model from the PC expansion yields a good approximation of the solutions.

Analysis of variance and correlations
In this section, we present and analyze the impact of the uncertainty on some of the computed variables.In particular, Figure 17 shows the Sobol indices for some of the variables of interest in the fractures γ 3 and γ 8 .Since γ 3 is closer to the inflow boundary than γ 8 , the associated Sobol indices behave similarly with respect to the ones of the previous test case.The effect of the inflow is more evident at the lowest tip of the fracture γ 3 with increased impact of θ inflow for γ and γ u γ compared to the other variables.For γ w γ the relation expressed by the Sobol index is less clear.Since fracture γ 8 is more distant form the inflow boundary, the high temperature front at the end of the simulation does not fully reach it.Some effect are still visible since warmer water has been transported by the fractures, especially for γ w γ where the activation energy E becomes less important than η γ compared with the other two variables under investigation.
In Figure 18 we present the partial variances of the porosity with respect to the uncertain data at final simulation time.We notice a small impact of η γ , and a much more significant relevance of E and θ inf low .Note that the effect of these two uncertain parameters on reaction speed is opposite.Moreover, we notice that, on the bottom of the domain, the effect of the temperature is more  pronounced since it is close to the inflow boundary, while for the activation energy E the impact is more predominant away from the inflow and close to the fractures.This can be motivated by the fact that the water and solute get more channelized into the fractures and transported upward.Since the fractures do not touch the outflow boundary, the solute flows again into the rock matrix and then alters the value of the porosity by creating more precipitate.The covariances between some of the computed variables are reported in Figure 19.The correlation between φ Ω and θ Ω is expected: below the warm water front, the increased temperature facilitates the chemical reaction and thus lowers the porosity; above the front, the solute is lower than the equilibrium value and the temperature is lower, hence precipitation may not occur once fractures have been sealed and stopped supplying reactant to the upper part of the domain.The correlation between φ Ω u Ω and φ Ω w Ω is not included here since they are, as expected, equal to -1 in the whole domain.The same applies for γ u γ and γ w γ .

Probability density functions
Finally, we present the PDFs of some of the variables at 25% and 75% of fractures γ 3 and γ 8 for level 2 of the sparse grid.We notice a phenomenon similar to the one observed in the previous test case.For γ 3 we might obtain unphysical values of γ u γ when the PDF is computed by the PC expansion.This issue does  realistic approximation of the physical process, we have considered a mixeddimensional model for the fractures where the latter are represented as object of lower dimension.We have introduced appropriate equations and derived coupling conditions with the surrounding porous media.Additionally, we have applied a splitting strategy to numerically solve the problem by using standard discretization schemes.
In a real scenario, several parameters might be affected by uncertainty.To quantify their effect on the system, we have considered a polynomial chaos expansion constructed by resorting to spectral projection methods on sparse grids.This strategy proved to be very effective, since it provides high quality approximations at low computational cost.Indeed, a limited amount of simulations are needed to construct a surrogate model which can then be used to perform multiple-simulations.Moreover, from the PC expansion one can easily compute useful statistical quantities such as partial variances and Sobol indexes to investigate the impact of input parameters on the model unknowns and gain insight into the complex model couplings.This technique has been applied to two numerical examples by increasing the geometrical complexity of the fracture network.The results obtained showed the validity of the proposed approach.

Figure 2 :
Figure 2: Sparse grid for two parameters of level 2 on the left, level 4 on the centre and level 6 on the right.

Figure 3 :
Figure 3: Domain Ω and fracture γ for the example of Subsection 6.1.

Figure 4 :
Figure 4: Reference solutions with mean value of the uncertain parameters.On the left pressure at time t = 0.1T and on the centre for t = T , on the right the solute for t = 0.1T .Test case of Subsection 6.1.

Figure 5 :
Figure 5: Error convergence for increasing level of the sparse grids, on the left at time t = 0.1T and on the right at final time t = T .Test case of Subsection 6.1.1.

Figure 6 :
Figure 6: On the left porosity in the media computed with the original model and on the centre with the polynomial chaos expansion, on the right the error between them.On the top at time t = 0.1T and on the bottom at final time t = T .Test case of Subsection 6.1.1.

Figure 7 :
Figure 7: Solutions along γ: in blue computed with the original model, in red with the polynomial chaos expansion and in green computed by the original model on each sparse grid node (sparse grid level 2).On the top at time t = 0.1T and on the bottom at final time t = T .Test case of Subsection 6.1.1.

Figure 8 :
Figure 8: First order Sobol indexes for different unknowns along γ, on the top at time t = t 1 and on the bottom at final time t = T .The local coordinate of the fracture starts from the bottom boundary.Test case of Subsection 6.1.2.

Figure 9 :
Figure 9: On the left the variance of porosity conditioned to the activation energy E, on the centre conditioned with η γ expansion, and on the right conditioned with the temperature inflow θ inf low .On the top at time t = 0.1T and on the bottom at final time t = T .Test case of Subsection 6.1.2.

Figure 10 :
Figure 10: Covariances between different solutions in the porous media.On the top at time t = 0.1T and on the bottom at final time t = T .Test case of Subsection 6.1.2.

×10 2 ×10 2 ×10 3 ×10 3 Figure 11 :
Figure 11: Probability distribution function of γ u γ for level 2 at two points in γ, on the top at time t = 0.1T and on the bottom at final time t = T .Test case of Subsection 6.1.3.

Figure 12 :
Figure 12: Probability distribution function of γ u γ for level 5 at two points in γ, on the top at time t = 0.1T and on the bottom at final time t = T .Test case of Subsection 6.1.3.

Figure 13 :
Figure 13: Reference solutions with mean value of the uncertain parameters.On the left pressure at time t = 0.1T and on the centre for t = T , on the right the solute for t = 0.1T .Test case of Subsection 6.2.

Figure 14 :
Figure 14: On the right the computational domain with boundary conditions and on the left the error convergence for increasing level of the sparse grids.Test case of Subsection 6.2.1.

Figure 15 :
Figure 15: On the left porosity in the media computed with the original model and on the centre with the polynomial chaos expansion, on the right the error between them.Test case of Subsection 6.2.1.

Figure 16 :Figure 17 :
Figure 16: Solutions along two γ i and in blue computed with the original model, in red with the polynomial chaos expansion and in green computed by the original model on each sparse grid node.On the top γ 3 and on the bottom for γ 8 .Level considered 2. Test case of Subsection 6.2.1.

Figure 18 :
Figure 18: On the left the variance of porosity conditioned to the activation energy E, on the centre conditioned with η γ expansion, and on the right conditioned with the temperature inflow θ inf low .Test case of Subsection 6.2.2.

Figure 19 :
Figure 19: Covariances between different solutions in the porous media.Test case of Subsection 6.2.2.