Simulation of Reactive Transport in Fractured Porous Media

Numerical simulations of reactive transport in fractured porous media require the solution of coupled physical and chemical processes that depend on the fractures. Such coupled processes are described by a system of nonlinear partial differential-algebraic equations, while strong heterogeneities characterise fractures. This paper presents an approach to simulate single-phase flow and non-isothermal reactive transport with mineral dissolution and precipitation in fractured porous media. Our numerical solution strategy is based on two ingredients. First, the model equations consist of coupled partial differential equations for the fluid flow, heat transfer and solute transport and nonlinear algebraic equations representing the chemical reactions. Second, fractures are explicitly represented and treated as lower-dimensional objects. The partial differential equations are discretised using finite-volume methods, and at each time step, we solve a nonlinear system of equations using Newton’s method. With numerical simulations, we illustrate our model’s ability to accurately describe the two-way interaction between coupled multi-physical processes and two- and three-dimensional porous media with intersecting fractures.


Introduction
Understanding the impact geochemical reactions, with some combination of fluid flow, solute transport and heat transfer, have on the subsurface is important since it has a wide range of applications. These include carbon sequestration (Addassi et al. 2021), geothermal energy extraction (Akaku et al. 1991;Todaka et al. 2004;Chen et al. 2018), nuclear waste 1 3 management (Socié et al. 2021), microbiofilm processes (Rittmann and VanBriesen 1996) and ion exchange and chromatography (Appelo 1996).
The subsurface has, in many cases, a complex fracture network. Fractures often serve as the preferred flow paths. They thus transport more rapidly chemical species that trigger chemical reactions. The chemical reactions alter the porous medium locally, e.g. through dissolution/precipitation of reactive minerals, impacting fracture conductivity (Salimzadeh and Nick 2019), or through accelerated degeneration of the geological structure (Socié et al. 2021). Such modifications affect the porous medium globally. Hence, fractures strongly influence the geochemical alteration of the porous medium.
This tight interplay between reactive transport processes and fractures makes the development of mathematical models and simulation tools difficult. Additionally, the coupling of reactive transport with fluid flow and heat transfer is characterised by a strong chemical and physical coupling. The flow and transport equations are described by coupled partial differential equations (Steefel et al. 2005;Cheng and Yeh 1998;Xu and Pruess 2001). If local equilibrium is assumed, the chemical reactions are represented by nonlinear algebraic equations (Yeh and Tripathi 1989;Kirkner and Reeves 1988).
Numerical strategies for solving the coupled chemical and physical equations are divided into the sequential approach and the global implicit approach (Yeh and Tripathi 1989;Steefel and MacQuarrie 1996;Zhang et al. 2016). In the sequential approach, the equations are split and solved sequentially, either once (non-iterative) or via iteration. These options are attractive since they are easy to implement and allow the usage of the most suited method for each equation. However, they require small time steps, even for moderately complex problems (Saaltink et al. 2001;de Dieuleveult et al. 2009), and the iterative option may require several iterations to converge (Saaltink et al. 2001). In the global approach, the equations are not separated but solved fully coupled. While this method is more challenging to implement, the advantages are that it often requires fewer iterations and permits larger time steps (Saaltink et al. 2001;Steefel et al. 2015). The performance of the coupling schemes is compared in, e.g. Steefel and MacQuarrie (1996), Saaltink et al. (2001) and Carrayrou et al. (2010).
Incorporating fractures in solution strategies implies the need for a representation that captures the transport and reaction effects within the fractures. One option is to use the multiple interacting continua model of Pruess and Narasimhan (1985). The idea of this model is to appropriately subgrid the matrix near the matrix-fracture interface to capture, e.g. chemical concentrations that propagate more rapidly through the fractures. This model is implemented in the simulator TOUGHREACT (Xu and Pruess 2001), which has been applied to simulate the geochemical interaction with fractured media in various applications (Xu et al. 2006;Dobson et al. 2003).
An alternative option is to explicitly represent the fractures in conjunction with the matrix (MacQuarrie and Mayer 2005; Berre et al. 2019). This option is commonly referred to as a discrete fracture-matrix model. In such models, the matrix and fractures are assigned individual equations. Further, the fractures are considered as objects one dimension lower than the matrix (see, e.g. Steefel and Lichtner 1998a;Fumagalli and Scotti 2021;Mac-Quarrie and Mayer 2005), a representation often called mixed-dimensional representation (cf. Fumagalli and Scotti 2021). In the lower-dimensional modelling of fractures, the fracture aperture is treated as a model parameter and not a geometric constraint (Berre et al. 2019;Fumagalli and Scotti 2021). This is beneficial since the aperture may change during the simulations without requiring regridding.
Examples using this approach include Steefel and Lichtner (1998a), who presented a multicomponent reactive transport formulation with mineral dissolution and precipitation in fractured media. In their model, the matrix-fracture coupling was through a continuity concentration condition and a zero-gradient concentration condition midway between the fractures. In a companion article (Steefel and Lichtner 1998b), they used this formulation to simulate the infiltration of hyperalkaline groundwater into fractures.
Another example is Fumagalli and Scotti (2021), who proposed a mathematical model and sequential non-iterative scheme for non-isothermal reactive transport for the geochemical model of Agosti et al. (2016). Here, the matrix-fracture interface was modelled using the dimension-reducing interface modelling of Martin et al. (2005).
This work presents a mathematical model and numerical approach to simulate singlephase flow and non-isothermal multicomponent reactive transport with mineral dissolution and precipitation in fractured porous media. The basis of our model is a discrete fracture-matrix mixed-dimensional framework, using the methodology of Martin et al. (2005) to define the dimensionally reduced governing equations for the fractures and the fracture-matrix interactions. The fluid flow, heat transfer and solute transport are modelled by conservation partial differential equations (PDEs) and discretised by finite-volume methods. The chemical reactions are represented by nonlinear algebraic equations, while the mineral dissolution and precipitation are formulated as a complementarity problem. To advance in time, we solve the governing equations fully coupled.
The article is organised as follows: section 2 introduces mixed-dimensional geometry and describes the governing model equations and the chemical system. We describe our numerical approach in Sect. 3 and present simulation results in Sect. 4. The numerical tests will show the convergence of our model and its ability to handle the interplay between challenging fracture networks and chemical and physical processes. Finally, we give some concluding remarks in Sect. 5.

Mathematical Model
This section describes the mathematical model for the reactive transport in fractured porous media. We begin, for simplicity, by presenting the governing model equations in the context of the porous matrix in Sects. 2.1-2.3. After introducing the mixed-dimensional geometry in Sect. 2.4, we expand the model PDEs in the mixed-dimensional setting in Sect. 2.5. Finally, changes in factors that couple geochemical transport to fluid flow and energy transport are addressed in Sect. 2.6.

Geochemical Equations
We consider a geochemical system that consists of aqueous, sorbed and precipitated chemical species and assume at any time chemical equilibrium. We divided the species into M 1 components and M 2 secondary species (Lichtner 1996;Kirkner and Reeves 1988;Yeh and Tripathi 1989), where the components form a linearly independent basis for the geochemical system. The chemical reactions between the components and secondary species are written in the canonical form (Lichtner 1996): (1)

3
The single species on the left-hand side is the secondary species, and the species on the right-hand side are the components. Finally, nk represents the stoichiometric coefficient for the kth component in the nth reaction.
Following the geochemical modelling of Yeh and Tripathi (1989), Kirkner andReeves (1988) andde Dieuleveult et al. (2009), we distinguish between the aqueous, sorbed and precipitated species by writing the above reactions for the various species. The reactions for the aqueous and sorbed species read: Here, c = (c k ) ∈ ℝ N c and s = (s k ) ∈ ℝ N s are vectors representing the chemical species of the aqueous and sorbed components. Further, = ( n ) ∈ ℝ N and = ( n ) ∈ ℝ N are vectors representing the chemical species of the aqueous and sorbed secondary species. Finally, the matrices S ∈ ℝ N ×N c , A ∈ ℝ N ×N c and B ∈ ℝ N ×N s contain the stoichiometric coefficients, with M ij being element (i, j) in the matrix M.
The equilibrium reactions are described by the mass action laws, which provide a link between the components and secondary species. Setting all activity coefficients to unity so that the activity and concentrations are equal and using the same notation for the chemical species and its concentration, the mass action laws for the aqueous and adsorbed reactions read: where K = (K n ) ∈ ℝ N and K = (K n ) ∈ ℝ N are the vectors of equilibrium constants.
The precipitate species dissolves or precipitates when they react with the aqueous components: where = ( n ) ∈ ℝ N is a vector of precipitate species and E = (E nk ) ∈ ℝ N ×N p are the stoichiometric coefficients.
If the precipitate species is present, we write the mass action law for the corresponding reaction as: where K = (K n ) ∈ ℝ N are equilibrium constants, and the activity of the precipitate species has been taken to be 1. If, on the other hand, the precipitate species is completely dissolved, the right-hand side of (3) is less than the left-hand side: Table 1 illustrates a schematic organisation of the algebraic relationships between the components, secondary species and stoichiometric coefficients in terms of the vectors c, s, , and , and the matrices S, A, B and E. With these variables, we can write a mass balance law for each component (Yeh and Tripathi 1989;Kirkner and Reeves 1988) where U and W are the total concentrations for the mobile and immobile components.

Flow Equations
We consider a single-phase flow in a porous medium, where the fluid flow is modelled by Darcy's law, together with the mass conservation for the fluid: In the above equations, v is the Darcy flux, K is the permeability, p is the pressure, is the porosity and f is the fluid density, with the subscript f referring to the fluid. Finally, is the fluid's dynamic viscosity and is, for simplicity, assumed to be constant.
We can insert Darcy's law into the mass conservation equation, yielding which we use to obtain the pressure.

Transport Equations
The transport of the solute species is assumed to happen due to advection but not dispersion. Denoting a concentration of a chemical species by u, we thus model the temporal evolution of the individual species by the transport equation where R u is a reaction rate term. We can compute a linear combination of the transport equations, under the assumption that the flow is independent of time and the chemical species (cf. de Dieuleveult and Erhel 2010). Furthermore, since the chemical reactions are at equilibrium, the reaction rate terms are eliminated when computing the combination. The outcome is the following transport equations, applied to the total concentrations where C = c + S T is the aqueous part of U. For the derivation, we refer to Yeh and Tripathi (1989); Kirkner and Reeves (1988). The mechanisms for energy transport are convection and thermal conduction, so we disregard external heat sources accounting for, e.g. heat production from chemical reactions. Moreover, local thermal equilibrium is assumed so that the temperature of the fluid and host rock are equal. With these assumptions, the temperature is modelled by the energy conservation equation: r and m = f + (1 − ) r are the heat capacity per unit volume and thermal conductivity, respectively. The subscript r refers to the rock. The specific heat capacity and thermal conductivity are denoted by b and , with the subscripts indicating if it is a fluid or rock property. Finally, r is the density of the rock. We assume these parameters are given and constant.

Mixed-Dimensional Geometry
We described the governing model equations in the porous matrix in the preceding subsections. In this and the subsequent subsection, we expand the model for fractures and their intersections. The collection of the model equations for the matrix, fractures and intersections results in a discrete fracture-matrix mixed-dimensional non-isothermal reactive transport model. We first introduce some notation for the mixed-dimensional representation of a D -dimensional ( D = 2, 3 ) fractured porous medium Ω . Following Martin et al. (2005) and Boon et al. (2018), we divide Ω into subdomains for the rock matrix, the fractures and fracture intersections. The subdomain for the matrix is a subset of ℝ D , the subdomains for the fractures are a subset of ℝ D−1 and the subdomains representing fracture intersection are a subset of ℝ D−2 . Additionally, in the case Ω ⊂ ℝ 3 , three fractures may mutually intersect at a single zero-dimensional point (intersection of intersections). Each subdomain is denoted Ω i , where the subscript i varies over all subdomains and is used to identify variables belonging to the ith subdomain. See Fig. 1 (left) for a 2D illustration of this decomposition.
Two subdomains whose dimensions differ by one communicate through an interface Γ j , where subscript j varies over the interfaces and is used to identify quantities belonging to the jth interface. Using the notation of Keilegavlen et al. (2021), we denote the higherdimensional subdomain by Ω h and the lower-dimensional subdomain by Ω l . The part of j Ω h and Γ j coincide geometrically. We also define Ŝ i and Š i as the sets of interfaces from Ω i towards its higher-and lower-dimensional neighbours, respectively. For example, for a subdomain Ω i representing a fracture, the set Ŝ i will represent the interfaces of the fracture to the neighbouring matrix, while the set Š i represents the interfaces of the fracture to possible neighbouring fracture intersection subdomains. For a subdomain Ω i representing the matrix, the set Ŝ i will be empty, while the set Š i will indicate interfaces of the matrix to neighbouring fractures.
To complete the description of the connection between neighbouring subdomains, we introduce projection operators for transferring variables between Ω h , Ω l and Γ j . The operators are depicted in Fig. 2. First, we have the trace operator, tr, which maps quantities in Ω h onto its internal boundary. Π is a projector that maps from a subdomain onto an interface. The superscript indicates whether the subdomain is of a higher or lower dimension, and the subscript points to the interface's index. Lastly, the operator Ξ maps from the interface onto a neighbouring subdomain, with a similar interpretation of the super-and subscript. Fig. 1 Left: Illustration of a 2D porous matrix ( Ω 1 ) with five fractures ( Ω 2−6 ), represented as red lines and two fracture intersection points ( Ω 7−8 ), depicted as black dots. Right: Visualisation of the connection between a higher-dimensional subdomain, Ω h , and a lower-dimensional subdomain, Ω l , via an interface, Γ j . The part of Ω h connected to the interface is an internal boundary j Ω h . The internal boundary, the interface and the lower-dimensional subdomain coincide geometrically but are shown apart for illustration Fig. 2 Illustration of the projection operators between two subdomains and a single interface. tr is a trace operator that maps quantities from inside Ω h onto its boundary, j Ω h . The Π l j and Π h j operators map from a subdomain onto the interface. Lastly, Ξ h j and Ξ l j map from the interface to the higher-and lower-dimensional subdomains, respectively

Governing Equations in Fractured Media
Utilising the notation from the previous subsection, we now present the governing PDEs in the fractures and their intersections; the geochemical equations (3)-(5) are as in the matrix.
For a subdomain Ω i of dimension d < D (i.e. the fracture and intersections), the fluid flow equations read In these equations, the nabla operator and the permeability K i operate in the tangential plane of the subdomain. Further, the scaling by the specific volume V i accounts for the dimension reduction (Stefansson et al. 2021). It is given by where a i is the hydraulic aperture of the fracture (Zimmerman and Bodvarsson 1996). We also point out that there is no tangential flow direction for zero-dimensional points, resulting in a vanishing Darcy flux (Boon et al. 2018).
The source term i accounts for the coupling with the higher-dimensional neighbour subdomain. With v j denoting an interface fluid flux, the source term on a single interface is Here, the scaling by the specific volume ensures that the dimension of the interface flux matches with the fluxes in the higher-dimensional neighbour (Stefansson et al. 2021). Since Ω i has several interfaces that couple it with its higher-dimensional neighbour (recall Fig. 1 (right)), the source term is a sum of fluxes over the interfaces, that is The higher-dimensional neighbour communicates with Ω i , via the interface fluid flux, through a Neumann boundary condition: The interface flux is modelled by the Darcy-type equation (Martin et al. 2005) Here, K j is an effective permeability in the normal direction of the lower dimensional subdomain, accounting for the conductance of fractures (Berre et al. 2019), and p h and p l are the pressure in Ω h and Ω l , respectively.
The equations describing the solute transport and heat transfer in the fracture and intersection subdomains, and at the interfaces have the same structure as the flow equations. For the solute transport equations, a sum of advective interface fluxes, j , enters the equations as a source term so that the system of equations will be: In Ω h , the Neumann boundary condition for the coupling takes the form: The advective interface flux follows the interface fluid flux direction as follows: Finally, the energy conservation equation in the mixed-dimensional setting reads: with the internal boundary conditions for j ∈Š i .
The interface convective and conductive fluxes, w j and q j , are calculated as

Coupling of Fluid Flow, Reactive Transport and Heat Transfer
The last set of equations we need is constitutive laws for the factors that couple the fluid flow, reactive transport and heat transfer. Our modelling framework allows using different constitutive laws to model a varying fluid density, such as linear dependence on temperature, pressure and ion concentration (Bringedal et al. 2014) or a combination of the temperature-dependent density of pure water and concentration, molecular weight and density of dissolved chemical species (Cheng and Yeh 1998).
In the present work, we model the fluid density of a slightly compressible fluid according to the following: Here, f ,0 is a reference fluid density, and the coefficients ĉ f and ̂f are the fluid compressibility and fluid thermal expansion, respectively. Finally, p 0 and T 0 are initial pressure and temperature.
The porosity is calculated from .
where n is the nth mineral volume fraction, and V n is the volume of the nth mineral. Further, nr is the mineral volume fraction of a non-reactive mineral and is only present in the matrix.
For alteration of the matrix permeability, we use a Kozeny-Carman equation, where K 0 and 0 are the initial permeability and porosity. Other choices are also possible, see, e.g. MacQuarrie and Mayer (2005). The fracture permeability is calculated by the cubic law We estimate fracture aperture as where a 0 is the fracture aperture when no reactive minerals are present, and A n is the nth mineral surface area. In this work, we set the surface areas to be constants for simplicity. However, constitutive laws where the surface areas depend on, e.g. porosity (Steefel and Lichtner 1998a) can also be used. At fracture intersections, we calculate the aperture as the average of its higherdimensional neighbours (Stefansson et al. 2021): In our model, the normal permeability at the fracture-matrix interface can be set independent from the tangential permeability. This flexibility can be used, for instance, to let the normal permeability determine which minerals are precipitated on the fracture walls. In the present simulation, we take the simpler option of setting the normal permeability equal to the permeability of its lower-dimensional neighbour: Finally, the equilibrium constants depend on temperature. We model this temperature dependency according to van't Hoff's relationship where K 0 is a reference equilibrium constant at the initial temperature T = T 0 and R is the ideal gas constant. Finally, ΔH 0 is the standard enthalpy of the particular chemical reaction and is assumed to be independent of temperature. However, as remarked by Cheng and Yeh (1998), we pointed out that this equation is only applicable if the standard enthalpies are available, and some data-based equations can be used to model K if this information is not available. We use the standard enthalpy values in Chang and Goldsby (2014).

Numerical Approach
We now present our fully coupled solution strategy to solve the system that consists of the PDEs (6)-(19) and the nonlinear geochemical equations (3)-(5), together with the constitutive laws (21)-(26).

Spatial Discretisation
The PDEs are solved using the method of lines, i.e. we first discretise in space. For spatial discretisation, we use cell-centred finite-volume methods. The advective flux terms are discretised by a first-order upwind scheme so that, e.g. the advective part of Eq. (11) is calculated as where sign(v kl ) is the sign of the flux moving from cell k to l at the previous iterate. We have a similar expression for the convective term in Eq. (13), where also the fluid density is computed from the previous iterate. For simplicity, the elliptic parts (the Darcy flux and conduction term) are discretised by the two-point flux approximation method (Aziz and Settari 1979). Finally, the fluid density in Eq. (9) is mapped from the cell centres to the faces by an upstream weighting.
For the lower-dimensional PDEs, we discretise the flux terms using lower-dimensional versions of the schemes described above Stefansson et al. 2021). Moreover, the interface fluxes are calculated using the projected subdomain variables. In particular, the advective interface fluxes (21) and (26) are discretised by an upwind scheme ).

Numerical Modelling of Geochemistry
To evaluate the conservation equations (5) and (6), we use the logarithmic formulation described in de Dieuleveult et al. (2009). We introduce the variables lc = ln(c) and ls = ln(s) . These variables guarantee that the component concentrations are always positive. Further, the mass action laws (2) and (2) can be rewritten as where Slc, Alc and Bls are matrix-vector multiplications. One issue with the logarithm formulation is when concentrations equal zero. We resolve this problem by approximation zero by exp{−30} . Lastly, the solubility condition (3) and (4) can be formulated as a complementarity problem and solved, in vectorised from, as (Kräutle 2011) The minimum function is a semismooth function, and the equation must be solved by a semismooth Newton method. We solve it with an active set strategy by defining the two sets (Kanzow 2004) so the minimum equation can be solved as where F A n and F I n are diagonal matrices with entries

Fully Implicit Scheme
To advance forward in time, we couple the semi-discretised PDEs with the geochemical equations (3)-(5) and the constitutive laws (21)-(26) at every computational cell. The resulting system of differential-algebraic equations is solved fully coupled, using the backward Euler method for temporal discretisation. Hence, at every time step, we solve a system of nonlinear equations in the form Here, G represents the residual form of the equations and Y n is a vector representing p, U, W, T, , a, K, K , at time step n. The nonlinear system is solved using a damped Newton's method (Dennis Jr and Schnabel 1996). The kth Newton iteration gives rise to a linear system where J is the Jacobian of G(Y n k ) and is calculated by forward automatic differentiation. The linear system is solved using the software superLU (Li 2005).
This strategy is relatively simple and straightforward but can be computationally demanding for problems where the number of mesh points is high and numerous chemical = K exp(Slc), reactions (and consequently solute transport equations) are present. One possible approach to resolve that issue could be to reduce the number of chemical transport equations (8) by applying the technique developed in Kräutle andKnabner (2005, 2007) and applied, e.g. in Hoffmann et al. (2010). Such an approach can also be applied to model solute transport and chemical reaction equations in the mixed-dimensional framework. Finally, the time step is adjusted following Saaltink et al. (2001): If the number of Newton iterations needed to converge is lower than a prescribed number, the time step is increased by a specified factor. Conversely, if the number of Newton steps is higher than a prescribed number or does not converge, the time step is decreased by a fixed factor. These are complemented with an upper and lower bound on the time step, i.e. dt min ≤ dt ≤ dt max .
We have implemented the presented model using the open-source simulator PorePy  to generate computational grids and for spatial discretisation.

Numerical Results
In this section, we test our numerical solution strategy on a synthetic test problem. The chemical reactions we consider were introduced by Bringedal et al. (2014) and are displayed in Table 2. This is a simple set of chemical reactions where the two minerals calcite ( CaCO 3 ) and anhydrite ( CaSO 4 ) might dissolve or precipitate and affect the permeability. The first three reactions are present since they affect the concentrations of three species in the last two reactions.
Using the Gaussian elimination procedure described in Steefel and MacQuarrie (1996), we find that the total concentrations are where H 2 O is treated as neither a component nor a secondary species and is assigned unit activity (i.e. it is incorporated in the equilibrium constant of the third reaction). We note that where are the stoichiometric coefficient matrices for the chemical reactions.
In the simulation examples, we consider a fractured porous medium filled with 20 mol∕m 3 of calcite and no anhydrite. We then inject a fluid, which causes dissolution of the calcite. The dissolved mineral releases Ca 2+ ions, which will react with present SO 2− 4 ions and trigger precipitation of anhydrite. An interesting point about our model is that, due to the chemical equilibrium assumption, the precipitation reaction happens simultaneously with the dissolution reaction. Hence, the dissolution front and the precipitation front coincide.
We estimate the mineral volumes using the mineral densities from Schön (2015), and the surface areas are set to two-thirds power of the computational cell volumes. The reference equilibrium constants are from Plummer et al. (1988) and are presented in Table 2.

Two-Dimensional Simulation-Convergence Study
We consider a convergence study in an isothermal and a non-isothermal context in the first two simulation examples.
The porous medium is initially at a state with a pressure p 0 = 1000 Pa and a matrix permeability K 0 = 10 −11 m 2 . We alter the initial state by injecting a fluid, based on a pressure of 7000 Pa, from the left-hand side of the domain. For the aqueous species, the initial and inflow values for SO 2− 4 and OH − are 10 and 1500 mol∕m 3 , respectively. The remaining initial and inflow values are so that the mass action laws and solubility conditions are fulfilled. The top and bottom walls are impermeable, while on the right-hand side, we impose the initial concentration and pressure values as boundary values.
The geometry is a 2 × 1 m two-dimensional domain, consisting of four fractures and one intersection, see Fig. 3 and is discretised with triangles. Starting with a coarse grid consisting of 187 cells for the matrix, 15 fracture cells and one intersection cell, we refine the grid five times. The solutions on the finest grid with 95972 2D cells, 318 1D cells and one 0D cell are used as reference solutions.

Isothermal Case
In this subsection, we verify our model in the isothermal context, i.e. we have T = T 0 at all times. Consequently, we do not update the equilibrium constants, and the fluid density varies only due to pressure (recall Eq. 27). Additional simulation parameters are listed in Table 3. Figure 4 shows the reference solution of the pressure and concentrations of Ca 2+ , calcite and anhydrite at the time 1300 s. In the simulation, the fractures are the preferential flow paths, and thus we see a more rapid transport of the aqueous concentration through the leftmost fractures. We see analogous fronts for the dissolution and the precipitation. Further, Fig. 5 shows the same reference solutions at the time 7000 s. The impact of the fractures as the preferential flow paths is again visible. The aqueous concentration propagates rapidly   Finally, Figs. 8 and 9 show the logarithm of the errors for pressure and the chemical species in the individual subdomains under the grid refinement at the time 7000 s. The error is calculated by projecting the coarse grid solution onto the reference grid and measuring the l 2 -norm of the difference between the reference solution and coarse solution. Lastly, the error is normalised by the number of cells in the reference subdomain (Stefansson et al. 2021).
In general, we observe first-order convergence. However, we also observe that the convergence is somewhat slow, and in the fractures, particularly the second one, there are minor kinks. We also note from Fig. 5 that all the calcite within the fractures is dissolved. This implies that the errors for this variable are always zeros, and while we cannot take the logarithm of the errors, it still is convergent. Modelling of reactive transport is difficult (Carrera et al. 2022) since it requires the solution of coupled PDEs and nonlinear algebraic equations (de Dieuleveult et al. 2009;Carrera et al. 2022). The obtained solution depends on several aspects, such as discretisation, the number of chemical species present (Addassi et al. 2021) and, in our simulation, the varying fracture aperture. These aspects strongly influence the result. Therefore, we consider the results in Figs. 8 and 9 to serve as a convergence verification of the model.

Non-isothermal Case
In this subsection, we carry out the convergence study in the non-isothermal context. The addition of the thermal effects increases the complexity of the simulation because it adds more couplings, particularly to the updates of the equilibrium constants. The varying equilibrium constants introduce an extra thermal front to the chemical front for the chemical concentrations.
The simulation parameters are the same as in Table 3. Regarding the thermal effects, the initial temperature is 573.15 K, and the inflow temperature is 543.15 K. The top and bottom borders are insulating barriers, and we use the initial temperature at the right border. Further parameters for the energy conservation equation and other temperaturerelated parameters are given in Table 4, and the standard enthalpy values are presented in Table 5. Figures 10 and 11 depict the reference solution of temperature and concentrations of Ca 2+ ion, the calcite and the anhydrite, after 1300 s and 5000 s, respectively. As in the isothermal case, the fractures are the preferred flow paths. Hence, the temperature and the ion concentration are transported faster through the fractures. Moreover, the chemical and thermal fronts are visible in the outcomes of the Ca 2+ ion concentration: The blue middle part represents the chemical front, while the leftmost part illustrates the temperature front. The  results for matrix and fracture permeabilities compared to their initial state and aperture clogging are not shown since they are indistinguishable from the isothermal case. Figures 12 and 13 show the logarithm of the errors in the matrix and the fractures for the variables at the time 5000 s. We again observe first-order convergence, albeit slow. Nevertheless, based on the discussion in the previous subsection, we again consider this to be a convergence verification of our model.
We finally observe in Fig. 11 what appears to be a non-physical oscillatory-like behaviour for the anhydrite and remaining calcite in the region that is affected by the cooling, in particular within the intersection fractures and in the matrix around these fractures. We zoom in on this region in Fig. 14. These instabilities were not present in the isothermal simulation case. Therefore, we believe they come from nonlinearities induced by the chemical reactions, in which the equilibrium constants vary due to an altering temperature in a porous media with different dynamically updated properties due to chemical reactions.

Three-Dimensional Simulation
In this section, we test our model on a three-dimensional problem. We consider the unit cube with two two-dimensional intersecting fractures; see Fig. 15.
The simulation parameters and initial and boundary values are as in Sect. 4.1, with the modification that the inlet and outlet for this problem are the planes x = 0 and x = 1 , respectively.
As in the 2D simulations, the reactive front is characterised by the calcite dissolving and the anhydrite precipitating. Figures 16 and 17 depict the calcite concentration and the anhydrite concentration, with the matrix concentrations represented as surface contours, after 700 and 1500 s, respectively. Similar to the 2D cases, the fractures are the preferred flow path, and the reactive front propagates more quickly in the fractures compared to the matrix (Fig. 16) and migrates into the matrix when the end of the fracture network is reached (Fig. 17). We also observe that within the fracture planes, the reactive front moves the fastest close to the fracture intersection.
This simulation illustrates that our modelling framework readily applies to 3D domains with intersecting fractures. We emphasise that this also applies to our implementation; beyond altering the computational grid, no code alterations were necessary to migrate from 2D to 3D.

Conclusion
This work presented a fully coupled approach to simulate single phase-flow and nonisothermal reactive transport in fractured porous media with mineral dissolution and precipitation. The model equations consist of coupled nonlinear partial differential-algebraic equations combined with several constitutive laws. Our solution strategy uses cell-centred finite-volume methods to discretise the PDEs and Newton's method to solve a nonlinear system of equations. Two-dimensional simulations have shown a convergence verification of the model and its ability to handle the two-way interaction between multi-physical processes and fractures. We also tested the latter point on a three-dimensional problem, demonstrating the full capabilities of our model. Author Contributions SIB wrote, edited and reviewed the entire manuscript, designed and ran all simulations, made all figures and analysed results. EK and IB supervised the work by reviewing the manuscript and took part in the result analysis.
Funding Open access funding provided by University of Bergen (incl Haukeland University Hospital). This work was financed by the Norwegian Research Council Grant Number 308733. The second author also acknowledges funding from the VISTA program, The Norwegian Academy of Science and Letters and Equinor.

Data Availability
No data sets were analysed in this study.

Code Availability
The run scripts for all simulations are available at https:// github. com/ taj004/ Fully_ coupl ed_ THC.

Conflict of interest The authors have no disclosed conflicts or competing interests.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.