Upscaling of Non-isothermal Reactive Porous Media Flow with Changing Porosity

Motivated by rock–fluid interactions occurring in a geothermal reservoir, we present a two-dimensional pore scale model of a periodic porous medium consisting of void space and grains, with fluid flow through the void space. The ions in the fluid are allowed to precipitate onto the grains, while minerals in the grains are allowed to dissolve into the fluid, and we take into account the possible change in pore geometry that these two processes cause, resulting in a problem with a free boundary at the pore scale. We include temperature dependence and possible effects of the temperature both in fluid properties and in the mineral precipitation and dissolution reactions. For the pore scale model equations, we perform a formal homogenization procedure to obtain upscaled equations. A pore scale model consisting of circular grains is presented as a special case of the porous medium.


Introduction
Geochemical reactions can affect the permeability in a geothermal reservoir. As the injected water and the in situ brine have different temperatures and chemical composition, reservoir rock properties can develop dynamically with time as the fluids flow through the reservoir. Mineral dissolving and precipitating onto the reservoir matrix can change the porosity and hence the permeability of the system. As mineral solubility is affected by the cooling of the rock and by the different ion content in the saturating fluids, large changes in permeability can occur. The interaction between altering temperature, solute transport with mineral dissolution and precipitation and fluid flow is highly coupled and challenging to model appropriately as the relevant physical processes jointly affect each other (Bringedal et al. 2014). The effect of changing porosity and hence permeability through the production period may have severe impact on operating conditions. The ion content of the injected cold water is normally different than the original groundwater, affecting the equilibrium state of the system. Also, as the solubility of minerals generally depends on temperature, the cooling itself will affect the equilibrium state. As reported from field studies and simulations, porosity and permeability changes due to precipitation and dissolution of minerals as silica, quartz, anhydrite, gypsum and calcite can be observed (McLin et al. 2006;Mroczek et al. 2000;Pape et al. 2005;Sonnenthal et al. 2005;Wagner et al. 2005;White and Mroczek 1998). Modeling of the mineral precipitation and dissolution is important in order to understand the processes and to better estimate to which extent the chemical reactions can affect the permeability of the porous medium.
When investigating porosity and permeability changes, it is important to understand the underlying processes at the pore scale. The pore geometry affects the reaction rates as the reactive surface is changed, and the permeability is depending on how the geometry changes. To achieve equations that quantify how the reaction rates and permeability depend on the pore scale effects, we start with a pore scale model and derive the Darcy scale model by homogenization. Pore scale models incorporating mineral precipitation and dissolution have been studied earlier in (van Duijn and Pop 2004;van Noorden et al. 2007), and the corresponding Darcy scale models have been investigated further in (van Duijn and Knabner 1997;Knabner et al. 1995). These papers assume that the pore geometry is not changed by the chemical reactions, which is a valid assumption when the deposited or dissolved mineral layer is thin compared to the pore aperture. Investigations honoring the porosity changes may be found in (Kumar et al. 2011;van Noorden 2009a, b), where mineral precipitation and dissolution have been considered on either circular grains or in a thin strip. In these papers, the position of the interface between grain and void space is tracked, giving a problem with a free boundary. Similar models can also be obtained for biofilm growth (van Noorden et al. 2010), for drug release from collagen matrices (Ray et al. 2013), and on an evolving microstructure (Peter 2009). These models do not include any temperature dependence.
The present work builds on (Bringedal et al. 2015), where mineral precipitation and dissolution are considered in a thin strip. There, the freely moving interface between the mineral layer and the void space is modeled as an unknown function depending on time and the location along the pore axis. Here we extend the results to a more general geometry with many solid grains periodically distributed in the medium. The geometry is as considered by van Noorden (2009a), but as we include temperature effects in fluid flow and chemical reactions, the resulting model is much more complex. Including temperature effects requires two new model equations at the pore scale: energy conservation in void space and in the grain space. The coupling between these two equations raises several issues in the upscaling process that are to be elaborated on. For geothermal systems, the temperature dependence is essential. We assume there is no phase change, but mention that Duval et al. (2004) upscale two-phase flow with phase change using volume averaging. Recently, Radu et al. (2013) proposed a mass conservative mixed finite element scheme for reactive transport with varying porosity for saturated/unsaturated porous media. Other applications of non-isothermal reactive transport models can be combustion models, and we refer to Lu and Yortsos (2005) who consider a pore-network model including filtration combustion.
The structure of this paper is as follows. In Sect. 2, we describe the pore scale model, while in Sect. 3, we perform formal homogenization on the model equations and obtain upscaled equations. The paper ends with a summary with some comments on applications and a presentation of a special case of the model equations in Sect. 4, together with some concluding remarks.

Pore Scale Model
We consider a two-dimensional domain Ω with boundary Γ . The domain is inside a square x = (x 1 , x 2 ) ∈ (0, L) 2 for some positive number L. The domain is perforated, consisting of a connected pore space Ω (t) and grain space G (t), and the boundary between them Γ (t). As we have two length scales, the domain size L and the typical pore size l, we use the superscript = l/L to emphasize dependency on both length scales. The grain space consists of some non-reactive solid surrounded by minerals. These minerals are the result of a precipitation process and may dissolve. The pore and grain space are given in a periodic manner, where the grains are assumed not to touch each other, see Fig. 1. The grains are assumed not to touch each other to simplify the upscaling as less precautionary steps are then needed, but also to avoid clogging. Figure 2 shows the zoomed-in pore structure. We consider a square Y having local variables y = (y 1 , y 2 ) ∈ (0, l) 2 , where l is much smaller than L. The region Y consists of grain space G i j (t), void space Ω i j (t) and the moving boundary between them Γ i j (t), where the outward unit normal of the void space is denoted n . The indices i, j denote which subdomain in the total domain Ω we consider, and the set of indices {i, j} is such that the whole domain is covered. As the subdomains Y are given in a periodic manner, we have periodicity in the variable y. The microscopic variable y and macroscopic variable x are related through x = ( i, j) + y.
We use the level set function S (x, t) to describe the position of the interface Γ i j (t) and is given as This choice of S (x, t) assures that the gradient of S points into the grain space and is parallel to n , hence where the notation | · | means the length when applied to a vector and the size when applied to a set. The level set function tracks the position of the interface, hence also the size and geometry of void and grain space as time evolves. As will result from below, the evolution of S depends on unknowns of the model, leading to a porous medium where the size and geometry of the void space are in evolution. In particular, this evolution leads to a non-periodic structure, although the grains are assumed periodically distributed. The boundary Γ i j (t) moves with some normal velocity v n due to the mineral precipitation and dissolution. The Rankine-Hugoniot condition guarantees conservation of quantities across a moving boundary (Fasano 2005): where u is the preserved quantity (e.g., mass or energy) and j is the flux of this quantity. The use of square brackets means the jump of the quantities and is the difference between the quantities at each side of the interface. The condition states that the net loss (or gain) rate of the preserved quantity u equals the difference in fluxes transporting u across the interface (Fasano 2005).
We apply the level set equation together with conservation of ions, mass, momentum and energy to form a complete set of equations at the pore scale and refer readers to, e.g., Patankar (1980) for justification of the conservation equations and to (Osher and Fedkiw 2003) for a detailed treatment of the level set equation. We prescribe boundary conditions at the internal boundary Γ (t). For a complete model to be used for computer simulations, boundary conditions at the external boundary Γ are required, as well as initial conditions, but as these are not necessary for the upscaling process, they will not be specified. The model equations are briefly presented in this paper, and we refer to (Bringedal et al. 2015) for more details on the model.

Interface Evolution and the Level Set Equation
As discussed, the medium is covered by translations of a square as presented in Fig. 2, with periodicity across the external boundary ∂Y . The solid part G i j (t) consists of a non-reactive part denoted B i j located in the center of the square and is surrounded by a mineral part that can dissolve and precipitate. A mineral molecule consists of two ions u 1 and u 2 . We assume for simplicity that u 1 and u 2 have the same initial and boundary conditions and will hence be equal. We denote their concentration with u in the following. The interface Γ i j (t) moves as mineral molecules dissolve or precipitate, and the normal velocity is proportional to the difference between the reaction rates (Knabner et al. 1995;van Noorden 2009b): where ρ C is the molar density of the mineral and is assumed constant, and f p and f d are the precipitation and dissolution rates of the chemical reaction, respectively. The rates are (Chang and Goldsby 2014;van Noorden 2009a) f where k 0 is a positive rate constant, E is the activation energy, R is the gas constant, T f is the fluid temperature, and K m (T f ) is the equilibrium constant for the mineral. Dissolution takes place as long as there are minerals present: The union of the non-reactive solid parts is denoted B = ∪B i j . Hence, dist(x, B ) is the distance from the present point to the non-reactive solid and can be thought of as the thickness of the mineral part. The function w(dist(x, B ), T f , u) ensures that the dissolution rate cannot exceed the precipitation rate when there are no minerals left, and is given by We will denote f = f p − f d as the net rate for increasing mineral width. This particular choice of reaction rates is inspired from (van Duijn and Pop 2004;Knabner et al. 1995;van Noorden 2009a, b), where the isothermal case is considered. Other rates can be adapted straightforwardly. Note that the reaction kinetics are assumed to be slow, as fast kinetics will require other considerations in the upscaling process; see, e.g., Kechagia et al. (2002) which shows a case where upscaling brakes down. We refer to (Mikelic et al. 2006) for upscaling of reactive flow through a thin strip for dominant Damköhler number.
Combining the above equation with (2), we get Note that the level set equation and the level set function are defined in the entire domain Ω. This is handled by defining a continuous extension of the reaction rates to Ω.

Conservation of Ions
There are two active ions in the fluid, both having molar concentration u. They satisfy the convection-diffusion in the void space: In the above equation, D is the diffusion coefficient which we assume to be constant, and q is the fluid velocity. The Rankine-Hugoniot condition (1) for conserving ions across the moving interface is As the minerals have zero flux, only the ion flux appears on the left-hand side. The difference on the right-hand side is the jump of the conserved quantity u. A mineral molecule consists of one of each type of the two ions; hence, the difference (ρ C − u) represents the jump across the interface for each of the ions.

Conservation of Mass
As the fluid consists mainly of water, the fluid molar density ρ f is assumed to not be affected by the chemical reactions, but depends on temperature. Hence, the mass conservation equation As we include the temperature effects on fluid density, density differences are accounted for and we cannot use the incompressible mass conservation equation as in van Noorden (2009a). At the boundary, ions can leave the fluid and become part of the grain space instead. The Rankine-Hugoniot boundary condition applied to mass is The difference on the right-hand side is the jump of the mass across the moving boundary. As a mineral molecule consists of the mass from two ions, the molar density ρ C is multiplied with the factor 2.

Conservation of Momentum
We assume the fluid is Newtonian, that the stress tensor is a linear function of the strain rates, that the fluid is isotropic, and that the body forces are such that the fluid is at rest at hydrostatic pressure. Hence, where p is pressure and μ is viscosity of the fluid. No-penetration conditions are assumed at the boundary, meaning that q has no tangential component at the interface, but is parallel to the normal vector n . Combining with Eq. (8), the new boundary condition becomes

Conservation of Energy
We distinguish between two temperatures, the fluid temperature T f and grain temperature T g . We assume no viscous dissipation; hence, energy transfer in the fluid phase can happen through diffusion and convection: In the grain space, flow is not possible, hence In the above equations, c f and c g are specific heats, and k f and k g are heat conductivities, of fluid and mineral, respectively, and are all assumed constant. The Rankine-Hugoniot condition for conservation of energy across the moving interface is and we also assume temperature continuity at the interface: The temperature continuity condition comes from the assumption of local thermodynamic equilibrium, and this second boundary condition at Γ (t) is needed to link the two energy conservation equations properly.

Non-dimensional Model Equations
To achieve non-dimensional quantities, we introduce t ref , Non-dimensional variables and quantities are denoted with a hat and are defined aŝ We emphasize dependence on the small variable by denoting our main variables with as a superscript. Observe that the dimensionless viscosity scales with 2 , which is a natural assumption leading to non-trivial upscaled flows. Since we will only use non-dimensional variables in the following, we skip the hat. Due to the scaling of y, the local square Y is now a unit square.
The non-dimensional level set Eq. (4) becomes The non-dimensional reaction rates are The Damköhler number k is assumed to be independent of , which means that the reaction rates are at the same timescale as the convective timescale. This assumption is applicable for most rock-fluid interactions as this branch of chemical reactions is known to be slow (Bundschuh and Zilberbrand 2011). Including fast kinetics would require extra care, as discussed by Kechagia et al. (2002), Mikelic et al. (2006).
The convection-diffusion Eq. (5) becomes with the boundary condition (6) now written as when (2) is inserted. Note that an underlying assumption is that the dimensionless diffusion coefficient D is not depending on ; hence, diffusion and convection occur at the same timescale, which is a valid assumption in the context of geothermal energy production when the injection rate is not too large. The mass conservation Eq. (7) transforms into The corresponding Rankine-Hugoniot boundary Eq. (8) has the dimensionless form The momentum balance Eq. (9) becomes while the boundary condition (10) is The non-dimensional form of the energy conservation Eqs. (11) and (12) is and These three constants are all assumed to not depend on , but are of order 1. Hence, heat diffusion in grain and fluid space is at the same timescale as the convective timescale, which is a valid assumption when the injection rate is not too large. The boundary condition (13) is written and the continuity condition (14) is The fluid density and viscosity are assumed to depend linearly on the fluid temperature T f , hence and for some positive constants β ρ f and β μ . The assumption of linear relationship is a common simplification, but using other continuous relationships between density/viscosity and temperature is straightforward.

Asymptotic Expansions
We incorporate explicit dependence of the microscopic variable y = (y 1 , y 2 ) into the main variables by introducing the asymptotic expansions and similarly for the other unknowns. Due to periodicity, all the functions S i are periodic in y. The separation between x and y enables to capture both slow, macroscopic variations (variations in x) and fast, microscopic variations (variations in y).

Preliminaries
Due to the relation between x and y, the gradient of a function Since n = ∇S /|∇ S |, we assume n = n 0 + n 1 + O( 2 ) and find the expressions for n 0 and n 1 by inserting the expansion for S : Above we used the Taylor expansion for the expression in the square root. This means, By introducing τ 0 , the unit vector that is orthogonal to n 0 , we can write n 1 as In the upscaling process, we need to be able to formulate boundary conditions at the zero level of the moving boundary; that is, for those values of y such that S 0 (x, y, t) = 0. The boundary conditions in the previous section are formulated at Γ (t), i.e., at every x where S (x, t) = 0. To overcome this challenge, we need to assume the existence of a parameterization of Γ i j (t) and that S 0 (x, y, t) and S 1 (x, y, t) have Taylor expansions, which also means assuming the boundary to have some regularity. These assumptions are necessary to perform the upscaling procedure, but excludes some pore geometries. At the boundary Γ i j (t), the level set function S is zero, so we assume there exists a parameterization k (s, t) such that S (k (s, t), t) = 0.
We further assume that k (s, t) has the asymptotic expansion where x = ( i, j). Using the asymptotic expansion for S , the periodicity of S i in y, and the Taylor expansions for S 0 and S 1 around (x, y) = (x, k 0 ), we obtain 0 = S (k (s, t), y) where D 2 ξ is the second-order differentiation matrix with respect to x and y and all derivatives are evaluated in (x, k 0 ). Collecting lowest order terms yields S 0 (x, k 0 , t) = 0, meaning that k 0 parameterizes locally the zero level set of S 0 . Second lowest order terms are We search for a k 1 that is aligned with n 0 ; hence, k 1 = λ 1 n 0 , where λ 1 = − 1 |∇ y S 0 | (S 1 + k 0 · ∇ x S 0 ). Third lowest order terms are We convert the boundary conditions given for x ∈ Γ (t) such that these are given at y ∈ Γ 0 (x, t) where Γ 0 (x, t) = {y | S 0 (x, y, t) = 0}, and which is parameterized by y = k 0 (s, t). Assuming the boundary condition is given as we use that x = k (s, t) at Γ (t) and then expand K (x, y, t) around (x, k 0 ): where we have inserted k 1 and k 2 when applicable. We use the two following results proved in (van Noorden 2009a).
Lemma 1 (Lemma 3.1 van Noorden 2009a) Let g(x, y, t) be a scalar function such that g(x, y, t) = 0 for y ∈ Γ 0 (x, t), then it holds that

Lemma 2 (Lemma 3.2 van Noorden 2009a) Let F(x, y, t) be a vector-valued function such
A result similar to Lemma 2 will be used for the energy equations. As only model equations in the void space are included in (van Noorden 2009a), the following result, considering coupled model equations in the void space and grain space, is developed and proved: Proof The proof is quite similar to the proof of Lemma 3.2 in (van Noorden 2009a), but requires some extra care on how the grain space is treated. We split S 1 into its positive and negative parts and focus on the positive part denoted for some nonnegative number δ. We use opposite sign for the S 1 term to assure that G δ , t), are not the same. These boundaries are, however, equal for δ = 0 as they both reduce to Γ 0 (x, t). We parameterize the two boundaries with k + Y (s, δ) and k + G (s, δ) such that t) = 0 for 0 ≤ s ≤ 1. By differentiating these two equations with respect to δ using the chain rule and product rule, and then evaluating the resulting equations at δ = 0, we obtain The unit normal vectors of Γ δ+ Y (x, t) and Γ δ+ G (x, t) are denoted n δ+ Y and n δ+ G . For convenience, both are oriented to be pointing into the grain space; hence, both are equal to n 0 when δ = 0. The integrals of ∇ y · F Y and ∇ y · F G over Y δ + (x, t) and G δ + (x, t), respectively, are zero, and the derivatives of the integrals with respect to δ are also zero. Hence, We now evaluate each of the terms in both integrals at δ = 0. For the last terms in the integrals, we use that ∂ δ (|∂ s k + Y |)| δ=0 = −∂ δ (|∂ s k + G |)| δ=0 since the variations in δ were made with opposite sign. Since n δ+ Y | δ=0 = n δ+ G | δ=0 = n 0 , we get that the last terms in the two integrals cancel out as n 0 · (F Y − F G ) = 0 at Γ 0 (x, t). For the middle terms in the two integrals, we use that n 0 · ∂ δ k + Y | δ=0 = −[S 1 ] + /|∇ y S 0 | and n 0 · ∂ δ k + G | δ=0 = [S 1 ] + /|∇ y S 0 |. Finally, calculating the derivatives of the unit normals with respect to δ, gives that ∂ δ n δ+ |∇ y S 0 | τ 0 . Inserting this into the above equation results in Repeating the same steps for [S 1 ] − and subtracting the two parts complete the proof.
Note that even though we consider a two-dimensional pore scale model, the model equations defined in Sect. 2 could easily be defined in a three-dimensional setting. However, the previous lemmas depend on the tangent vector of the interface Γ 0 (x, t), which would need to reformulated into a tangent plane. Also, the unit tangent vector τ 0 would not be uniquely defined in a three-dimensional setting.

The Level Set Equation
For the asymptotic expansion, we need Lipschitz continuous reaction rates. Therefore, we replace (3) by its Lipschitz approximation: for some small δ > 0. As proved in (Kumar et al. 2013), in the limit as δ approaches zero, one ends up with the original (3). Thus, by f δ (T f , u , y), we mean the reaction rate f where w (dist(x, B ), T f , u ) is replaced with w δ (dist( y, B )). The level set Eq. (15) becomes We collect the lowest order terms and then let δ approach zero, giving back the original w(dist( y, B ), T f , u ); leading to

Conservation of Ions
We insert the asymptotic expansions for u and q into (17) and into the boundary condition (18) and apply Eq. (29) to obtain a boundary condition valid on Γ 0 (t). The lowest order terms are Recall that u 0 is periodic in y. This means u 0 cannot depend on y; hence, u 0 = u 0 (x, t). The second lowest order terms are where we have introduced the notation F u = D∇ x u 0 + D∇ y u 1 − q 0 u 0 . Finally, the third lowest order terms are We integrate (33) over Y 0 (x, t) and apply the boundary condition (34): The reactive term f (T f 0 , u 0 ) can be rewritten using the level set Eq. (30) and its relation to the normal velocity. Since G 0 (x, t) has volume (1 − |Y 0 (x, t)|), we can rewrite the reactive term using Reynolds transport theorem such that The last six integrals add up to zero; We are therefore left with

Conservation of Energy
Upscaling the two energy conservations (23) and (24) follows similar steps as in the previous section, but requires some extra care as the grain space plays an explicit role. Following (27), the leading order term in the fluid density reads Inserting the relevant asymptotic expansions into (23) and (24) and into the internal boundary conditions (25) and (26) while applying (29), the lowest order terms will be As T f 0 and T g0 are periodic in y, the only solution to this problem is that T f 0 and T g0 do not depend on y, which together with the continuity condition implies T f 0 = T g0 = T 0 (x, t). The second lowest order terms are where Finally, the third lowest order terms are We integrate (41) over Y 0 (x, t) and (42) over G 0 (x, t) and sum the integrals. In the integrals involving a divergence with respect to x, we interchange the order of integration and differentiation. In the integrals involving a divergence with respect to y, we apply Gauss' theorem. We rewrite boundary terms using (43) and insert the expressions for n 1 and λ 1 where applicable. We then get As before, the last six integrals add up to zero; The reactive term can be rewritten using the level set equation as earlier. We therefore end up with where ρ f 0 is given by (36).

Conservation of Mass
While van Noorden (2009a) could incorporate the incompressible mass conservation equation into the upscaling of the momentum equation, the compressible mass conservation Eq. (19) must be explicitly accounted for in the upscaling process due to the varying density. We insert the asymptotic expansions into (19) and (20) and apply (29) as before. To show the similarities with conservation of ions and energy, we introduce the notation F m = −ρ f 0 q 0 . The lowest order terms are Since ρ f 0 does not depend on y, the first equation can also be read as ∇ y · q 0 = 0 in Y 0 (x, t).
The second lowest order terms are We integrate (48) over Y 0 (x, t), replacing the order of integration and differentiation where the integrand involves ∇ x and apply Gauss' theorem when the integrand involves ∇ y . We apply (49) and insert expressions for n 1 and λ 1 at the same time, obtaining As earlier, we rewrite the reactive term using the level set equation and Reynolds transport theorem. The six remaining boundary integrals add up to zero; We are therefore left with

Conservation of Momentum
We insert the necessary asymptotic expansions into (21) and (22) and apply (29). Following (28), the leading order term in the fluid viscosity is We collect the lowest order terms from the model equation and the boundary condition, obtaining As ∇ y p 0 = 0 in Y 0 (x, t), we conclude that p 0 = p 0 (x, t). The second lowest order terms from the model equation are It is known from (46) that ∇ y · q 0 = 0 in Y 0 (x, t); hence, the last term disappears. Since μ f 0 is independent of y and ∇ y · (∇ y q 0 ) T = ∇ y (∇ y · q 0 ), we can rewrite the above equation into with the fluid viscosity given by (50).

Cell Problems
Here we eliminate u 1 from (35). We make use of the second lowest order Eq. (31) and boundary condition (32). Since this problem is linear, u 1 can be written as a linear combination of the derivatives of u 0 with respect to x, for unknown functions v 1 (y) and v 2 (y). Rewriting (31), using ∇ y · q 0 = 0 and inserting the above expansion for u 1 yield Doing the same for the boundary condition (32), recalling that n 0 · q 0 = 0 from (47), we get From this equation, we obtain n 0 · (e j + ∇ y v j (y)) = 0 at Γ 0 (x, t) for j = 1, 2.
The two functions v 1 (y) and v 2 (y) are determined by Eqs. (53) and (54) together with requiring periodicity in y. This model problem is called the cell problem for u. Inserting the series expression for u 1 into the model problem (35) results in for unknown functions Π j (y) and vector functions ω j (y). Inserting the expressions for p 1 and q 0 into (52) yields The boundary condition (51) gives us that Through ω j (y), we get that the integral of the lowest order velocity can be written as where K = {k i j } is the matrix given by where ω j i (y) is the i'th component of the vector function ω j (y).

Summary and Discussion
We summarize the derived effective equations, without making any assumptions on the grain geometry and hence the shape of the level set function. Since we only use the lowest order expansions, we skip the subscript 0 from the variables. We have five unknowns: S (x, y, t), u(x, t), T (x, t),q(x, t) and p(x, t). All but the first depend only on spatial variable x and are defined for all x ∈ Ω, while the level set function S(x, y, t) also depends on the microscopic variable y ∈ [0, 1] 2 . The five upscaled equations are (x ∈ Ω, t > 0) The matrices A u , A f , A g and K have components given by while the functions v j (y), θ j (y) and ω j (y) are solutions to the cell problems (x ∈ Ω, t > 0) together with the periodicity in y for v j (y), θ j (y), Π j (y) and ω j (y), j = 1, 2.
We observe that the three upscaled conservation equations for ions, mass and energy show similarities in the first term; the time derivative consists of a volume-weighted sum of the preserved quantity in void and grain space. Hence, we take the time derivative of the total quantity occurring in the unit square. The upscaled diffusive terms include a matrix with components arising from cell problems, instead of a volume factor as in the case of a thin strip (Bringedal et al. 2015). This way we take into account the geometry of the pore space and not only the size. Same happens for Darcy's law where the permeability tensor K depends on pore geometry through the cell problems. While the thin strip case in (Bringedal et al. 2015) reduced to a one-dimensional model, the present model remains two-dimensional also in the upscaled version, which is honored through the diffusion and permeability tensors.
Compared to the reactive transport model TOUGHREACT (Xu et al. 2012), which is well used for geothermal reservoirs, the upscaled model equations found here are similar to the model equations that are applied in TOUGHREACT. The conservation equations for mass, ions and energy have the same structure with porosity-weighted averages in the time derivative, convective term with an average volume flux, and a porosity-weighted diffusive term that includes heat transfer in the grains. Hence, the upscaled model we have found behaves as expected compared to well-known models, but introduces also the pore scale effects: While TOUGHREACT (Xu et al. 2012) uses a simple equation for the porosity evolution due to mineral precipitation and dissolution, we include an explicit equation through the level set equation. Further, TOUGHREACT (Xu et al. 2012) uses only scalar diffusion and permeability, while we find that the permeability and diffusion will in general be tensors. TOUGHREACT (Xu et al. 2012) uses similar reaction rates for kinetic mineral precipitation and dissolution as considered here, but relies on estimates on the reactive surface, which in our upscaled model is handled explicitly through the level set formulation. The upscaled model contains a discontinuous reaction rate, and Agosti et al. (2015a, b) show well posedness for a similar model and suggest a numerical strategy for the implementation of such a reaction rate.
The porosity is not explicitly used in the upscaled system of equations, but the void and grain space volumes are present and are defined through the level set function. The derivatives of the void and grain space volumes are handled through the level set equation which connects these derivatives with the reaction rates. The upscaled equations are computationally cheaper than the original model problem as the two scales x and y are separated, and the microscopic variable y only appears through the cell problems and the level set equation. Redeker and Eck (2013) proposed an adaptive solution strategy for a model with two-scale dependence through a similar separation of the micro-and macroscale as considered here, and they showed that the scheme was convergent. Note that the upscaled model is not capable of describing cases with clogging as we assumed the presence of a connected pore space in Sect. 2. Hence, in numerical simulations, caution must be made to avoid situations where clogging could occur and abort the simulation if clogging were to take place.
Compared to the work of van Noorden (2009a) who considers the isothermal case, the present model includes compressible flow due to temperature-dependent fluid density and honors heat transfer in both fluid and grain space, which is expressed through a single energy conservation equation in the upscaled model due to the assumption of local thermal equilibrium at the pore scale. Temperature effects are also found in the varying fluid viscosity in Darcy's law. We can note that at isothermal conditions, the present model reduces to the one by van Noorden (2009a).
The upscaled system of equations introduces some useful results in how the permeability and diffusion tensors of the coupled system evolve due to the chemical reactions that can be incorporated into simulator codes. TOUGHREACT (Xu et al. 2012) incorporates permeability changes based on the work by Verma and Pruess (1988) by utilizing a simple power law relationship between the scalar permeability and the porosity. Only scalar solute diffusion and heat conduction are considered in TOUGHREACT, but as shown by the present upscaling process, these will in general be tensors. The upscaled system of equations formulated here includes permeability and diffusion tensors based on cell problems, giving a more detailed and accurate description of the transport processes.
Using mass conservative mixed finite element methods (MFEM), Frank (2013) implemented a Stokes-Nernst-Planck-Poisson system based on upscaled pore scale equations, taking into account the permeability and diffusion tensors depending on pore scale effects through cell problems. Frank (2013) considered several grain shapes in the fixed geometry case and applied circular grains for varying porosity. Numerical computations using MFEM on circular grains with varying radii were also considered by Ray et al. (2013) for drug release from collagen matrices. Ray et al. (2013) derived approximate solutions of the corresponding cell problems and found good correspondence with experimental results. van Noorden (2009a) implemented his upscaled pore scale equations for reactive transport on circular grains with altering grain radius, fitting the effective diffusion and permeability with a rational polynomial function. Both Frank (2013) and van Noorden (2009a) consider incompressible flow, but extending to compressible flow should be straightforward by replacing the mass conservation equation. As the upscaled temperature equation has the same structure as the upscaled solute transport equation, it can be implemented using similar steps as performed by Frank (2013), Ray et al. (2013) or van Noorden (2009a. If we make assumptions on the pore geometry and hence the shape of the level set function, our model equations can be further simplified. For special choices of the level set function, e.g., circular geometry, we obtain a system of equations only depending on the macroscopic variable (van Noorden 2009a). Specifically, assuming that the grains are circular and centered in the middle of the unit square, we can use the level set function S(x, y, t) = R 2 (x, t) − (y 1 − 1/2) 2 − (y 2 − 1/2) 2 , where R(x, t) is the radius of the grain. Using this level set function and pore geometry, we reformulate the above equations in terms of R(x, t). All the five unknowns R(x, t), u(x, t), T (x, t),q(x, t) and p(x, t) are now defined for x ∈ Ω and do not depend on the microscopic variable y. The equations read now (x ∈ Ω, t > 0): The reaction rate f uses the distance between R and R min , where R min is the radius of the non-reactive part, to calculate the width of the mineral layer. The matrices A u , A f , A g and K depend only on R(x, t) as the integration area is determined by the radius alone. The cell problems are defined as before, but where Y 0 (x, t) = {y ∈ [0, 1] 2 | (y 1 − 1/2) 2 + (y 2 − 1/2) 2 > R 2 (x, t)}, G 0 (x, t) = {y ∈ [0, 1] 2 | (y 1 − 1/2) 2 + (y 2 − 1/2) 2 < R 2 (x, t)} and Γ 0 (x, t) = {y ∈ [0, 1] 2 | (y 1 − 1/2) 2 + (y 2 − 1/2) 2 = R 2 (x, t)}. Note that assuming radial symmetry simplifies the level set equation, no y dependence is needed. The model equations are also simplified through the tensors A u , A f , A g and K being cheaper to compute.
To summarize, we have considered a model for fluid flow combined with mineral precipitation and dissolution in a perforated domain. The model is non-isothermal and includes changes in the geometry of the pores, arising due to dissolution and precipitation. This leads to a model involving free boundaries at the pore scale, which are described by a level set formulation. We have applied a formal homogenization procedure and obtained upscaled effective equations for an idealized porous medium. The upscaled model describes the average behavior of the coupled system at the Darcy scale, but still honors the pore scale effects through effective parameters and cell problems. For a general pore geometry, the upscaled equations still depend on the microscopic variable through the level set equation. Assuming particular pore geometries may simplify the problem, but may result in less realistic models.