Analytical Solutions of a One-Dimensional Linear Model Describing Scale Inhibitor Precipitation Treatments

The deposition of mineral scales such as barium sulfate and calcium carbonate in producing oil wells is a well-known problem in the oil industry, costing many millions of dollars per year to solve. The main preventative measure for managing downhole scale is to inject chemical scale inhibitor (SI) back into the producing well and out into the near well reservoir formation in a so-called “squeeze” treatment. The scale inhibitor is retained in and subsequently released from the formation by the two main mechanisms of “adsorption” (Γ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Gamma $$\end{document}) and “precipitation” (Π\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Pi $$\end{document}). As the well is brought back onto production, the scale inhibitor then desorbs (adsorption) or re-dissolves (precipitation) and the low SI concentration that is present in the produced water effectively retards the scale deposition process. A complete model of SI retention must have a full kinetic Γ/Π\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Gamma {/}\Pi $$\end{document} model embedded in a transport model for flow through porous media. In this paper, we present a subcase of this model involving only kinetic precipitation (Π\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Pi $$\end{document}). The simple quasi-linear problem with an infinite source of precipitate is straightforwardly soluble using conventional methods for a precipitate described by a solubility Cs\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C_\mathrm{s}$$\end{document} and a dissolution rate κ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\kappa $$\end{document}. However, the problem with a finite amount of precipitate is more complex and novel analytical solutions are presented for the (transient) behavior for this case. The mathematical difficulty in this latter system arises because, when the precipitate is fully dissolved close to the system inlet, a moving internal boundary develops along with some related flow regions defined by the parameters of the problem. The problem is solved here by making certain assumptions about the internal moving point at position αΠ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha _\Pi $$\end{document}, where Π=0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Pi = 0$$\end{document}, and we derive an expression for the velocity of this point (dαΠ/dt\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathrm{d}\alpha _\Pi {/}\mathrm{d}t$$\end{document}). From this, we then build the solution for all possible regions which may develop (depending on the problem parameters). Understanding the behavior of this idealized system gives us some practical formulae for precipitation squeeze design purposes. It also serves as an important set of reference solutions in the search for analytical solutions of more complex cases of the Γ/Π\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Gamma {/}\Pi $$\end{document} model which we will investigate in forthcoming papers.


Introduction
Chemical scale inhibitors (SIs) are commonly applied in oil fields to prevent the formation and deposition of "scale" minerals such as barium sulfate (BaSO 4 ) and calcium carbonate (CaCO 3 ). These scale inhibitors are applied in "squeeze" treatments that involve stopping a producing oil well (which is usually producing oil and water) and then pumping a chemical SI solution back into the well. After a short "shut-in" period, the well is put back onto production again. The SI chemical then desorbs or dissolves into the produced water stream at a low concentration (typically from 1-50 ppm) over a long period, which is sufficient to prevent or slow down the formation of mineral scale. It is well established that scale inhibitors (SIs) are retained within porous media by the two main mechanisms of "adsorption" and "precipitation". A large number of publications have appeared over the last three decades describing these mechanisms in some detail. For example, SI adsorption can be described by a generalized adsorption isotherm, (C), and precipitation may be modeled by a solubility function, (C), and a dissolution rate constant (denoted as κ or r 4 ) (Sorbie et al. 1992;Yuan et al. 1994b;Malandrino et al. 1995). Some researchers have related the SI adsorption mechanism based on the specific mineralogy of the (sandstone) rock (Sorbie and Gdanski 2005). Other workers have described SI retention by both an adsorption mechanism at low [SI] and a precipitation/dissolution mechanism at higher concentrations based on the solubility of the various Ca_SI salts that are formed (Kan et al. 1992(Kan et al. , 2004. All of these mechanisms are appropriate in certain circumstances, and the issue is not which approach is correct, but when each mechanism is correct and should be applied. Work on unifying the various approaches has led to a complete formulation of a generalized kinetic coupled adsorption/precipitation ( / ) model (Sorbie 2010). This model is fully dynamic and consistent; that is, it can model the full non-equilibrium (kinetic) adsorption/precipitation process such that this correctly limits in a consistent manner to the equilibrium case. Experimental work has broadly validated these more complex coupled adsorption/precipitation ( / ) models (Ibrahim et al. 2012a, b). Some time ago, Hong and Shuler (1988) published analytical solutions for equilibrium adsorption scale inhibitor flow through porous media described by a Freundlich adsorption isotherm. A considerable amount of more recent work studying precipitating SI systems has appeared in recent years for various type of scale inhibitor species. Shaw and Sorbie (2015) have presented results on the stoichiometry and modeling of mixed Ca/Mg /phosphonate scale inhibitor systems, in which the Ca/Mg salt of the SI is sparingly soluble, and details of this precipitation/dissolution behavior are discussed by Valiakhmetova et al. (2017). The precipitation behavior of polymeric types of SI inhibitor, such as PPCA (poly phosphino carboxylic acid), has also been studied (Farooqui and Sorbie 2016), where again the Ca salt of the polymer precipitates. The subtlety of polymer/Ca precipitates is that the molecular weight plays an important role (Farooqui and Sorbie 2016, and references therein).
The current work presents analytical solutions for precipitation squeeze treatments, describing the precipitation/dissolution process by a simple kinetic model. The quantity of precipitate locally is denoted , which is taken to be in the same units as solution concentration (e.g., M). If there is an infinite source of precipitate, the model effectively reduces to a single transport equation and is very straightforward to solve. More complex behavior is observed when the amount of precipitate is finite, so that the dependent variables C and are coupled together, and it is this case we solve here. For the purposes of this paper, we will assume that a uniform finite precipitate level 0 is initially present in the system and that a constant flow rate Q is used. An interesting composite solution of the model is obtained when the process is allowed to evolve long enough for the precipitate to be used up entirely. This solution constitutes a base case that may be extended to a more complex version of the model in which varying flow rates are used (Stamatiou and Sorbie 2018b). The base case is also instrumental in understanding numerical solutions of the general adsorption/precipitation model and indeed for finding an analytical solution of this model in the case with kinetic precipitation and equilibrium adsorption (Stamatiou and Sorbie 2018a).

Simple Bulk Kinetic Model for Scale Inhibitor (SI) Precipitation
In situ precipitation of scale inhibitor (SI) has been induced deliberately in order to improve its retention characteristics within the porous medium, i.e., the reservoir rock matrix (Malandrino et al. 1995), and this has been modeled previously (Yuan et al. 1994b;Malandrino et al. 1995) both in the laboratory and in the field. A complete consistent model for kinetic coupled adsorption/precipitation has been presented previously (Sorbie 2010). However, here we will take a much simpler approach in order to clarify what is happening in kinetic precipitation processes in 1D rock cores (at the laboratory scale) in order to obtain some useful working formulae to describe the process.

Bulk Precipitation
The simple bulk precipitation model assumed is that which has been used in the software codes for many years (Malandrino et al. 1995). The solid precipitate forms to some level, , and there a corresponding solubility (C s ) in the bulk liquid phase is reached; both C s and are expressed in molar (M) units. Note that if further precipitate is added ( is increased), then C s would remain constant, so long as there is precipitate present. The kinetic aspect of the precipitation/dissolution process is really relevant only in the dissolution part of this model, and here we describe this process by the following simple rate law: where C = C(t) is the mobile phase SI concentration, C s is the solubility of the precipitate as explained above, and κ is the dissolution rate constant. This simple kinetic model of precipitation in a bulk system given by Eq.
(1) has the well-known analytical solution:

Kinetic Precipitation in a 1D Linear Core System
We now consider taking our simple bulk kinetic model and putting it into a transport model for flow through a 1D porous medium, as may occur in flow through a packed bed or in a core flood, as shown schematically in Fig. 1. For a volumetric flow rate Q (say in cm 3 /s) through a linear core of length L (cm), cross-sectional area A (cm 2 ) and porosity φ, the fluid residence time and the fluid velocity v (cm/s) are related as follows: The sequence of steps in an experiment involving a simple precipitating system can be imagined as shown in Fig. 2. The key problem is to calculate what is observed in stage 5. At the initial condition (t = 0) defined by stage 4 in this figure, it is clear that C = C s all along the core; i.e., there is precipitate distributed all along the system and the mobile phase concentration is at the solubility of this precipitate, C s . Similarly, there is a uniformly distributed finite amount of scale inhibitor precipitate present in the system, denoted 0 . The task is to calculate the evolution of the concentration profile as the displacement progresses and we inject only brine containing no scale inhibitor (C = 0).
Assuming no dispersion takes place, the mobile phase SI concentration C = C(x, t) and the amount of SI precipitate = (x, t) can be coupled together in the following system of partial differential equations, using the bulk kinetic model in Eq. (1): The initial conditions for this system are C (x, 0) = C s , (x, 0) = 0 on the interval 0 ≤ x ≤ L representing the 1D porous medium. The boundary condition C (0, t) = 0 for t > 0 is needed as well, expressing the assumption that that the concentration at the reservoir inlet drops to zero immediately after the flow is started and remains so for all time. From Eq. (5) we deduce that, as longs as > 0, the rate of dissolution of the precipitate at x = 0 is the constant −κC s . The precipitate level here decreases linearly as (0, t) = 0 − κC s t. We see that (0, t) = 0 at time t * given by It will be shown that the solution satisfies (x, t) > 0 on (0, L] for all 0 ≤ t ≤ t * , from which it follows that the precipitate runs out for the first time at x = 0 (when t = t * ). The Heaviside term H ( ) in Eq. (5) then implies a sudden change in the dissolution rate from κC s to 0. We will explore what effect this has on the solution in the region of where t > t * . The key idea will be the invariance of a certain relationship between C and . This will lead to the emergence of a traveling wave solution in the t > t * region, joining on to the existing solution profiles of C and . It will be shown that the resulting composite solution of the system obeys the principle of mass conservation, thus confirming that our conjecture is physically correct. This is also verified by comparison with numerical solutions of the initial / boundary value problem using a first-order finite difference scheme.
Another approach to solving the Cauchy problem for Eqs. (4) and (5) would be to notice that it is a hyperbolic system of two partial differential equations. It is already in canonical form, and the distinct eigenvalues 0 and v lead to two families of characteristic curves on which the solutions can be developed (Rhee et al. 1989). However, this is complicated because of the non-homogeneity of the system and the difficulty in dealing with the function H . Even if we were to replace H with one of its smooth approximations (Bracewell 2000), the resulting pair of ordinary differential equations is hard (if not impossible) to solve explicitly. In contrast, the solution method presented in this paper is very straightforward and intuitive.

Solution for 0 ≤ t < t *
We first consider the period when there is some precipitate at every location in the reservoir. At t = 0, we have = 0 for all 0 ≤ x ≤ L. Since ∂ /∂t is finite, there must be some positive time interval 0, t such that (x, t) > 0 for all 0 ≤ x ≤ L, t ∈ 0, t . This

Method of Characteristics
The solution strategy used to solve Eq. (7) is the method of characteristics (Alinhac 2009), which we will outline here briefly. Consider a general quasi-linear PDE with initial conditions on a domain D ⊂ R 2 : We may visualize a solution η = η (x, t) of Eq. (8) as a surface in R 3 : Now consider the vector field V : D × R → R 3 given by The key observation is to note that To see this, it suffices to note that at any point Q ∈ S, the vector normal to the surface is given by Next, let us consider the integral curves of V . In terms of a parameter s ∈ R, these are the curves γ = γ (s) that satisfy the ordinary differential equation (ODE) A fundamental result in the theory of ODEs states that for a given point (x 0 , t 0 , z 0 ) ∈ R 3 , there exists δ > 0 and a unique solution γ : Hirsch et al. 2004 for example). Thus, for any point Q ∈ S there is an integral curve γ of V whose derivative at Q is γ (0) = V Q , tangent to S. Consequently, the integral curve starting at Q remains on the surface S. (This can be shown by choosing suitable local coordinates.) We now reverse these ideas. The initial data η (x, 0) = η 0 (x) in Eq. (8) can be thought of as a curve 0 ⊂ R 3 . From what was said above, for each point Q 0 ∈ 0 we can construct an integral curve of V starting there. By taking the union of this and all other integral curves starting in a neighborhood of Q 0 ∈ 0 , we actually obtain a surface S ⊂ R 3 made up of integral curves. (This process is referred to as "weaving" in Alinhac 2009.) Clearly, S constructed in this way contains 0 and is also such that V is tangent to it at any point. Finally, it can be shown that, at least locally, S is the graph of a smooth function.
In practice, the method for solving Eq. (8) works as follows: using a parameter r ∈ R, we can describe the initial data curve 0 and look for the surface S = {(x(r, s), t (r, s), z(r, s)) , (r, s) ∈ } obtained by solving the system of differential equa- The constants of integration are functions of r only and are determined using 0 .

Analytical Solution of the System When > 0
Comparing Eq. (7) with Eq. (8) we have a (x, t, z) = v, b (x, t, z) = 1, c (x, t, z) = κ (C s − C) and we therefore need to solve the system In terms of arbitrary functions f , g, h, these equations yield x (r, s) = vs + g (r ) Let us choose h ≡ 0 and g(r ) = r , so that Finally, the general solution to Eq. (7) is The problem now is to specify the function f for the situation we wish to model. From Eq. (18), C (x, 0) = C s + f (x). In order to satisfy the initial condition C (x, 0) = C s on [0, L], we therefore require f (x) = 0 for x ≥ 0. We also have the boundary condition The solution C is obtained by substituting Eq. (19) into (18). Thus, This concentration profile is characterized by an upward shock front of height C s exp −κ x/v moving along [0, L] at the fluid velocity v. Behind the front the solution profile is equal to the steady-state solution of Eq. (7) obtained when setting ∂C/∂t = 0.
In order to find a corresponding precipitate profile, we substitute Eq. (20) in Eq. (5) and obtain ∂ /∂t = 0 for all x ≥ vt, which means that = 0 here. For x < vt, we have ∂ /∂t = −κC s exp(−κ x/v). Integration yields the solution in terms of an arbitrary function: We will choose g in such a way that the profile = (x, t) is continuous everywhere, which implies that Eq. (21) should take the value 0 at x = vt. Substituting t = x/v, this condition may be expressed as This choice of the function g results in the following precipitate profile: Fig. 3 Equations (20) and (23) describe solution profiles for t < t * Differentiating Eq. (23) with respect to x, we obtain This shows that the first time the precipitate runs out must be at x = 0, t = t * = 0 /κC s . Until this time, the concentration and precipitate profiles on 0 ≤ x ≤ L are entirely defined by Eqs. (20) and (23). Figure 3 shows these profiles at some time t < t * . The horizontal line represents a fixed precipitate value P ∈ [ 0 − κC s t , 0 ]. Let x P denote the point of intersection of (x, t) with the line = P, so Using the Lambert W-function (Corless et al. 1993), Eq. (25) can be solved explicitly for x P as a function of t. However, for reasons that will become clear shortly, we are mainly interested in the rate of change of x P . By defining G (x P , t) = (x P , t) − P and applying the implicit function theorem (Krantz and Parks 2012) to the equation G (x P , t) = 0, we find Since P is a constant, ∂G/∂t = ∂ /∂t and ∂G/∂ x P = ∂ /∂x P . Introducing the notation C P = C (x P , t) and using Eqs. (24) and (25), we obtain Substituting these into Eq. (26) yields the rate of change of x P in terms of the precipitate level P and the corresponding concentration value C P : 123 Analytical Solutions of a One-Dimensional Linear Model …

Fig. 4 Solution profiles at t = t *
For instance, at some fixed time τ < t * the precipitate level at x = 0 is P = 0 − κC s τ with corresponding concentration level C P = 0. Equation (29) now implies that the "horizontal" rate of change in the precipitate profile at the inlet is v/(1 + κτ ).

Solution for t ≥ t *
The precipitate profile at time t * satisfies (0, t * ) = 0 and (x, t * ) > 0 for x > 0. The dissolution rate at the reservoir inlet suddenly dropped from κC s for t < t * to 0 at t = t * , and the system of Eqs. (4) and (5) changed to its alternate form ∂C/∂t = −v∂C/∂ x, ∂ /∂t = 0 here. This simpler system is no longer satisfied by Eqs. (20) and (23). In order to describe what is going on, we introduce the point (α , 0) on the precipitate profile such that ≡ 0 if x ≤ α and > 0 if x > α . Similarly, we let (α C , 0) be the point on the concentration profile such that C ≡ 0 if x ≤ α C and C > 0 if x > α C . These definitions for the "roots" are physically motivated: precipitate and concentration cannot be negative quantities and, once depleted in a particular location by the flow, they cannot reappear there without external influence. We also observe that we must have α C ≤ α . Indeed, if α C > α , then C ≡ 0 and > 0 on the interval α < x ≤ α C . But since C ≡ 0 is not a solution of Eq. (4) when > 0, this is a contradiction. Figure 4 shows the situation at t = t * , when α = α C = 0. On 0 < x ≤ L, Eqs. (20) and (23) still define the solution profiles and therefore Eq. (29) still describes the speed x P (t) of any precipitate value 0 < P ≤ 0 in terms of the corresponding concentration value C P = C (x P , t). Let us now make the assumption that this intrinsic relationship also applies to the motion of the point (α , 0), despite the fact that the system of equations takes the form ∂C/∂t = −v∂C/∂ x, ∂ /∂t = 0 at x = α . Then, substituting P = 0 and C 0 = C (α , t), we obtain the relation Because 0 ≤ C 0 ≤ C s and 0 > 0, Eq. (30) implies that α (t) < v.
In order to find the rate of change of α C , we observe that any solution of the transport equation ∂C/∂t = −v∂C/∂ x must be of the form C = C (x − vt). Since such a solution is propagated forward at velocity v, we deduce that α C (t) = v. Given that α (t * ) = α C (t * ) = 0, it would therefore seem that α C is about to overtake α . However, as noted above, there is the constraint α C ≤ α . This condition instantaneously "slows down" α C , forcing it to have the same speed as α , from which it follows that α C = α during some time interval t * , t * + δt . Repeating this argument at t = t * + δt and subsequent stages, we deduce that α C = α for all t ≥ t * .
In line with the notation used previously, we let x 0 := α , the root of the precipitate profile. Since x 0 = α C by the above argument, we find that C 0 = C (x 0 , t) = 0 in Eq. (30), and hence dx 0 dt Using the initial condition x 0 (t * ) = 0, we find x 0 (t) = U (t − t * ). We now look for a solution C * of the system of Eqs. (4) and (5) satisfying the condition C * (x 0 , t) = 0 for all t ≥ t * . Since > 0 for x > x 0 , C * is of the form of Eq. (18). Substituting x = x 0 (t), the "moving boundary condition" is Introducing the variable z = U t − U t * − vt and rearranging terms, we find Putting C * into Eq. (5), we obtain As before, integrating with respect to t determines the precipitate up to an arbitrary function h = h(x): We can use the condition * (x 0 , t) = 0 to find h(x) = 0 . Comparing this with Eq. (34), we see that * is just a scalar multiple of C * : Direct substitution shows that C * and * coincide with Eqs. (20) and (23) . It is therefore reasonable to suggest that the solution be defined by Eqs. (34) and (37) in the "new region" U (t − t * ) < x < v (t − t * ) and by Eqs. (20) and (23) for all x ≥ v(t − t * ). This idea is illustrated in Fig. 5. It is consistent with the requirement that any discontinuities in the solution or its derivatives must propagate along characteristics (see Rhee et al. 1989). In Sect. 6, we will explain why this requirement does not apply to x = x 0 (t).
We note that C * , * satisfy ∂C/∂t = −U ∂C/∂ x, ∂ /∂t = −U ∂ /∂x, respectively, and that they are consistent with Eq. (29): substituting P = * (x P , t) and C P = C * (x P , t) = C s P/ 0 yields dx P /dt = U . We also observe that the velocity of a precipitate value P changes in a continuous way as it enters the new region. To see this, note that the front end of this region reaches P at time τ = t * + x P /v. Using Eq. (23), we then have P = Fig. 5 New solution region at some time t > t * 0 − κC s t * exp (κt * − κτ ) and C P = C s − C s exp (κt * − κτ ). Substituting these values in Eq. (29) yields dx P /dt = U , which shows that P already has velocity U when it falls into the new solution region. On the other hand, the velocity of the corresponding concentration value C P "jumps" from 0 to U at t = τ .
Dropping the notation involving C * and * , the solutions for t ≥ 0, 0 ≤ x ≤ L are concisely written as follows: and It remains to check that these composite solutions conserve the total amount of scale inhibitor initially present in the system, which is L (C s + 0 ). The amount of chemical that is added to the mobile phase after some time t > 0 is found by evaluation of the definite integral of C(x, t) between x = 0 and x = vt. For the solution to be physically corrected, this must be equal to the total amount of precipitate used up. In other words, we need It is straightforward to show that Eqs. (38) and (39) indeed satisfy Eq. (40). An alternative way to prove the conservation of mass is to consider the effluent concentration flux vC eff (t) := vC (L , t), which is the amount of chemical passing the reservoir outlet x = L per unit time. Using Eq. (38) we have Evaluation of the definite integral of vC eff (t) between 0 and t * + L/U yields The question that arises is whether the argument can be reversed and if we can arrive at Eq. (38) by imposing the mass balance requirement in the first place. If we knew a priori that U is constant, this task is simple. We would start out by assuming the existence of traveling wave solutions C (x − U t) and (x − U t) on the interval U (t − t * ) < x < v (t − t * ) and then use Eqs. (20), (23) and (40) to derive the correct value of U . The trouble is that, in general, U is not constant, and traveling wave solutions do not exist (see Stamatiou and Sorbie 2018a, b). We would therefore require a more general derivation which proves that the "invariance" of Eq. (29) is a direct consequence of mass balance.

Interpretation of the Solution in Terms of Characteristics
We will now describe Eq. (38) in terms of characteristics.
The Cauchy data on ∂ are C (x, 0) = C s and C (0, t) = 0. If t < t * , Eq. (43) reduces to a single, non-homogenous equation on with general solution given by Eq.  20), thus defining initial conditions on ∂ 1 and ∂ 2 . We also have the boundary condition C (0, t) = 0 for t > 0, which allows the solution on 1 to be determined completely. In order to find the solution on 2 , we need to have knowledge of C on the boundary curve x 0 = x 0 (t). This is illustrated in Fig. 6. The characteristics of Eq. (43) are the lines x − vt = β. Since the equation is homogenous in 1 , the solution here is constant along characteristics. With C (0, t) = 0, this implies that C ≡ 0 in 1 , which also determines the required boundary condition on ∂ 2 as C (x 0 , t) = 0. Equation (18) can now be specified to obtain Eq. (38), as was done in Sect. 5. The solution profile at some τ ≥ t * consists of four regions.
The boundary data on ∂ 1 are "sent" along characteristics to the boundary curve x 0 = U (t − t * ), thus providing the information needed to solve Eq. (43) on 2 . This is only possible because U < v here. Since this will not be the case in general (see Stamatiou and Sorbie 2018a, b), it is worthwhile examining what would happen if U > v. For simplicity, we will assume that U is a constant, although we could equally well carry out the analysis in this section for a more general velocity, U = U (x, t, C) say.
Suppose that the solution is known at some time τ (see Fig. 7). Let D = (x 0 , τ) ∈

Fig. 7
Moving boundary in the case U > v the moving boundary "builds" a traveling wave solution in 1 from the solution in 2 . As a consequence, C no longer vanishes at x 0 (i.e., α C = α ). Indeed, the process described in Fig. 7 results in a (nonzero) traveling wave falling behind the boundary, so that we have α C < α .

Numerical Example
Just one numerical example is presented here, since this one contains all the regions that are predicted to occur by the analytical solution. We consider the case with full solubility level C s = 1 and initial precipitate level 0 = 0.5 on a reservoir of unit length (L = 1). Let the dissolution rate parameter and fluid velocity be κ = 1 and v = 0.5, respectively, so that the Damkohler number is κ L/v = 2. With these parameters, we find t * = 0.5 and U = 1/3. Figure 8 shows the concentration and precipitate in situ profiles for this case at time t = 1.5. The exact solutions given by Eqs. (38) and (39) are plotted together with two numerical solutions obtained using a first-order finite difference scheme. While the coarse calculation ( x = 0.01) clearly shows the effects of numerical dispersion, the finer calculation ( x = 0.001) achieves very good agreement with the exact solution. Figure 9 shows the effluent concentration flux curve for this example. The point x 0 reaches the outlet at t = t * + L/U = 3.5 and C ≡ 0 afterwards. We note that increasing or decreasing the Damkohler number κ L/v in these plots would lead to a higher or lower concentration profile, respectively. Analogous to this "vertical" effect of the Damkohler number, the ratio

Conclusion
In this paper, we have explored the properties of a simple model for the kinetic precipitation of scale inhibitor. It was assumed that a finite amount of scale inhibitor ( = 0 ) was uniformly distributed along a rock core of length L, containing a mixture of water and scale inhibitor at equilibrium solubility (C = C s ). Fresh water (C = 0) was then injected into the core at a constant flow rate Q. This process can be described by the non-homogenous, coupled system of Eqs. (4) and (5) involving the source term κ (C s − C)· H ( ). The Cauchy problem with C (x, 0) = C s , C (0, t) = 0 and (x, 0) = 0 was solved by construction of a weak solution consisting of four components. The key feature of this solution is a moving boundary x = x 0 (t) dividing the domain into regions 1 (where ≤ 0) and 2 (where > 0). The boundary velocity U was determined using the assumption that the relationship expressed in Eq. (29) is invariant when going from > 0 to ≤ 0. In the particular case studied in this paper, we found that U is constant and less than the fluid velocity v. But this is not true in general. In Stamatiou and Sorbie (2018a, b), we will encounter a more complex example with non-constant boundary velocity U = U (x, t, C) such that U < v in some regions and U > v in others. This will lead to the kind of interaction between the solutions in 1 and 2 we described in Sect. 6 for the case of constant velocities.
It was shown that the weak solution C(x, t), (x, t) is physically correct in the sense that the total amount of scale inhibitor is conserved. This was done by explicit integration of the concentration flux through x = L. Although this approach is certainly sufficient, it begs the question as to whether the correct velocity U can in fact be derived by requiring the solution to satisfy mass conservation. This remains an unsolved question.
An example with specific parameter values was calculated in which near perfect agreement between the numerical and analytical solutions was achieved.
The present case is very interesting in its own right. It examines the non-trivial effect of a discontinuity in one of the equations defining the system. However, the solutions derived in this paper can be considered as the base case solutions to a more complex problem in which the assumption of a single flow rate is relaxed. In laboratory experiments, the flow rate is often changed after a shut-in period during which there is no flow at all. In order to describe this with the present linear model, we must introduce a time-dependent expression for the fluid velocity v. This changes the setup of the problem and, as we shall explore later (Stamatiou and Sorbie 2018a), leads to a rather more complex analytical solution. Furthermore, we will also show that the finite precipitate ( 0 ), variable rate model can also be solved in the presence of an equilibrium nonlinear (Langmuir) adsorption (Stamatiou and Sorbie 2018b).
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.