Analytical Solutions for a 1D Scale Inhibitor Transport Model with Coupled Adsorption and Precipitation

In a previous publication (Sorbie and Stamatiou in Transp Porous Media 123:271–287, 2018), we presented a one-dimensional analytical solution for scale inhibitor transport and retention in a porous medium through a kinetic precipitation mechanism. In this process, a chemical complex of the scale inhibitor precipitates within the porous matrix and it then re-dissolves through a kinetic solubilisation process. Considering the re-dissolution of this precipitate in a one-dimensional linear system such as a reservoir layer or indeed in a laboratory core/pack flood, the flowing aqueous phase gradually dissolves the precipitate which is then eluted from the system. The most novel aspect of this previous analytical solution arose from the fact that, at a certain point in time (or pore volume throughput), the precipitate in the system was locally fully re-dissolved, forming an internal moving boundary between where no precipitate remained (closer to the system inlet) and where a precipitate was present (further into the system up to the outlet). In the current paper, we extend this work by presenting analytical solutions for the case where precipitation/dissolution occurs simultaneously with an adsorption/desorption interaction between the scale inhibitor and the rock surface, described by the nonlinear Langmuir isotherm. When examining this more complex problem in the flow scenario where the local precipitate is completely dissolved, several interesting analytical solution structures are obtained as a result of the internal moving boundary. Which of these structures occurs is rigorously categorised according to the solubility, the initial levels of precipitate and adsorbate, as well as the shape of the Langmuir isotherm. After the mathematical development of the analytical solutions, they are applied to some example problems which are compared with numerical solutions. Finally, a number of different generic features in the scale inhibitor effluent concentration profile are predicted and discussed with regard to practical field applications.


Introduction
When water is produced from petroleum reservoirs, it may form certain mineral "scale" deposits, such as calcium carbonate and barium sulphate, because of the aqueous geochemistry of the produced brine, or brine mixtures. In order to prevent the formation of this mineral scale, a production well may be treated with chemical scale inhibitors (SI) in a so-called squeeze treatment. To give a long "squeeze lifetime" in terms of the production of scale inhibitor in the produced brine, the chemical SIs must be retained in the reservoir formation. The main two mechanisms of SI retention within the porous medium are adsorption (denoted Γ ) and precipitation (denoted Π ) and, in many cases, this also occurs through a coupled adsorption/precipitation ( Γ∕Π ) mechanism.
A major problem for the petroleum industry is the formation and deposition of scale minerals on the downhole equipment and reservoir rock surfaces in contact with formation brine. The most common scaling minerals are calcium carbonate (CaCO 3 ) and barium sulphate (BaSO 4 ). Precipitation of these salts occurs when their concentration exceeds their solubility in the formation brine (Stiff and Davis 1952;Miles 1970). These deposits form on the tubing and in the near-wellbore formation and thereby cause a significant reduction of volumetric flow rates, resulting in a decline of oil production (Meyers et al. 1985). One-off chemical treatments to clean the well can restore productivity, but often only for a short period of time. For a more permanent resolution of the problem, the well is treated with scale inhibitors. These chemicals prevent the formation of scale minerals at either the nucleation or crystal growth stages of their deposition (Amjad and Demadis 2015). Scale inhibitor concentrations as small as 5 ppm can be sufficient to accomplish this (King and Warden 1989). For the scale inhibitor (SI) to be effective in reducing scale problems, it must be present in the brine occupying the pore spaces of the rock formation surrounding a producing well. This is achieved in a socalled squeeze treatment (Kerver and Heilhecker 1969;Vetter 1973), in which a volume of high concentration scale inhibitor is injected into the near-well reservoir formation. Production of the well is stopped, and a SI solution is injected. The well is then shut in for a period of time to allow the SI/rock interaction in the pore spaces to reach equilibrium. The SI chemical is retained in the rock formation mainly through two mechanisms: adsorption (denoted Γ ; Trogus et al. 1977;Ramirez et al. 1980;Sorbie et al. 1991Sorbie et al. , 1992 and precipitation (denoted Π ; Kahrwad et al. 2009;Sorbie 2010Sorbie , 2012Vazquez et al. 2010;Zhang et al. 2000). Once production is re-started following the shut-in, the SI chemical then desorbs or dissolves into the formation brine at a concentration which is sufficient to prevent the formation of scale crystals. More recently, there has been some revival of interest in the kinetics of inhibitor adsorption (Khormalia and Moghaddam 2007) and on salt precipitation in porous media (Safari 2014), and this still remains a subject of active research.
The desorption/dissolution of SI into the flowing aqueous phase is a kinetic process, and the treatment could be designed (using modelling) in order to increase the lifetime of a squeeze treatment and hence minimise costs. Mathematical models are invaluable in designing efficient SI squeeze treatments. The rate laws describing the adsorption/ desorption and precipitation/dissolution processes must be embedded in a transport equation for flow in porous media. In the present work, we will consider (simplifications of) the following three-dimensional Darcy-scale transport model including kinetic adsorption ( Γ ) and precipitation ( Π): Here, = x 1 , x 2 , x 3 ∈ ℝ 3 is the Cartesian coordinate vector labelling a point in the porous material, ( , t) denotes the (Darcy-scale) fluid velocity field, C( , t) is the SI concentration in the mobile phase, Π( , t) the level of SI precipitate in the pore spaces and Γ( , t) the amount of SI adsorbed onto the solid surface. The quantities C, Π and Γ are per unit volume of porous material. The dispersion coefficient D and porosity are assumed to be constants in this model. Precipitation/dissolution rates are given by Eq. (2). This formulation involves a rate constant , which can be related to the temperature via the Arrhenius equation (Yuan et al. 1993). The SI solubility C s in general also depends on field conditions such as temperature and pH (Malandrino et al. 1995), but here we simply assume that there is a critical temperature T cp such that C s is infinite for T < T cp and constant for T ≥ T cp . Equation (2) further ensures that (1) precipitation can only take place if C > C s and (2) no dissolution takes place if Π = 0.
Finally, the adsorption/desorption rate mechanism is described by Eq. (3). Next to a rate constant r a , it depends largely on a two-parameter Langmuir isotherm, given by Both rate Eqs. (2) and (3) appear as source/sink terms in the advection-dispersion equation (2), describing the transport of SI chemical in the bulk fluid. Quite often, with regard to field applications, model equations are written in spherical or cylindrical polar coordinates, which are particularly suitable for describing a near-wellbore geometry (see, for instance, Akanji and Falade 2019). However, in this paper, we will simplify the three-dimensional Cartesian equations (2)-(3) into a one-dimensional form appropriate for the analysis of core-flood experiments (Sect. 2). The reduced set of equations leads to a free-boundary problem which is analytically soluble (Sects. 3-5). A variety of cases are identified. Essential mathematical concepts such as the method of characteristics, weak solutions and shock conditions are utilised, and the reader is referred to Alinhac (2009), Ockendon et al. (2003, Holden and Risebro (2002), Lax (1957) and Smoller (1994) for more details.

Simplified One-Dimensional Model Equations
In a typical core-flood experiment, a homogeneous rock core of length L, cross-sectional area A core and porosity is saturated with a scale inhibitor (SI) solution of concentration C = C i (stage 1 in Fig. 1). We will assume that the core is oriented in the x 1 -direction. Precipitation occurs when the temperature of the core is raised above T cp and if C i > C s . At the same time, scale inhibitor adsorption/desorption takes place, and the interaction of these processes eventually leads to an equilibrium state with uniform levels C = C s , Π = Π 0 , Γ = Γ 0 = Γ eq (C s ) throughout the core. After a shut-in (no flow) period during which this equilibrium is reached (stage 2 in Fig. 1), SI-free water ( C = 0 ) is injected into the core at constant volumetric flow rate Q (stage 3 in Fig. 1). This translates into a constant linear fluid velocity v = Q∕(A core ⋅ ) ( cm s −1 ) in the x 1 -direction, so that Eqs.
(2) and (3) become one-dimensional and we can write x 1 ≡ x . Furthermore, we shall assume isothermal conditions, negligible dispersion ( D = 0 ) and equilibrium adsorption/desorption ( r a → ∞ ). This last condition implies that, for a given mobile phase concentration C(x, t), the adsorption level instantaneously becomes Γ eq (C(x, t)) . Then, by the chain rule, Γ∕ t = Γ � eq (C) C∕ t and Eqs.
(2)-(3) reduce to where the factor (1 − )∕ has been accommodated in the Langmuir isotherm coefficient A. The rock core can be represented mathematically as Ω ∶= ( In order to solve Eqs. (5)-(6) on the domain Ω , we use the following initial/boundary conditions on Ω: The boundary condition reflects the physical assumption that the concentration at the inlet drops to zero at the start of stage 3 in Fig. 1 ( t = 0 ) and remains so for all time.
If there is no precipitate present (i.e. Π 0 = 0 ), then Π∕ t = 0 and Eq. (5) can be written as C∕ t + V L (C) C∕ x = 0 , where the term V L (C) defines the Langmuir speed (see Fig. 2) of a concentration value C and is given by implies that the discontinuity in the data at x = 0 develops into a centred rarefaction wave, with each concentration value c ∈ [0, C s ] travelling at constant velocity V L (c) . This is illustrated in Fig. 3 and contrasted with the "pure dissolution" (i.e. Γ 0 = 0 ) solution obtained in a previous paper (Sorbie and Stamatiou 2018). Equilibrium desorption causes a retardation of concentration values with respect  to the "water front" moving through the system at velocity v. In the region V L C s t ≤ x ≤ vt , a bank of C = C s is sustained. At points behind this bank, the injected water has succeeded in bringing down the concentration level to some c < C s , triggering instantaneous desorption to Γ = Γ eq (c) , which in turn slows down the transport of the concentration in the x-direction. In contrast to adsorption/desorption, a precipitation/dissolution process aims to restore concentration to full solubility level C s . We saw that this leads to a discontinuity at the water front, followed by a steady-state region in which dissolution balances the horizontal transport of concentration (see Fig. 3). The question addressed in the current paper is how the two mechanisms combine (i.e. Π 0 ≠ 0 , Γ 0 ≠ 0 ) and interact.

Solution Method
From Eqs. (5)-(8) we see that Π(0, t) = Π 0 − C s t and therefore the SI precipitate at the inlet is completely dissolved at time t * , given by Construction of the analytical solution shows that Π( For t > t * , there is a boundary curve, denoted as x = Π (t) , which divides Ω into the subdomains Ω + (where Π(x, t) > 0 ) and Ω 0 (where Π(x, t) = 0 ). This curve will be shown to satisfy a differential equation of the form Therefore, the boundary curve is entirely determined by the concentration level on it, which itself is "supplied" by the solution in Ω 0 or Ω + (see Fig. 4). We will see that if d Π ∕dt > V L (0) at t = t * , there is a time > t * at which the boundary becomes a straight . This corresponds to a characteristic projection of Eq. (5) in the domain Ω 0 . For t * ≤ t < , the concentration on the boundary curve is given by the solution in Ω + , which in turn defines the solution in Ω 0 . For t ≥ , when the boundary coincides with a characteristic in Ω 0 , the concentration on it remains constant. This will lead to a travelling wave solution region in Ω + . Solutions for the various time regions are discussed in turn below for t < t * , while Π > 0 at the inlet, and for t > t * , when the quantity of precipitate drops to Π = 0 at the inlet.

Solution for t < t *
In this section, we first consider the initial period before the precipitate ( Π ) runs out at the inlet at time t * .

Concentration Profile
Equations (5)-(6) with Π > 0 reduce to a single PDE of the form a u∕ x + b u∕ t = c and can be solved analytically by the method of characteristics (see Alinhac 2009;Ockendon et al. 2003 for details). This method constructs the solution u from integral curves of the vector field (a, b, c) ∈ ℝ 3 defined by the coefficients of the PDE, which themselves may depend on u. In this context, it will be convenient to introduce so-called characteristic coordinates r = r(x, t), s = s(x, t) . These coordinates are chosen in such a manner that s is a parameter running along the integral curves of the vector field, while r is a parameter describing a space curve Σ 0 which encodes the initial/boundary data accompanying the PDE. The collection of integral curves intersecting Σ 0 then make up the solution surface u = Z(r, s).
For the specific problem of Eqs. (5)-(8), let us imagine that the boundary condition given by Eq. (8) is replaced with an initial condition prescribed on the negative x-axis (outside the domain Ω representing the rock core). We then consider the space curve Σ 0 = (r, 0, f (r)) ∈ ℝ 3 , with the parameter r running along the x-axis and the function f(r) describing arbitrary initial conditions. We can now solve the problem on the entire upper half of the x, t-plane and subsequently choose f(r) to satisfy the initial/boundary conditions in Eqs. (7) and (8).
The characteristic system of ordinary differential equations defining the integral curves corresponding to Eq. (5) with Π > 0 in terms of the parameter s is From Eq. (12) and Σ 0 , we obtain the general solution in characteristic coordinates: (1 + BZ) 2 (15) Z(r, s) = C s + f (r)e − s Moreover, Eq. (13) yields x(r, s) = vs + r . Introducing the notation Γ 0 ∶= Γ eq C s and Γ 0 � ∶= Γ eq � C s , we substitute Eq. (15) into Eq. (14) and integrate to find We now choose the arbitrary function g(r) such that t(r, 0) = 0 , so From Eq. (14), we see that t∕ s > 0 and hence We now express the initial/boundary conditions given by Eqs. (7) and (8) in the r,s-coordinates in order to determine the function f (r) . Using Eq. (18), we find that the initial condi-

the required function is
In order to single out a point on the discontinuity at r = 0 we introduce an auxiliary parameter ∈ [0, 1] and let f (0) = ( − 1)C s . This will allow us to distinguish between the characteristic projections emanating from the origin. At t = 0 , we have r = x and therefore Z(r, 0) = C(x, 0) = C s + f (x) . This is shown in Fig. 5. The feed condition C(0, t) = 0, t > 0 is accommodated by an initial condition C(x, 0) = C s − C s e − x∕v on −∞ < x < 0 . We will proceed to consider what is happening on the larger domain (−∞, L] × [0, ∞) and subsequently restrict attention to Ω.
Using x = vs + r we can eliminate s from Eqs. (15) and (16) by defining C (x, r) and t(x, r) as follows: According to the different components of the piecewise function f (r) , we distinguish solutions regions I, II and III (see Figs. 6, 7). The solutions in these regions will be denoted by C 1 ,C 2 ,C 3 in the x, r coordinates and by C 1 , C 2 , C 3 in the x, t coordinates.
(20) Region I r > 0 and hence C 1 (x, r) = C s . The characteristics t =t 1 (x, r) = (x − r)∕V L (C s ) are straight lines. We can invert to find r = x − V L C s t and therefore the explicit solution C 1 (x, t) = C s . There is a widening bank of constant concentration C = C s behind the advancing water front. This is purely due to equilibrium desorption.
Region II: here, r = 0 and hence C 2 (x, 0) = C s + f (0)e − x∕v . The auxiliary parameter ∈ [0, 1] is now used to pick a value of f(0), and we may write the solution as and The characteristics are the curves t =t 2 (x, ) , with the curve 0 ∶= t =t 2 (x, 0) and the line 1 ∶= t =t 2 (x, 1) bounding region II. It will sometimes be convenient to write t 2 (x, 0, ) and C 2 (x, 0, ) to emphasise the correspondence with r = 0. Given a point (X, T) in region II such as illustrated in Fig. 6, we can use Eq. (23) to find the unique value of such that t 2 (X, ) = T . The concentration level C 2 (X, ) at (X, T) is then found using Eq. (22). This parametric description is the closest we can get to an explicit solution of Eq. (5) in region II. As stated in Fig. 6, dissolution and desorption occur simultaneously in this zone of mixed behaviour.
Region III: here, r min < r < 0 , and hence C 3 (x, r) = C s − C s e − x∕v . Thus, for a fixed value of x, the concentration remains constant, while t(x, r) varies with r. In other words, the concentration in this region is independent of t and we can write We recognise this as the steady-state solution of Eq. (5). Since C 3 ∕ t = 0 , no (equilibrium) desorption occurs in this region. The characteristic projections are the curves This establishes the analytical solution for the mobile phase chemical concentration, C(x, t), for t < t * , as comprised of three components.

Precipitate Profile
We now use the concentration profile, C(x, t), derived above to construct an expression for Π on Ω . If Π > 0 , Eq. (6) yields This integral is straightforward to compute in regions I and III. With C 1 (x, t) = C s it follows that Π 1 (x, t) = Π 0 , since the initial condition Π(x, 0) = Π 0 , 0 ≤ x ≤ L must be satisfied. In region III, where the concentration is given by Eq. (24), we find for some unknown function y(x) that will be determined using the expression for the precipitate in the adjoining region II. An explicit formula C 2 (x, t) is not available here, but we can make the substitution t =t 2 (x, ) in Eq. (26) and let C 2 x,t 2 (x, ) =C 2 (x, ) , given by Eq. (22). Denoting the precipitate in these coordinates by Π 2 (x, ) , it may be shown that where the anti-derivative p(x, ) is We now choose the arbitrary functions y(x) and q(x) in Eqs. (27) and (28) in such a way that the precipitate profile is always continuous. Thus, we need Π 2 (x, 1) = Π 0 , which determines q(x) = Π 0 − p(x, 1) and hence Finally, by equating Π 2 (x, ) and Π 3 (x, t) on the curve 0 separating regions II and III (see Fig. 6) we find Substituting this into Eq. (27) and cancelling terms then yields the following expression for the precipitate profile in region III: For a fixed value P, consider the level curve t, x P (t), P on the surface and its projection t, x P (t) onto the x, t-plane, as indicated in Fig. 8. In region III, x P satisfies Π 3 x P , t = P and, denoting the corresponding concentration value C 3 x P by C P , the inverse function theorem can be used to show that Equation (33) holds the key to the construction of the entire solution. It describes the velocity of a precipitate value P in terms of the corresponding concentration level C P . It can be verified that this relationship also holds in region II, with C P = C 2 (x P , t) . A strictly mathematical proof of this fact is still at large, but it would probably involve the use of the chain rule in conjunction with the parametric solution components C 2 (x, ) , Π 2 (x, ) . Ideally, we would like to prove the stronger statement that Eq. (33) must apply to any new solution region which emerges, as long as the precipitate profile is required to be continuous. In the absence of a rigorous proof of this statement, we will henceforth assume that Eq. (33) is invariant between solution regions that are joined by a continuous precipitate profile.
4 Solution for t ≥ t * ; Case 1: 5 0 ≥ BC s 0 0 We now consider the time period after the precipitate "runs out" at the inlet; i.e. t > t * . For this discussion, we will employ the curve x = Π (t) introduced in Sect. 2.2. With Π 3 (0, t * ) = 0 , Eq. (6) implies that Π∕ t = 0 , which causes Eq. (5) to change its form to C∕ t + V L (C) C∕ x = 0 . Analogous to the argument given in Sorbie and Stamatiou (2018), it will be convenient to introduce also the point (5) in Ω + , we must always have C ≤ Π , which will effectively "slow down" the movement of C .
Note that Π (t * ) = C (t * ) = 0 and we consider the propagation speeds of the points Π and C to decide what is happening for t > t * . In the case of C , we observe that For the motion of Π , we use the assumed invariance of Eq. (33) between solution regions. Thus, for t ≥ t * , we suppose that In particular, C Π , t * = C C , t * = 0 and hence Comparing Eqs. (34) and (36), we see that Let us now suppose that Π 0 ≥ BC s Γ 0 (Case 1). Together with the constraint C ≤ Π , this implies the emergence of a joint root (35)]. This suggests there is some new region in Ω , next to region III, in which the concentration and precipitate components are given by travelling wave solutions. To specify these solutions, let z = x − U 1 t and C = c(z) in Eq. (5). This yields the ODE Equation (38) may be solved to give z up to an arbitrary constant 1 : Putting back z = x − U 1 t and using the condition that c = 0 at x 0 = U 1 (t − t * ) , we determine the constant Then, the new solution component, C 4 say, is implicitly defined as follows: To determine region IV in which Eq. (41) applies, we consider its intersection with Eq. (25). Substituting c = C 3 (x) into Eq. (41) and dividing by U 1 , we obtain We recognise Eq. (42) as the characteristic t =t 3 x, r TW with parameter value r TW defined by the relation Figure 9 shows region IV bounded by the (blue) line x 0 = U 1 (t − t * ) and the characteristic projection t =t 3 x, r TW emanating from t = t * . The characteristics in region IV are not shown, but it is possible to obtain them by using a different parameterisation of the Cauchy problem in which s = 0 on the curve x 0 = U 1 (t − t * ) instead of at t = 0 . The condition C x 0 = 0 then translates into an equivalent condition in terms of r, enabling us to determine the (new) function f (r).
Note that, for a fixed concentration value c, Eq. (41) describes its path in region IV. The initial position of c (in region III) is X c ∶= −v −1 ln 1 − c∕C s and its velocity is U 1 for t ≥ T c ∶=t 3 X c , r TW . This observation can also be used to construct Π 4 : the precipitate value corresponding to c is P c ∶= Π 3 X c , T c and the paths traced out by c and P c in the x,t-plane coincide for t ≥ T c .
In summary, if Π 0 ≥ BC s Γ 0 , the solution in Ω + consists of four regions I, II, III and IV in which the concentration is given by the components C s , C 2 , C 3 and C 4 , respectively. Region I is due to adsorption and region III due to dissolution, whereas both mechanisms are active in regions II and IV. In Ω 0 , Eq. (5) becomes C∕ t + V L (C) C∕ x = 0 , describing the case of pure adsorption/desorption. The characteristics here are straight lines of slope V L (0) > U 1 , determined by the data C(0, t) = 0 , and they run into the moving boundary x 0 = U 1 (t − t * ).
Example 1 This Case 1 ( Π 0 ≥ BC s Γ 0 ) solution is illustrated by a numerical example. Let A = 1, B = 10, C s = 0.1, v = 1, = 1, L = 1 . We then calculate Γ 0 = Γ eq C s = 0.05 and the Langmuir velocities V L (0) = 0.5 , V L C s = 0.8 . Note that BC s Γ 0 = 0.05 , so Case 1 (just) occurs if we choose Π 0 = 0.05 . Figure 10 shows the solution profiles consisting of C 1 ,C 2 , C 3 and Π 1 ,Π 2 , Π 3 at t = t * = 0.5 , with the adsorption curve plotted in red. Since Π 0 = BC s Γ 0 , we have � Π (t * ) = V L (0) and the travelling wave solution component emerges with this speed (see Fig. 11). We observe that C∕ x = ∞ at x = Π = V L (0)t , which is a (43) Fig. 9 Solution regions and characteristic projections for Case 1 result of the tangency of the boundary curve and the characteristic t =t 3 x, r TW at (0, t * ) . This blow-up is also visible in the effluent concentration plot in Fig. 12. Note that a larger value of Π 0 would result in a lower travelling wave speed and a more horizontal solution component, the derivative remaining finite everywhere. On the other hand, a lower value of Π 0 leads to qualitatively different solution (see Case 2 below). The effluent concentration flux can be used to prove explicitly that our constructed solution conserves the total amount of chemical present in the system. The profile can be divided into four parts from left to right, as shown in Fig. 12. Here, we have And, defining the constant ∶=C 2 (L, 0) = C s − C s e − L∕v , we compute Furthermore, denoting Eq. (41) as t = t TW (x, c) , observe that The integral on the right hand side evaluates to and, from Eq. (23), With Π 0 < BC s Γ 0 , we have Π � (t * ) > V L (0) and, during some time interval after t * , the boundary curve (32)]. Its slope, V Π = d Π ∕dt , is then given by Eq. (35) with C Π , t = C 3 Π . For the purpose of the following discussion, we will let y = C 3 Π , so that Note that, in terms of y, the equation Π 3 Π , t = 0 can be re-written as t = t Π3 (y) , where the function t Π3 ∶ ℝ → ℝ is Furthermore, we define the function T 3 ∶ ℝ → ℝ by evaluating t 3 (x, r = 0) (or, equivalently, t 2 (x, = 0) ) at x = Π = C 3 −1 (y) , so Finally, we also consider V Π = d Π ∕dt in terms of y: This will describe d Π ∕dt until y is the minimum of the following two values: (a) a ∈ 0, C s such that V Π a = V L a and V Π (y) > V L (y) for all y ∈ 0, a (b) b ∈ 0, C s such that t Π3 b = T 3 b and t Π3 (y) > T 3 (y) for all y ∈ 0, b Note that V Π (0) > V L (0) and t Π3 (0) = t * > 0 = T 3 (0) . We will now determine the values a , b and establish when a ≤ b and a > b . From Eqs. (9) and (53), it follows that The quadratic on the right hand side has roots Furthermore, it can be shown that (53), we deduce that V Π (0) > V Π C s = 0 and that V Π is continuous on 0, C s . It then follows that V Π has a maximum on 0, C s and, by Eq. (56), V L = V Π here. In order to determine whether this point is y a − or y a + , we write Π 0 = (1 + ) 2 B −1 C s −1 Γ 0 for some ∈ ℝ . If > 0 , the denominator in Eq. (55) is negative and y a + < 0 . The above argument now guarantees that y a − ∈ 0, C s . With < 0 , the denominator is positive and we need to examine whether y a − ≥ 0 or y a − < 0 . Substituting for Π 0 in the numerator, we find By assumption of Case 2, we also have BC s Γ 0 > Π 0 = (1 + ) 2 B −1 C s −1 Γ 0 . Hence, 1 + < BC s and y a − > 0 . Moreover, since y a − < y a + , it must be that y a − ∈ 0, C s . Finally, if the denominator in Eq. (55) is zero (i.e. = 0 ), we can take the limit → 0 to find y a − = (BC s − 1)∕2B > 0 , while the other root y a + is undefined. Thus, for all parameter choices satisfying Π 0 < BC s Γ 0 , y a − is the lowest value such that V Π = V L and we write a = y a − (see Fig. 13). We need to compare a = y a − with the solutions of the equation t Π3 (y) = T 3 (y) . From Eqs. (51) and (52), it follows that where equality/inequality occur simultaneously and In this case, it follows from Eq. (58) that t Π3 (y) > T 3 (y) for all y ∈ 0, C s . The same is true if D < 0 . Thus, for all parameter choices satisfying Π 0 > − , there is no b such that t Π3 b = T 3 b . Now suppose that Π 0 ≤ − . Then, Π 0 < Γ 0 and we have 0 < y b − < y b + < C s . It may be verified that For this particular choice of parameters, V Π = d Π ∕dt is equal to the Langmuir velocity of y = C 3 Π exactly at the time that solution region III disappears, i.e. when the boundary curve x = Π (t) enters region II. Since D = 0 , the boundary curve is tangent to the characteristic curve t =t 3 (x, 0) at the point x a = C 3 −1 y a − , a = t Π3 y a − (see Fig. 14). Fig. 13 The curves V Π and V L meet at y = a = y a − , given by Eq. (55) Since y a − satisfies Eq. (54), it must be that y a − > F , which is equivalent to Finally, this and Eq. (58) imply that Since t Π3 (0) > T 3 (0) and t Π3 , T 3 are continuous on 0, C s , the curves t = t Π3 (y) and t = T 3 (y) intersect on 0, y a − and we must have 0 < y b − < y a − ≤ y b + . In summary, if Π 0 > − , then b does not exist. Region III continues to exist for all time and we have a = y a − . Define a ∶= t Π3 a and x a ∶= C 3 −1 a . The boundary curve x = Π (t) now lies above the characteristic t =t 3 (x, 0) in the x, t-plane and V Π = d Π ∕dt varies according to Eq. (35) with y = C 3 ( Π , t) until t = a . We will refer to this scenario as Case 2a. It is to be distinguished from Case 2b, which occurs if Π 0 ≤ − . The boundary curve then intersects with the characteristic at the point At time t = b , solution region III disappears and for t ≥ b , the motion of Π is determined by the parametric solution in region II. By assumption of the invariance of Eq. (33), the velocity V Π = d Π ∕dt is now given by Eq. (35) with y = C 2 ( Π ( ), ) . We can therefore still expect to have V Π = V L when y = a = y − a . This will be verified in terms of the auxiliary parameter in Sect. 5.2.

Solution for Case 2a: 5 0 >ˇ−
For t > a , we have C( Π , t) = a , which is the concentration level "supplied" by the pure adsorption solution on 0 < x < Π (t) . Using the invariance of Eq. (33), we consider the constant velocity U 2 ∶= V Π a = V L a and look for a travelling wave solution of Eq. (5) in Ω + . Such a solution is of the form of Eq. (39) with U 1 , 1 replaced by U 2 , 2 . The constant of integration, 2 , is now determined by the condition that c = a at x = Π (t) = x a + U 2 t − a . We then obtain the solution Substituting c = C 3 (x) and dividing by U 2 , it may be verified that, for t > a , the travelling wave solution intersects with component C 3 along the characteristic t =t 3 x, r TW2 , where the parameter value r TW2 is defined by a =t 3 x a , r TW2 . This curve and the line x = x a + U 2 t − a enclose the new solution region IV, as shown in Fig. 15. The characteristics in this region are not shown, but may be obtained as described for Case 1. The boundary curve x = Π (t) is plotted in blue, and the characteristic t =t 3 x, r TW2 is tangent to this curve at x a , a . As mentioned before, the motion of Π for t * ≤ t ≤ a also results in a pure adsorption/desorption solution in Ω 0 : Equation (66) describes the path of a concentration value c ∈ 0, a . This lies at x = −v −1 ln 1 − c∕C s in region III until t = t Π3 (c) , when c = C 3 Π . For t ≥ t Π3 (c) , c moves at its own Langmuir velocity V L (c) . In Fig. 15, it is illustrated how a characteristic emanates from each point on the boundary curve between (0, t * ) and x a , a . The effluent concentration flux for Case 2a consists of five regions (see Fig. 19 in Example 2) and can be used to prove that the solution satisfies mass balance, as was done for Case 1.  Figures 17 and 18 show the solution profiles at t = a and t = 2.5 > a , respectively. The travelling wave solution component C 4 emerges at t = a and lies between x = 0.5845 and x = 0.8591 in Fig. 18. Its speed is U 2 = V L a = 0.3937 . Component C 3 is always present, which is more apparent in the effluent concentration flux plotted in Fig. 19. . 16 Solution profiles of Example 2 at t = t * = 1

Solution for Case 2b: 5 0 ≤ˇ−
The characteristic t =t 3 (x, 0) =t 2 (x, 0) now intersects the boundary curve x = Π (t) at x b , b , and region III disappears (see Fig. 20). For t > b , the boundary is determined by the solution in region II and can be found in terms of the auxiliary parameter ∈ [0, 1] , using the equation Π 2 Π , = 0. Noticing that t =t 2 (x, ) and employing the chain rule, we can determine the velocity of Π in terms of as follows: The concentration value at x = Π ( ) is C 2 Π ( ), , and we want to find out when To simplify this analysis, we introduce Since Π 0 < − < Γ 0 , we always have a + > 0 and it may be verified that Π ( ) < 0 for all ∈ a + , 1 , so that we can limit our attention to . This makes perfect sense in terms of the solution in region II: at t = b , the boundary curve hits the characteristic corresponding to = 0 . The boundary curve then intersects with subsequent neighbouring characteristics until = a − , after which Π ( ) stops being real-valued. Using Eqs. (67)-(72) (and after a lot of algebra), it can actually be proved that C 2 Π a − , a − = a = y − a [from Eq. (55)] and V Π a − = V L a . This verifies the invariance of Eq. (33) between solution regions II and III. As before, we now define x a ∶= Π a − , a ∶=t 2 0, a − . The characteristic t =t 2 x, a − is tangent to the boundary curve at x a , a . As in Case 2a (Fig. 15), a new region IV will exist for t ≥ a in which the solution is a travelling wave of velocity − . This region is enclosed by the curve t =t 2 x, a − and the boundary x = Π (t) = x a + U 3 t − a .
During the time interval t * , a , the motion of Π determines a nonzero solution in Ω 0 . For t * ≤ t ≤ b , this pure adsorption solution is given by Eq. (66). A similar relation applies for b ≤ t ≤ a : given an arbitrary concentration value c ∈ b , a , we solve c =C 2 Π c , c for c to obtain the time t 2 Π c , c . For t ≥t 2 Π c , c , the velocity of c is V L (c) . This is expressed in the following relation: The effluent concentration flux for this case again consists of five regions (see Fig. 23 in Example 3) and can be used to prove that the solution satisfies the principle of mass conservation.  71) and (72), we find a − = 0.0188 , Π a − = 0.356 . These can be used to determine a =C 2 Π a − , a − = 0.03124 and verify that this agrees with a = y − a . Then, we find a =t 2 Π a − , a − = 1.563 . Figure 21 shows the formation of the desorption tail in terms of characteristic projections. Figure 22 is a close-up of the solution profiles at t = 2 > a . The travelling wave component has velocity U 3 = V L a = 0.397 and is between x = Π = 0.5284 and x = 0.5511 on this plot.

Summary and Discussion
In this paper, we have derived analytical solutions for the scale inhibitor model with kinetic (non-equilibrium) precipitation and equilibrium adsorption described by Eqs. (5) and (6). We solved the Cauchy problem for this system with C( (5) is a non-homogenous, quasilinear PDE that incorporates the effects of precipitate dissolution and equilibrium desorption into the mobile phase. We were able to solve this equation using the method of characteristics in combination with the introduction of an auxiliary parameter , which enabled us to describe the Fig. 21 In Case 2b ( Π 0 ≤ − ), the boundary curve intersects region II and defines a solution in Ω 0 consisting of two components evolution of the discontinuity in the Cauchy data at x = 0 , t = 0 . This solution exhibits a mixture of the shock discontinuity found in the case of "pure precipitation" ( Γ eq = 0 ) and the rarefaction wave solution for the case of "pure equilibrium adsorption" ( = 0 ). The introduction of the parameter also allowed for the construction of the precipitate profile, Π(x, t) , and we showed that Π > 0 on Ω for all t < t * = Π 0 ∕ C s . An expression for the velocity of an arbitrary precipitate value P in terms of the corresponding concentration level C P was found [see Eq. (33)]. The assumption that this relationship is "invariant" between all solution regions leads to the equation of motion of the point Π [Eq. (35)]. If � Π (t * ) ≤ V L (0) ("Case 1"), we saw that a travelling wave solution emerged immediately and the concentration and precipitate profiles behind it were identical zero. On the other hand, if � Π (t * ) > V L (0) ("Case 2"), the solution was characterised by a nonzero pure adsorption tail in the region where the precipitate was used up. This tail continued to form until the velocity of Π became equal to the Langmuir velocity of the concentration value at x = Π . It was only at this stage that a travelling wave component began to emerge. The solutions thus constructed were tested for mass conservation by integration of the effluent concentration flux profile and further validated by comparison with numerical solutions in a few example cases.

Effluent Concentration Profiles
This work has primarily been concerned with the construction of the in situ solution profiles, without paying special attention to the finite length L of the rock core. However, in a laboratory setting, the principal measurable quantity is the effluent concentration C(L, t). The variation of this level with time can tell us which regime (desorption or dissolution) is more prominent. It also determines the "lifetime" of a squeeze treatment, the amount of time that the effluent concentration is at least equal to a threshold level C t below which scale inhibition becomes ineffective. In order to address these matters in terms of the various parameters, we briefly recall the three qualitatively different in situ profiles found in Sects. 4 and 5. They can be distinguished purely in terms of the (uniform) initial levels of precipitation ( Π 0 ) and adsorption ( Γ 0 = Γ eq C s ). To this end, we write Γ eq (C) = BΓ max C∕(1 + BC) and introduce the functions Here, Γ max is the maximum amount of scale inhibitor that can be retained on the rock surface through adsorption. If Π 0 ≥ P 1 Γ 0 , a travelling wave solution emerges at t = t * (Case 1). On the other hand, if If P 2 Γ 0 < Π 0 < P 1 Γ 0 , then the point ( Π a , a ) is in region III, so a at that time lies on the steady-state component (Case 2a). Finally, if Π 0 ≤ P 2 Γ 0 , then ( Π a , a ) is in region II (Case 2b). Prior to this, region III disappears at t = b such that C( The possible effluent concentration profiles are summarised in Fig. 24. We note that all cases with initial precipitate level Π 0 > P 2 Γ 0 have a "plateau" of constant concentration C s − C s e − L∕v as a result of the breakthrough of the steady-state component in the in situ profiles. For Π 0 ≤ P 2 Γ 0 , this only occurs if Π ( b ) = v −1 ln 1 − b ∕C s > L . It should be emphasised that Π 0 , Γ 0 , Γ max determine which qualitative case (1, 2a, 2b) occurs. The values a , b only have an additional explicit dependence on B. Once these are fixed, variations in , v and L cause the emergence and length of particular solution regions in the effluent profiles. In the next sections, we will identify situations in which a desired threshold concentration C t is to be maintained for as long as possible, taking into account practical constraints on the variables of the system.

Slow Versus Fast Dissolution
The rate equations in the full scale inhibitor deposition/retention model are governed by the adsorption/desorption rate parameter r a and the precipitation/dissolution rate parameter . There are four limiting behaviours (see Fig. 25), ranging from the fully kinetic case, when both rate parameters are finite, to the full equilibrium case, when both rate parameters are infinite. In practical applications, the rate parameters need to be considered in relation to the length of the system (reservoir or rock core) and the flow rate, which both determine the impact of the two mechanisms on the mobile phase concentration. For instance, in case of a core-flood experiment on rock core of length L, a low flow rate Q implies a Fig. 24 Possible effluent concentration profiles low fluid velocity v = Q∕A ( A = cross-sectional area, = porosity), which corresponds to a long fluid residence time L/v. A high but finite rate parameter can then lead to behaviour that is very similar to an equilibrium process. This is no longer the case if the flow rate is increased significantly. Thus, the rate parameters must be considered relative to the residence time of the fluid, and kinetic/equilibrium-type behaviour can always be achieved by increasing/decreasing v sufficiently. For fixed L/v, we say that a desorption/dissolution process is "fast" if it is very close to equilibrium behaviour (typically very high values of , r a ). On the other hand, the process is called "slow" if it deviates noticeably from equilibrium (low or medium values of , r a ). The analytical solutions found in this paper are for fast adsorption/desorption processes only (i.e. equilibrium adsorption), but capture both slow and fast precipitation/dissolution. Figure 26 illustrates slow versus equilibrium (fast) dissolution in Case 2a. Both in situ profiles are sketched at a time a + Δt , when the travelling wave component (IV) is present. Now, as is increased, a becomes closer to t * , which itself decreases due to faster dissolution of the precipitate. Regions II and III become narrower, while IV widens. At the same time, the concentration in these regions gets closer to C s . In the limit → ∞ , we have a = t * = 0 and regions II and III disappear completely. The characteristic separating regions III and IV (see Fig. 15) then collapse onto the characteristic x = V L C s t, and the travelling wave develops a discontinuity at x = V L a t . In order to emphasise the importance of the time a for the effectiveness of a squeeze treatment, we compare the effluent concentration profiles for "slow" versus "fast" dissolution. In the sketch in Fig. 27, it is assumed that the threshold concentration is C t < a and that Π ( a ) = v −1 ln 1 − a ∕C s < L , so that the travelling wave can be seen in the effluent profile. We observe that fast (equilibrium) dissolution sustains the concentration level C = C s before dropping off sharply to a and decreasing further down to the threshold concentration. The lifetime of the squeeze treatment now depends solely on the Langmuir speed V L (C t ) , as indicated by the red dashed line. If scale inhibitor is just as effective at threshold concentration as it is for higher concentrations, then the fast dissolution process is not economic at all in this case, because a pure adsorption/desorption squeeze treatment would achieve exactly the same lifetime. The initial mobile phase concentration for such a treatment just needs to be equal to C t , far below the solubility level C s . This evidently requires much less scale inhibitor than the coupled process, in which the injected concentration needs to be greater than C s in order for the precipitate to form (virtually always the case in practice). All this extra scale inhibitor is "wasted" if it just results in a higher return curve rather than delaying the breakthrough time of C t . The latter is achieved by the slow (kinetic) dissolution mechanism, mainly because of the resulting increase in t * . In order to optimise usage of the available precipitate, the plateau of concentration C s − C s e − L∕v should be as close to the threshold level C t as possible (just above, in fact). For a given rate parameter , such a low return curve can be obtained by increasing the flow rate. In terms of produced pore volumes, this results in the same lifetime while ensuring optimal usage of the scale inhibitor in the system. Moreover, we should bear in mind that the model with constant fluid velocity only applies to core-flooding experiments. In case of a producing well in the field, the fluid velocity is inversely proportional to the radial distance from the well. In order to predict the return curves accurately, we would of course need to derive an analytical solution for this radial model. However, qualitative predictions can already be made using the solution of the present linear model. Broadly speaking, the same type of solution components will appear in the radial setting, but they will be stretched considerably in the region close to the producing well, where the fluid velocity is very high. Also, if the same parameters are used as input for both models, the effect of dissolution in the radial case will be less than in the linear case, due to the shorter fluid residence time. In order to achieve similar concentration levels in both return curves, the flow rate in the radial case will then need to be lowered. This is even more true if precipitate is formed only at a certain distance from the producing, due to low temperatures in the near-well region and higher temperatures further into the reservoir.

Lifetime Increase Due to Precipitation
We noted above that if the steady-state component appears in the effluent profile, the available precipitate is used in the most efficient way when the flow rate is adjusted to make C s − C s e − L∕v = C t . If no adsorption were to occur, this yields a squeeze lifetime (in pore volumes) of exactly This lifetime is increased if adsorption is then considered. For example, in Case 1 ( Π 0 ≥ P 1 (Γ 0 ) ), the lifetime of the adsorption/precipitation squeeze treatment is found by calculating when the characteristic separating the steady-state and travelling wave components intersects with the line x = L . The characteristic is given by equation (42), from which it follows that where Γ 0 = Γ eq (C s ) , Γ � 0 = Γ � eq (C s ) and C t = C s − C s e − L∕v . The last term shows the addition in lifetime due to adsorption. The percentage increase with respect to T pptn is actually limited due to the inter-dependence of the parameters involved. For instance, an increase in the adsorption capacity Γ max will cause both Γ 0 and Γ � 0 to increase. However, for larger increases, we will have to increase Π 0 accordingly in order to stay in Case 1. This evidently caps the proportional increase in lifetime, since it increases T pptn . Similarly, changes in B or C s might increase Γ 0 , but simultaneously decrease Γ � 0 . Although only a rigorous analysis will reveal what combinations of parameters maximise the "extra lifetime term" in Eq. (79), we can make some general observations regarding this issue. Notice that, on the one hand, the lifetime tends to v∕L 1 + Γ � 0 pore volumes as → ∞ ( C t → C s ), corresponding to the retardation of the value C s due to desorption. On the other hand, it can be shown that, as → 0 , the last term in Eq. (79) becomes v∕Γ � 0 L + BC s Γ � 0 + BΓ 0 . Thus, for small , we have T ads/pptn ≈ T pptn ≈ Π 0 ∕ C s . In order to appreciate the benefits of inducing precipitation, we have to compare T ads/pptn with T ads = v∕L 1 + Γ � eq (C t ) , the lifetime achieved by a pure adsorption treatment. Figure 28 illustrates the variation of the ratio of T ads/pptn and T ads with the threshold concentration C t . Large increases in lifetime are observed particularly if the threshold concentration is much lower than the solubility. A similar benefit-analysis can be carried out in Cases 2a and 2b (also sketched in Fig. 28). Here, the advantage of inducing precipitation is less prominent for higher threshold values. This is already obvious from the effluent profiles in Fig. 24. In Case 2b for instance, the effects of desorption can outweigh dissolution so much that the steady-state component is never able to break through. This is because either the amount of precipitate is very low or the isotherm is initially very steep and then levels off, causing Γ 0 to be close to Γ max . Compared to pure adsorption, the squeeze lifetime is still improved, but less so than in Case 1 or 2a. Eventually, when is decreased or v is increased (and hence C t lowered) to such an extent that Π ( b ) < L , the steady-state will appear in the effluent profile, resulting in larger percentage improvements.

Lifetime Dependence on Desorption
In the previous sections, we considered how dissolution improves squeeze lifetime. In that discussion, it was assumed that the flow rate can always be adjusted to make C t = C s − C s e − L∕v , so that the lifetime is defined by the length of the steady-state plateau in the effluent concentration profiles (except possibly in Case 2b). However, in practice, it might not be possible to increase the flow rate to this extent. We could have a situation in which, at the fastest possible flow rate, the effluent concentration level C s − C s e − L∕v is still much higher than C t (as was shown in Fig. 27). The lifetime then depends significantly on the nature of the desorption mechanism, which is governed by the Langmuir isotherm. A steep rising isotherm with a high adsorption capacity Γ max will cause the lower concentration values in the pure desorption tail (present in Cases 2a and 2b) to propagate slowly. Dissolution improves this further, by delaying the advance of C t . However, as we saw previously, for very steep isotherms, the improvement on the lifetime achieved by a pure adsorption treatment is limited. In such circumstances, inducing precipitation is only worthwhile if a lot of scale inhibitors can be injected. This will increase Π 0 , making the improvement due to dissolution more significant. We then benefit from the "delay" caused by dissolution as well as the slow Langmuir velocity of C t .
There is another aspect involved with flow rate changes. In the discussion, thus far the adsorption/desorption process was always fast (at equilibrium), whereas precipitation/dissolution could be fast or slow (kinetic). This applies if the desorption rate parameter r a is very high in relation to the fluid residence time L/v. Thus, large increases in the flow rate could lead to adsorption/desorption becoming kinetic too, so that the bottom left corner of the "phase diagram" in Fig. 25 applies. In this case, the system can only be solved numerically, which we have done but space considerations prevent the inclusion of these results here. However, our findings highlighted the importance of the ratio ∕r a . For instance, if precipitation is slow with = 10 and adsorption is fast with r a = 100 , then making the flow rate ten times faster effectively decreases the rate parameters to = 1 and r a = 10 . Both processes are now slow in the sense that the solution profiles clearly show some deviation compared to the r a → ∞ case, and we found that the effect on squeeze lifetime will be rather limited. This may change for high values of ∕r a , when desorption is slow compared to dissolution. Then, if the threshold concentration is sufficiently low, the squeeze lifetime may be increased hugely due to very long pure adsorption tails.