Streamline simulation of a reactive advective flow with discontinuous flux function

Reactive transport in porous media with dissolution and precipitation has important applications in oil and gas industry and groundwater remediation. In this work, we present a simulation method for reactive flow in porous media of two salts that share an ion. The method consists of a front-tracking solver that uses the Riemann solutions of the underlying set of hyperbolic partial differential equations. In addition to the discontinuities stemming from the nonlinearities of the flux function, the flux function for the corresponding advection reaction equation also admits discontinuities for a heterogeneous medium. Here, we solve the Riemann problem for the governing nonlinear hyperbolic system with a discontinuous flux function. We use mass balance across the interface and the non-decreasing sequence of velocity of waves to seek the unique solution for this problem. Furthermore, a model is provided for mixing of streamlines at the well to estimate the amount of precipitate. In the use of streamline methods, we have modified the definition of time-of-flight to allow the model to be easily utilised for the heteregeneous case. The simulation method is applied to model dissolution through injection of an unsaturated fluid. It is shown that the first dissolution shock, which causes induced precipitation due to the co-ion effect, results in accumulation of precipitate at the well.


Introduction
Dissolution and precipitation of minerals in reactive flows are important phenomena in many industrial and natural processes. Water-flooding operations in oil and gas industry are one of these processes where evaluation of precipitation and dissolution of minerals plays a major role in reducing potential risks. In water-flooding, oil is driven out of the production well by injecting sea water into the reservoir at different locations. The two water streams, consisting of resident formation brine and the injected sea water, respectively, become chemically reactive as they propagate through the reservoir due to their different compositions.
Masoud Ghaderi Zefreh mg65@hw.ac.uk 1 Institute of Petroleum Engineering, Heriot-Watt University, Edinburgh, UK 2 SINTEF Research Institute, Oslo, Norway These chemical reactions involve precipitation of minerals in the reservoir and more importantly at the wellbore [33].
In this work, we present a streamline method for flows with two chemical reactions involving dissolution and precipitation of two salts relevant to water-flooding. In principle, the two salts that initiated the discussion of this problem were barium sulphate and calcium sulphate. These salts are common precipitates that form around the wellbore. It should be noted that for simulating scaling at the wellbore, the porosity away from the well will at most change by a factor 0.15% in the most extreme case. This happens when a stream of 4000 ppm barium mixes with another stream of water that is reach in sulphate. Hence, assuming a constant porosity seems to be a reasonable assumption. This holds true everywhere except for around the well where the accumulation of precipitates changes the flow path dramatically [33]. Assuming constant porosity, along with other relevant assumptions, allows us to decouple the flow and transport in this problem. In addition, considering the fact that the flow is dominated by trasport in these type of problems [2], we can effectively model the problem using a hyperbolic system of partial differential equations. This class of problems is a good candidate to be used in streamline simulators. However, it should be also noted that numerical solutions for a class of more general problems including diffusion has been studied and developed [1,23,24]. The models in these papers have been obtained from the work of [22]. More comprehensive models have been developed accounting for the effect of porosity change on the flow and temperature on the chemical equilibrium [8]. Others have studied the effect of temperature on the viscosity and density of the fluid [9] and multiphase flow settings [31].
Streamline simulation is a divide-and-conquer strategy that solves flow (pressure/flux) and fluid transport in separate steps. It uses a Lagrangian grid with sub-grid resolution for the fluid transport and exploits loose couplings in the equations to enable large time steps. In principle, a streamline simulator for porous media has three main components that work together; 1. A pressure solver that calculates the pressure field by solving the elliptic or parabolic pressure equation and evaluates the velocity field using Darcy's law. This solver works on the Eulerian coordinates. 2. A 1D solver that solves the evolutionary 1D transport problem for a given time step on the Lagrangian coordinates. 3. A mapping from Eulerian to Lagrangian coordinates and vice versa. The former uses the flow information (pressure and velocity field) to obtain the Lagrangian coordinates. This transformation decouples the 2D/3D problem into a family of 1D problems along streamlines. Since the 1D problems are independent, we can pick a representative set of streamlines. These streamlines are then passed to a 1D solver that propagates compositions or other quantities. The results of the 1D solver are mapped back to the Eulerian coordinates through averaging the quantities from streamlines.
The choice of the 1D solver depends on the physical processes involved and affects the overall performance. Front-tracking algorithms are very well suited for convection dominant flows, because they use the exact solution of the corresponding Riemann problem and can resolve displacement fronts and other discontinuities exactly [12,20]. A Riemann problem is an initial-value problem for a hyperbolic partial differential equation where the initial data is a step function. The solution to a Riemann problem involves propagating waves or traveling discontinuities. These waves may smear out and form a continuous variation as they propagate (in the nonlinear case) or they may retain their sharpness while traveling in the system (in both linear and nonlinear case).
For the the corresponding Riemann problem of the problem at hand, we use the solution provided by Helfferich in a porous medium with constant porosity (homogeneous) [18]. We extend Helfferich's solution to the heterogeneous case by reformulating it in hyperbolic terminology. The flux function for the corresponding reactive advection problem depends on the porosity. Hence, we need to solve a system of nonlinear hyperbolic equations with discontinuous flux function in the heterogeneous case.
Hyperbolic partial differential equations with discontinuous flux function appear frequently in simulation of flow in porous media. The theory is well established for solution to generic 1D scalar cases [3,4,15,16] and multidimensional scalar cases [6]. However, a general approach to solve these type of challenging problems for system of equations is not available. Existence and uniqueness of solutions for linear systems have been studied and some numerical techniques for nonlinear systems are suggested based on linearisation [17]. Another type of approximation for systems with fixed interface is suggested in [7]. The solution in the general case can depend on the small scale assumptions. The main point of our problem is that all the waves travel in the same direction and a global solution can be derived. Our solution for the 2×2 system is based on conservation of mass across the interface and the simple principle of increasing speed of propagation for the sequence of solution waves which for two waves means that the leading wave must propagate faster than the trailing wave.
The novelties of the approach in this manuscript are three main ideas. The first one is on the mixing rule at the well. Generally, at the well, where the flow converges, one obtains a maximum mixing of incoming fluid. In finite volumes, it is needed to have a finer mesh close to wellbore to capture this mixing. However, the mixing model discussed in Section 3.2 does not require a finer spatial resolution. The second is the spatial parameterization used along streamlines; to incorporate heterogeneities, this is chosen to be different than the usual time-of-flight. Finally, we have solved a 2×2 Riemann problem with discontinuous flux functions. The general solution for hyperbolic systems with discontinuous functions do not exist in a closed form [17] The structure of this manuscript is as follows: Section 2 describes the physical and chemical model. In Section 3, we provide a brief discussion on streamline simulators and develop a model for mixing of streamlines at the well. In Section 4, we provide the Riemann solution for the onedimensional problem. We provide a review of the solution provided by [18] together with our extension of the solution for the heterogeneous case. A few examples for different settings are illustrated in Section 5. Concluding remarks are presented in Section 6.

Model equations
The problem we address here is the injection of brine into a reservoir that already contains another brine with a different composition. The difference in compositions causes the two brines to react as the injected brine pushes the initial one. The flow field is obtained independent of the chemical reactions and transport of the materials. We assume there is no difference in the density and the viscosity of the injected and the produced brine. For the use case of this problem, i.e. dissolution and precipitation of sulphate minerals (BaSO4 and CaSO4), the change of density of the water is negligible. However, the reader should refer to [29] for stablity of density-driven flows and [8,9,31] for full models that incorporate the density change. The Darcy's law in the absence of gravity reads where K is the permeability tensor and μ is the viscosity of the solution. Mass conservation for an incompressible fluid in an incompressible rock requires Equations (1) and (2) with relevant prescribed boundary conditions are solved to obtain the velocity field. For the reactive part, we consider two ions A and B that react with another ion C to produce two salts AC and BC which precipitate, At constant temperature and when ions are in equilibrium, the products of the ion concentrations are constant, i.e., with n α being the concentration of component α in the fluid (e.g., moles per litre). The total amount N α of component α in a unit volume of a porous medium is where φ is the porosity of the rock and m α denotes the amount of component α in form of precipitate (AC and BC in Eq. (2)) per volume of rock (mole per litre). We assume that the volume of precipitates is negligible and thus the porosity does not change due to dissolution and precipitation of salts. In addition, we assume that the dissolution or precipitation of salts has a negligible effect on the volume and density of water. These two assumptions mean that precipitation can not drive flow and therefore limit the coupling between flow and transport. Using the charge balance equation in Eqs. (4) and (5) and the fact that the amount of C precipitate equals the sum of that of A and B, the two independent primary variables N A and N B are chosen. We assume the precipitates are stationary. Hence, the amount of component α that is transported by the fluid is vn α . The governing equations for the corresponding reactive transport problem read with the initial compositions

Streamline simulation
We choose a combination of streamline and front-tracking methods to solve Eq. (7) numerically for several reasons. First, the inherent numerical diffusion in finite volume methods will often overestimate or underestimate the amount of precipitate. As we will see in Section 4, the Riemann solution for the 1D case involves contact discontinuities in many cases. For contact discontinuities, the characteristic wave speeds are equal on both sides of the propagating discontinuity. Unlike multiphase displacement fronts, these discontinuities do not have any inherent selfsharpening mechanism and are therefore more susceptible to numerical smearing when solved with a low-order finite volume scheme. Second, general numerical methods are not computationally efficient for Eq. (7a). The flux functions n α are not differentiable across the boundary curves on the hodograph plane (see Section 4.1) for constant values of φ. This non-differentiability imposes restrictions on time steps for implicit finite volume methods because standard Newton solvers fail to converge for large time steps. Nonsmooth variations of φ in space renders the flux function discontinuous which in turn introduces more difficulty in implementing finite volume methods. Third, the solution for the simplified case is already developed [18] and we can develop the more general case based on that. Availability of the solution for the 1D Riemann problem makes the problem well suited for front tracking. Finally, streamline methods can often be more computationally efficient than finite volume methods in large domains and highly heterogeneous media [14].

Flow in the reservoir
In this section, we provide a streamline method to simulate the reactive flow in a porous medium. In a streamline simulator, we convert the multidimensional transport equations into a set of simpler 1D equations. Each 1D problem explains the propagation of quantities along one streamline. A streamline at an instance of time is a curve that is everywhere tangential to the velocity field. Therefore, two streamlines can never cross and the transport along a streamline is a 1D problem [14].
Herein, we will only solve the flow equations, Eqs. (1) and (2), on rectilinear grids. This means that we can use a variant of Pollock's method [30] to trace streamlines. In streamline methods for typical problems, e.g., tracer transport or two-phase flow saturation calculation, the streamlines are traced using Pollock's method by calculating the time-of-flight for an arbitrary particle at initial position x i . The time-of-flight τ * is defined as the time that is required for a particle at x i to reach a certain position x The integral in the above definition is along the particle path. Here, we use a modified definition of time-of-flight that does not include φ, i.e., The reason we use this definition is twofold. First, we can reconstruct the 1D form of Eq. (7a) along a streamline with no explicit dependence on φ. The result is similar to a truly 1D problem. To reconstruct the 1D form, we expand the divergence operator in Eq. (7a) and use Eq.
The result is the 1D problem in time-of-flight coordinates ∂ ∂t Second, for streamline methods, the porosity only occurs explicitly in the temporal derivative. Here, in contrast, in Eq. (7a) both terms n α and N α are functions of φ. This means that the 1D wave speeds in Eq. (11) will involve φ, whereas in other streamline methods, the effect of φ on the wave speed is accounted for in τ and not seen in the 1D equations.
After the propagation along the streamlines, the solutions are mapped back on to the spatial grid. The overall value of quantities for each cell is calculated using the contribution of each streamline that passes through that cell. The result of the 1D solver for each streamline is averaged to yield the concentration on each cell in the multidimensional domain using the following equation where nsl represents the number of streamlines and τ i,j is the incremental time-of-flight of streamline i at cell j and V i is the amount of fluid that streamline i carries [21].
To evaluate the amount of fluid V i that each streamline carries in 2D settings, we use the method suggested in [28]. We distribute the injected fluid on the outflow faces of the inlet cells in the Eulerian grid. An inlet cell is defined as a cell through which fluid is injected into the domain. In this regard, each streamlines is associated with a conceptual stream-tube. The boundaries of these stream-tubes are constructed by bisecting the location of streamlines on the outflow faces of the inlet cell. The width of each stream-tube is proportional to the amount of flow of its corresponding streamline. Hence, the fluid is distributed to the streamlines according to the weights that are determined by the width each stream-tube.
If the flow field changes for any reason, it has to be reflected on the streamlines. Hence, after a simulation from time t 0 = 0 to t 1 , the quantities are mapped to the Eulerian grid (using e.g., Eq. (12)), the flow field is recalculated and the fronts are propagated from t 1 to t 2 . This introduces a special type of numerical diffusion that may smear the fronts. For this problem, however, we will not consider any scenarios that changes the flow field and hence it remains constant. Herein, we run the simulation always from t 0 = 0 to any time t and therefore, the results are free from any error caused by numerical dispersion.

Mixing at well
In many circumstances, such as during the water-flooding, the flow is dominated by advection, and thus mixing through diffusion is limited in the reservoir. However, at the well all streamlines converge and mixing takes place. If this mixing results in an excessive production of precipitates (scales), the well will clog and the flow stops. Scaling requires very expensive remediation strategies. In this section, we provide a model for mixing of streamlines at the well to estimate the concentration and the amount of produced precipitates. We assume that all streamlines converge at the well and only at the well. We assume that the streamlines mix and reach equilibrium in a certain volume V w . In this volume, precipitates form if the molar content of a species exceeds its saturation value. The fluid mixture without the precipitate then leaves the well towards the surface. The subsequent dynamics are not addressed here. The mass balance for the components at the well reads as where f in and f out are the molar inflow and outflow of ions to the well, respectively. According to our mixing assumption, we set δt such that at the next time step, the well contents are completely swept out and replaced by the new mixture of streamlines. We approximate f in and f out as constant values from t k to t k + δt. Given the total flow rate q we obtain the time step δt from the following relation For the molar inflow we have where q i denotes the flux of streamline i and n i (t k ) is the vector of concentrations at time t k associated with streamline i. Note that in Eq. (15) we have used n and not N , because the precipitates do not move and it is only the fluids that mix at the well. The molar outflow is expressed via the following equation wheren(t k ) is the concentration of the equilibrated outflow from the well at time t k . The summation in the Eq. (16) yields the total out-flux from the well q. Using Eqs. (5a), (15) and (16) in Eq. (13) with φ = 1 yields with m = (m A , m B ) T . Inserting the value of δt (Eq. (14)) in Eq. (17) yields the total amount of components at the next time step The average value of concentrations in the fluid at the next time stepn(t k +δt) and the amount of precipitates m(t k +δt) is calculated using Eqs. (5a) and (36) with φ = 1 and N(t k + δt) obtained from Eq. (18). Note that we have not used any spatial resolution for the well in deriving Eq. (18). In fact, we only post-process the results from the streamline simulator. Note further that V w is a model parameter that represents the mixing volume in the well. To obtain the history of concentrationn(t k ) and the amount of precipitate at the well m(t k ), we need to iterate through Eq. (18) with the estimated value for δt. An interpolation might be required since δt and the time step for the streamline simulator are not essentially equal.

Front-tracking and Riemann solver
The initial conditions for Eq. (11) are piece-wise constant, because a streamline passes through several cells that potentially have different concentration of components. We first develop the solution for the case with a single discontinuity in the initial compositions, i.e., Here, we use subscripts l and r representing left and right compositions, respectively. Equation (10) with initial compositions given by Eq. (19) is referred to as a Riemann problem. In Section 4.3, we combine Riemann solutions to construct the global solution for several discontinuities using front tracking. The solution to the Riemann problem follows from Helfferich [18]. In the next subsection, we give a summary of his analysis. The solution and the analysis applies to the case where porosity of the medium is constant. This solution is extended to involve heterogeneities of the system in the second subsection.

Helfferich's solution for constant porosity
The solution to any 1D Riemann problem depends only on the flux function and the two compositions (states), N r and N l , of the initial condition. Any other composition (state) that appears is transient and we refer to it as an intermediate state. The solution is self-similar and is obtained through the method of characteristics (also known as Riemann fans). Accordingly, the solution is characterised by a set of waves that separate different states [26]. Any state can be represented as a point in the state-space, which is often referred to as the hodograph plane [32]. Constructing the solution in the hodograph plane allows us to trace the waves by following composition paths from the two initial states. Figure 1 shows the corresponding hodograph plane for the problem at hand. Four regions representing four fluid conditions are marked in Fig. 1. Region I is unsaturated while region II is fully saturated with both salts. Regions III and IV are partially saturated (i.e., only one salt is saturated). For a single saturated component A, Eqs. (4a) and (6) with n B = 0 yield the single component saturation concentration as for A and similarlỹ for component B. When both A and B are saturated, we need to solve Eqs. (4a), (4b) and (6) for n A and n B to obtain the saturation concentrations aŝ The separating curve between regions I and III in Fig. 1 is obtained when no precipitates are present, i.e., m A = m B = 0. Eqs. (4a), (5a) and (5b) can now be combined to obtain Similarly, is obtained for the curve separating regions I and IV. The border between regions II and III and between II and IV are defined by the straight lines respectively. It is now possible to mathematically define regions I to IV using Eqs. (23) to (26). The structure of the self-similar solution depends on the saturation of the components in both left and right compositions, i.e., their location in the hodograph plane. The propagating concentration variations for this system are either contact discontinuities or shocks. Contact discontinuities are elementary waves that are sharp in the absence of diffusion and smear proportional to the square root of time in systems with diffusion. The latter occurs because these waves lack any self-sharpening properties. Shocks on the other hand tend to balance the diffusion mechanisms with the hyperbolic self-sharpening properties of the shock. Hence, a shock retains its shape after an initial smearing [26]. For our problem, shock waves are retarded chemical fronts that form as a result of a complete dissolution of a salt. Here, the number of shocks is equal to the number of salts being dissolved [10]. Contact discontinuities are transport waves that advance a concentration variation whenever there is no chemical reaction happening, i.e. when the solubility product is less than k sp . A non-stationary contact discontinuity propagates with the same speed as the fluid flows. Another form of contact discontinuities appear if a salt is partially dissolved. These contact discontinuities are stationary waves [18].
To use the method of characteristics, we analyse the Jacobian matrix of the quasi-linear form of Eq. (11) (see details in Appendix). In summary, the eigenvalues in the unsaturated zone (region I) are both 1/φ because no reaction happens if both states are in this region and the two PDEs decouple to two independent transport problems. In region II, both eigenvalues are zero. Finally, if only one salt is saturated (regions III and IV), we obtain λ 1 = 0, λ 2 = 1/φ. We analyse the behaviour of the system in regions III and IV with the help of integral curves. Integral curves are defined as the curves which are everywhere tangent to the corresponding eigenvector r p of an eigenvalue λ p . They are given by the following equation There are no rarefaction waves in this system because the eigenvalues are always constant along the integral curves [26]. We denote the two types of contact discontinuities as 1-contact and 2-contact (or C 1 and C 2 ). The former is stationary whereas the latter waves travel with the fluid velocity (i.e. speed 1/φ). If a wave has endpoints in different zones in the hodograph plane, the structure of the solution is different from the ones already discussed. In general, a wave with the left composition in a region with less precipitate is a shock (for example from regions I to III or from III to II). The velocity of a shock,v, for a discontinuity with states N and N * is found by utilising the mass balance law across the discontinuity resulting in the Rankine-Hugoniot condition In solving Eq. (28), we need an additional condition to choose the physical solution out of all the possible candidate solutions. Such additional conditions are often referred to as entropy conditions. For our solution, we use the Lax entropy condition [25], which states that a shock is admissible if it satisfies either of the inequalities where W is either a shock (S) or a contact discontinuity (C).
In the rest of this section, we provide the solution for all the cases.
-Cases 1, 2: Case 1 refers to N l and N r being in region I. As stated before, the solution for this case is only a simple wave N l C 2 − → N r , and the two PDEs decouple.
Case 2 refers to the case where both states lie in region II. The solution for this case is similar in that it is only one wave, but now this single wave is stationary because both eigenvalues of the Jacobian matrix are zero. -Cases 3, 4: In case 3, both left and right states belong to region III and in case 4 they both lie in region IV. The solution for case 3 is illustrated in Fig. 2a. The solution in notation of waves and an intermediate state where N i is the intersection of the horizontal line from N l and the 2-integral curve from N r in region III (I 1 (N l ) and I 2 (N r )). Case 4 has an equivalent structure.
Cases 5 to 9 refer to the situations where the upstream composition has more precipitate than the downstream composition. They all have a similar structure: a stationary wave followed by a fast transport wave, because precipitates do not move with the fluid. The notation for the solution reads N l -Cases 5, 6: The initial condition for Case 5 is N l ∈ II, N r ∈ III (and N r ∈ IV for Case 6). The solution for Case 5 in the hodograph plane is illustrated in Fig. 2b. Cases 10 to 14 describe the solution for cases where the upstream is relatively fresher than the downstream, i.e. the left composition belongs to a region with less precipitate in the solution. The solution for these cases involves at least one shock S 1 or S 2 .

The intermediate state for case 5 is the intersection of
-Cases 10, 11: Case 10 refers to the situation where N l ∈ III and N r ∈ II. The solution for this case is a stationary wave followed by a fast shock N l Point N i is the intersection of I 1 (N l ) and H 2 (N r ) and lies in region III. Case 11 has a similar structure with N i in region IV. -Cases 12, 13: In Case 12, N l lies in region I and N r in III. The slow wave for the solution to this case is S 1 and the fast wave is C 2 . Hence, we expect to The solution for this case involves two shocks, N l − → N r , thus a combination of cases 10 to 14. Point N i can be either in region III or IV or it can disappear. Helfferich concluded that if N l is below the line passing through N r and the triple point, then N i will be in region III to satisfyv (S 2 ) >v (S 1 ). Figure 3a shows the solution for this case. On the other hand, if N l lies above the line, the admissible value for N i is in region IV. For the case that N l is exactly on this line the two shocks have the same velocity, merge and state N i disappears.
Two more special cases are considered as cases 15 and 16.
-Cases 15, 16: In case 15, we analyse the solution for N l ∈ III and N r ∈ IV, which is a combination of cases 7 and 13. This case is illustrated in Fig. 3b. The solution involves a compound shock-contact wave and hence there are two intermediate states N i and N * i . The full solution is N l i is the intersection of I 2 (N r ) and the boundary curve A BS between regions II and IV. Since N * i is on the boundary, it belongs to both regions II and IV. Therefore, the rest of the solution follows as case 10, i.e. N i = I 1 (N l ) ∩ H 2 N * i . The solution for case 16 follows a similar structure to that of case 15.

Solution for variable porosity
In this section, we describe the solution for the Riemann problem when the porosity is spatially discontinuous in addition to discontinuities in the concentrations. Hence, the problem is given by Generally speaking, the flux functions n α depend on the primary variables N α with porosity acting as a parameter, i.e. n α = n α (N A , N B ; φ). Hence, Eq. (31) can be decoupled into two separate problems. We aim to obtain the solution by transforming Eq. (31) into two systems and analysing the behaviour on the interface in between where we have discontinuity of the quantities. This is similar to the scalar case in which one finds the global solution using the minimum jump entropy condition across the interface [5,15,16]. Equation (31) is decoupled into two problems with n l (·, ·) := n(·, ·; φ l ) and n r (·, ·) := n(·, ·; φ r ). Given that N(x, 0 − ) from Eq. (32) and N (x, 0 + ) from Eq. (33) satisfy additional condition [17,20], each of the problems Eqs. (32) and (33) has a unique solution according to the previous section. It remains to find the states N l and N r to concatenate these solutions together. In addition, ∂φ ∂t in Eq. (31) suggests that there is a stationary wave accounting   11) and (19) with N l = (0.05, 0.05) T , φ l = 0.1 and N r = (0.3, 0.4) T , φ r = 0.2 for the discontinuity of porosity 1 . Thus, we seek the endpoints of the aforementioned stationary wave,W 0 , N l and N r that account for the porosity change. This wave glues the weak solution of Eq. (32) to that of Eq. (33).
We first find admissible values for N l . The solution to Eq. (32) should not have any waves with positive speed because the wave to the right of this solution, i.e.W 0 , is stationary. It is not possible to have any waves with negative velocity either, because the eigenvalues of the Jacobian matrix are always non-negative. Hence, the solution to Eq. (32) can contain only stationary waves. This means that the solution to Eq. (32) is a stationary wave and is to the left of another stationary waveW 0 between N l and N r . They both start at the same position x = 0 and as a result, the two waves merge and the state in between, N l , disappears. The conclusion is N l = N l .
To determine the state N r , we use the conservation of flux across the interface, i.e. n lA (N l ) = n rA N r , n lB (N l ) = n rB N r .
In Eq. (34), we are stating that only the dissolved ions in the fluid travel. This is in accordance with our assumption that precipitates are stationary. Unfortunately, this condition is not enough to ensure a unique value for N r since neither of the four mappings N A → n A (N A , ·), N B → n A (·, N B ), (·, N B ) is surjective. We look at this problem case by case depending on the values of n A and n B for N l , i.e. location of N l in the hodograph plane. We identify four cases depending on the location of N r .
1. If N l is not saturated with any of the salts, N r is uniquely determined using Eqs. N l belongs to region III in Fig. 1), then N r should also belong to region III so that the flux is conserved. To show this, we pick N l = (N lA , N lB ) T ∈ III, then solving Eq. (34) for N r yields  Eq. (23)). Hence, the set D completely lies in region III and we have infinite possible candidates for N r . In the following, we first show that the global solution to Eq. (31) is independent of N r . We select an arbitrary value for N r in the set D so that we can construct a solution for Eq. (33). The solution to Eq. (33) with N r saturated in A involves a stationary wave (possibly) followed by nonstationary waves. To show this, we consider four cases for location of N r and we refer to the corresponding case number in the previous section for details.
(a) N r ∈ I: N r In summary, if N r belongs to region III, the solution to Eq. (33) always involves the stationary wave C 1 . Thus, the two stationary waves in sequence, (W 0 ) and the first wave in the solution of Eq. (33) (C 1 ), merge and the state in between (N r ) disappears. In other words, the solution for this case is independent of N r . However, to be able to construct the solution for Eq. (33), we choose N r to be the intersection of the boundary curve A BU and the set D. Figure 4 illustrates solutions to the above four cases (a to d) with φ l = 0.1 and φ r = 0.3.
3. If N l is saturated with B and unsaturated with respect to A, the solution follows the same structure as in the previous case. State N r is chosen to lie on the intersection of B AU and the corresponding vertical line. 4. Finally, when the left composition is completely saturated (i.e., N l belongs to region II) the situation is analogous to the previous two cases. There are infinitely many possibilities for N r such that n l (N l ) = n r (N r ) (basically any point in region II). However, the Riemann solution with left composition in this region always contains a stationary wave, c.f. cases 2, 5, 6 and 9 in Section 4.1. Hence, the solution is independent of the state N r . The simplest choice for N r is N r = φ rnA φ rnB (the intersection of all of the separating curves).
In summary, we have proved the admissible set for N l is {N l } and the admissible set for N r in the set D. In addition, for any values of N r in D we obtain the same solution, because N r does not appear in the solution.

Front tracking
Sections 4.1 and 4.2 described the solution for a single discontinuity whereas the initial condition for Eq. (11) has several discontinuities. Hence, initially we need to solve several Riemann problems. The solution waves from these Riemann problems might collide at later times in different Fig. 9 a History of the concentration of the sampled water (dashed for n A and solid for n B ) for homogeneous medium (φ = 0.2). b History of precipitate of A, m A (dashed blue) and precipitate B, m B (solid red) at well in mole/volume with V w = 0.1. As soon as B is dissolved, due to the co-ion effect, precipitate A starts to accumulate and reaches a maximum dependent to V w before it dissolves Fig. 10 a Porosity and b permeability field (−log(K))) of 20th layer of SPE10 model from MRST [11] locations. Every time two or more waves collide, a new Riemann problem appears and needs to be solved. The process of following each front to calculate the time and location of collisions is called front tracking [13,19]. Figures 5 and 6 illustrate the front-tracking idea in the characteristic plot. Characteristics are curves in x-t plane along which the solution remains constant. Figure 5 shows the characteristics for a problem with given initial condition. Each of the lines in Fig. 5 represents a wave that travels with constant speed equal to the inverse of the slope of that line. In the x-t characteristic plane, these collisions are represented by intersecting lines. Figure 6 illustrates an example for intersecting characteristic curves. In this example, an unsaturated solution is used to dissolve precipitates in a homogeneous system. The system initially contains a bank of precipitate B (x ∈ [0.05, 0.4]) followed by both precipitates (x ∈ [0.45, 0.6]) and precipitate A (x ∈ [0.6, 0.75]) with gradual transition between each bank.
In our implementation, we have used two doubly linked lists. The first list (X-list) contains all the discontinuities and is sorted by the location of each discontinuity. Initially at t = 0, we have only the X-list. Each discontinuity in the X-list is represented by the solution waves of the its corresponding Riemann problem. The collision time for each two neighbouring waves are computed and the collisions are stored in the second list (T-list). The T-list is sorted with respect to the collision times. All waves are propagated until the first collision (from T-list) occurs. The Riemann problem for the corresponding collision is solved and the solution is replaced by the colliding waves in the X-list. The T-list is updated based on the collision of the new waves with their neighbours. The algorithm stops when there are no collisions or when the next collision happens after the simulation time.

Illustrative results
In this section, we present the results of two simulations of injection of unsaturated fluid to dissolve two salts with different solubility product in a 1/4-five-spot-pattern. More specifically, we take k A = 1 and k B = 2, c.f. Fig. 1, with constant concentration of species for injection fluid N in = (0.05, 0.05) T . The boundary conditions for the pressure equation for both cases are similar (constant pressure at producer with a source term of 1 pore volume of fluid per unit time injected at the injection well). The two simulations are different in that the first is in a homogeneous medium whereas the second is for a heterogeneous rock. For both cases, the domain dimensions are 366 × 671 and the simulation time is set to 6 pore volumes with dt = 0.01 pore volume. For the homogeneous case the grid dimensions are 120 × 220. Figure 7 shows the propagation of injection fluid and evolution of dissolution fronts in a homogeneous medium with constant porosity 0.2 and permeability 100 milli Darcy using 120 streamlines. The initial condition of the medium is (0.4, 0.3) T . The 1D version of this problem is explained in case 14 of Section 4.1. As illustrated in Fig. 3a and explained in Section 4.1, the intermediate state N i has a higher concentration of N A . The reason is that the faster shock S 2 dissolves salt B and salt A starts to precipitate due to the co-ion effect. In the 2D case of this problem, the streamlines mix at the well resulting in even more precipitate of A. These precipitate do not leave the system until an unsaturated mixture arrives and dissolves the precipitates. In fact, if the arriving streamlines are not saturated and they have different compositions, they may induce more precipitates.
The accumulation of precipitates at the well through the mixing of streamlines is best illustrated in the plane spanned by the secondary variables n A and n B . Figure 8 shows an example. The values of n A and n B are bounded by zeros (horizontal and vertical axes) and the projection of B AU and A BU from the hodograph plane into this plane. The concentration of the components for any fluid (n) can be represented by a point in this domain including its boundaries. Now consider for example, the streamlines 1 and 5 in Fig. 7 after 1.5 pore volumes have been injected. Streamline 1 at the well is saturated with both components and the concentration of the components is n A ,n B T . Along streamline 5, the fluid is unsaturated with B but saturated with A. The two fluid concentrations are illustrated in the n A -n B plane. The domain is not convex and therefore any linear combination of n 1 and n 5 is outside the domain and implies an oversaturated fluid. As a result, more precipitates are produced and the precipitates accumulate and potentially clog the well. In a finite volume-based simulation, precipitates are distributed over the whole cells containing the well due to the numerical dispersion. This yields an underestimate of the amount of precipitate at the well. Figure 9 shows the produced water and the amount of accumulated precipitate at the production well.
As an example with heterogeneous porosity, we take the 20th layer of the SPE10 model from [11,27] (c.f. Fig. 10). As discussed in Section 4.2, the discontinuity in porosity causes the flux function to be discontinuous and thus the solution to the Riemann problem may involve an additional constant state. For the initial condition, we set the fluid concentration of species (n α ) equal for all blocks to establish an initial chemical equilibrium. The amount of initial precipitate for each block is set to be proportional to the porosity of that block. Hence, a cell with more void space would accommodate more precipitates initially. Figure 11 illustrates propagation of the fronts (n α ) in the whole domain at three times using 380 streamlines in 60×220 grid cells. Figure 12 shows produced water at the production well.

Conclusions
We have presented a streamline method for flow in porous media involving two chemical reactions with dissolution and precipitation. The two reactions share an ion and thus the common-ion effect determines the saturated concentration levels. The method is based on a front-tracking algorithm that uses the Riemann solution along each streamline. To apply the method for a heterogeneous medium, the 1D solution is extended to take into account discontinuities of the flux function for the underlying hyperbolic PDE. The unique entropy satisfying solution for the hyperbolic system is sought through honouring conservation of mass across the interface and the monotonic increase of wave speeds.
To predict scaling at the well, a model for mixing streamlines at the well is provided. The well model postprocesses the results from streamline simulator. Through this model, we have shown how streamlines with different concentrations arriving at different times result in high accumulation of precipitates at the well.
This problem is pathological to handle with finite volumes especially because of the characteristics of dissolution/precipitation function. These terms make the flux function non-differentiable (and discontinuous in the heterogeneous case). Other simple numerical methods may add too much numerical dispersion in the solution.
However, we see that most of the waves in 1D solution are contact discontinuities that are not self-sharpening. Hence an over dispersive method may be inaccurate.
Funding information This research is financially supported by the Energy Industry Doctorate Program Grant ETP136 of the Energy Technology Partnership, Scotland and Energi Simulation.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http:// creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Appendix: Secondary variables
In this section, we derive the expressions for the secondary variables in terms of the primary variables N A and N B . The secondary variables n α are used for the analysis of the corresponding eigen-problem and to identify wave types.
For the secondary variables n A and n B , we note that they are continuous but not necessarily differentiable across the borders. In general In region II, we have saturation with respect to both ions and therefore n α,II =n α .