Model-based characterization of permeability damage control through inhibitor injection under parametric uncertainty

Damage in subsurface formations caused by mineral precipitation decreases the porosity and permeability, eventually reducing the production rate of wells in plants producing oil, gas or geothermal fluids. A possible solution to this problem consists in stopping the production followed by the injection of inhibiting species that slow down the precipitation process. In this work we model inhibitor injection and quantify the impact of a set of model parameters on the outputs of the system. The parameters investigated concern three key factors contributing to the success of the treatment: i) the inhibitor affinity, described by an adsorption Langmuir isotherm, ii) the concentration and time related to the injection and iii) the efficiency of the inhibitor in preventing mineral precipitation. Our simulations are set in a stochastic framework where these inputs are characterized in probabilistic terms. Forward simulations rely on a purpose-built code based on finite differences approximation of the reactive transport setup in radial coordinates. We explore the sensitivity diverse outputs, encompassing the well bottom pressure and space-time scales characterizing the transport of the inhibitor. We find that practically relevant output variables, such as inhibitor lifetime and well bottom pressure, display a diverse response to input uncertainties and display poor mutual dependence. Our results quantify the probability of treatment failure for diverse scenarios of inhibitor-rock affinity. We find that treatment optimization based on single outputs may lead to high failure probability when evaluated in a multi-objective framework. For instance, employing an inhibitor displaying an appropriate lifetime may fail in satisfying criteria set in terms of well-bottom pressure history or injected inhibitor mass.


Introduction
Formation damage is a generic terminology referring to the impairment of the permeability of subsurface reservoirs by various adverse processes [1]. Being an irreversible process, if not treated on time and properly, formation damage prevents the efficient exploitation of reservoirs [1,2].
Porosity decrease and permeability impairment can be the result of many chemical and physical processes including mechanical deformation, swelling of clay minerals, mineral and organic precipitation/dissolution, and thermal deformation [3][4][5][6][7][8]. Formation damage control and remediation are among the most important issues to be faced for maintaining an efficient exploitation of hydrocarbon or geothermal reservoirs [1]. We focus in this work on the damage in the near-well region made by the mineral precipitation of dissolved substances in the formation fluid, as sketched in Fig. 1.
A classical approach to modelling of mineral precipitation is based on the saturation index (SI) [9][10][11][12]. In this approach the precipitation kinetics is proportional to the difference between the concentration of aqueous species and an equilibrium concentration [9,10]. Assuming a simple system where a single chemical species is considered the SI [−] is defined as [12,13] SI ¼ The equilibrium concentration c eq is influenced by the reservoir conditions and, in turn, affects the saturation index. The primary factors that affect the aqueous solubility and contribute to high values of saturation index, reported in literature, are e.g. pH [14], temperature [5,9,14], velocity [5] and pressure change [5,7,12,15], mixing with an incompatible fluid [1,7,13,14,16,17]. When the saturation index exceeds the value of one, the minerals precipitate, as a result, eventually accumulating on the rock-solid matrix or to the well structure. Under natural conditions the fluid in the reservoir is in equilibrium with the host rock: the saturation index is nearly neutral, namely, precipitates are not expected to form, and deposits are not expected to dissolve. Drilling of boreholes to exploit the subsurface fluid, injection to and extraction from the wells cause fluid displacement and consequently a departure from the equilibrium [16][17][18].
To prevent the occurrence of formation damage, inhibitors have been used successfully in conventional geothermal water production [7,[19][20][21][22]. The most widely used technique to deliver the scale inhibitors to the subsurface system is an injection treatment [1,21,23], during which the production stops and the well is treated like an injection wellthe fluid (compatible with the reservoir one) is continuously injected for a fixed amount of time from the well-bottom to the reservoir together with the dissolved inhibitor. A part of the inhibitor adsorbs on the solid rock surface preventing precipitation with various degrees of efficiency [18,21]. Thus, the inhibitor must have some affinity with the host rock to be adsorbed on its surface and modify the kinetics of mineral precipitation dissolution [7,18,19]. However, the precise mechanism by which the inhibitor affects this process is not understood clearly [1,18].
Our key objective is to consider the joint effect of the variability (or uncertainty) of diverse input properties on the effectiveness of the inhibitor injection treatment. We study in detail the effects of uncertain inputs falling into three distinct categories: i) the rock-inhibitor affinity, described through a Langmuir isotherm, ii) inhibitor efficiency quantifying the effects of the inhibitor on the mineral precipitation kinetics, iii) treatment design properties characterized by the inhibitor concentration and injection time, i.e., the duration of the injection phase of the treatment. To this end, we formulate a mathematical model for mineral precipitation and inhibition treatment by injection as a system of partial differential equations and constitutive laws, based on previous literature [1,2,8,[10][11][12][24][25][26][27][28]. Output of the model is a set of unknown variables such as pressure, permeability, fluid velocity and inhibitor concentrations. We then build a numerical simulator to obtain an approximation of these variables and ultimately increase our understanding of the involved processes. Our simulator relies on a set of simplifying assumptions, yet it retains coupling between fluid flow and transport processes and allows providing a characterization of uncertainty propagation in the system. As a first key objective of our work we diagnose the system response to these uncertain inputs through a global sensitivity analysis [29,30]. The considered problem entails nonlinear relationships, thus robust global sensitivity indicators are required to parameterize the input-output mapping. In this work we rely on momentbased sensitivity indices developed in [29,30]. As a second objective, we characterize the impact of input uncertainty onto various indicators of the performance of the treatment process. Field characterization of inhibitor characteristics is often based on return concentration histories observed by injection and subsequent extraction of the inhibitor through a single well [31], where optimization criteria are often set in terms of inhibitor mass and of required well operations. As the inhibitor has to be injected through the extraction well, each injection has a definite cost that should be minimized. Our specific interest is to identify existing relationships between target output variables, notably between the inhibitor lifetime observed in a return concentration breakthrough curve, pressure values observed at the well after a fixed time and the total mass of injected inhibitor within a given time window. The objective of our analysis is quantifying the mutual information between these outputs, i.e., numerically quantifying the relationship existing between these different performance indicators and model parameters. In this context our work allows quantifying failure and success probabilities under different scenarios, set in terms of inhibitor properties and treatment design. These efforts could (a) lead to a better understanding of the process, and (b) be used as a scientific guidance for development of remediation strategies for formation damage control in reservoirs.
The work is organized as follows. Section 2 presents the assumed mechanistic mathematical model and related assumptions (Sec. 2.1-2.2). Then we present the simulation setup and the related input (Sec. 2.3), the output variables of interest along with the employed sensitivity indicators (Sec. 2.4). Section 3 present the results, by first recounting the results obtained from the sensitivity analysis and then projecting those to analyse the impact of input uncertainty on a set of three performance indicators. Concluding remarks end the study.

Model definition and assumptions
A deterministic mathematical model of a reservoir generally requires solving a system of partial differential equations (PDEs). The number and types of the equations depends on the dominant governing processes. Our objective is to provide a simplified description of the physical and chemical processes, while preserving a reasonable applicability level, particularly as for what concern the feedback between inhibitor transport and well-bottom pressure. The main model concepts employed here are derived by former modelling studies dealing with solute transport in groundwaters and porous media, readopted to the specific problem by following the general theory of formation damage [1]. The model and the assumptions are described next.
We assume a radial domain r ∈ [r well , R] of radius R [m] where r well [m] is the radius of the well. We consider here a single well model, in other words, either there exists only one well in the field, or other wells existing in the same exploitation field are assumed to be far enough that the discharge and the drawdown of each well will not be affected by neighbouring wells. This assumption allows neglecting any interference to the dynamics of the model due to neighbouring wells. The height H [m] of the well screen is much smaller than the radius R [m] of the circular area under consideration, H ≪ R. The existence of radial symmetry allows considering a onedimensional radial flow on a circular plate in the subsurface as a depth-averaged model. This is a fundamental simplification of the problem, even if widely accepted in literature (see [1,25,26], for instance). Assuming full penetration of the well screen allows expecting that the well receives fluid from horizontal flow. From a physical standpoint we assume single-phase fully saturated flow and isothermal conditions, assuming the temperature fluctuations in deep subsurface reservoirs can be neglected. The reservoir properties are assumed to be homogeneous thus neglecting spatial fluctuations of porosity and permeability. In this model we consider a single precipitating substance and thus the selected constitutive laws (e.g., precipitation model) is linked to the kinetics of this matter as defined by the saturation index (Eq. (1)). In our model only precipitation is considered; dissolution of the species from the formation rock is not considered, i.e., we assume that the solution is always supersaturated in the near well region. We neglect the influence of the chemical substance on the fluid properties.
The adopted simplifications of the system physics and geochemistry are in line with key target of the work, i.e., the analysis of pressure and inhibitor concentration dynamics from a fundamental standpoint. For instance, the extension to non-isothermal flow may be envisaged to address the variability in fluid rock interactions due to temperature variations, which for example could be triggered by the injection of a fluid having a different temperature with respect to the reservoir [32,33]. Moreover, detailed geochemical rock-fluids interactions and the heterogeneity of the system properties are here neglected, as addressing these elements would necessarily introduce additional levels on uncertainty, thus preventing an in-depth analysis of the feedback between the different processes under consideration.

System of PDEs and constitutive Laws
The system is set in a one-dimensional spatial coordinate r ∈ When mineral precipitation occurs, effective porosity of the medium ϕ(r, t) [−] reduces [1,12,[26][27][28]: where ε p is the porosity difference (fractional bulk volume) due to precipitation.
Fluid mass and momentum balance. In cylindrical coordinates (r, θ, z) and considering one dimensional radial flow, the mass balance equation reads as [12].
Transport equation for dissolved matter. Our model considers that precipitation can be computed through a single chemical species. The advection-dispersion-reaction equation for the concentration of dissolved matter, c p [mol/m 3 ], reads as [12,25].
where v = u/ϕ [m/s] is the average velocity, a L [m] is the longitudinal dispersivity coefficient and the right-hand-side represents the sink term which is the precipitation rate in the pore spaces of the reservoir. The precipitation rate, R p [mol/m 3 /s] is defined as [5,10,[12][13][14]: R p is proportional to the reaction constant k p [mol/m 2 /s] and the specific surface area S [m 2 /m 3 ] of the pore space. Λ represents the saturation index defined as in Eq. (1). The exponent m takes various values for different precipitation regimes [13]. Note that concentration c p can be interpreted as a generic concentration of a dissolved mineral. Our model greatly simplifies mineral-water geochemistry with respect to a real scenario and neglects various environmental indicators which may influence the process, such as fluid chemical composition, mineralogical compositions, temperature, and pressure. This is justified by our objective to demonstrate the coupling between the different considered processes (fluid flow, inhibitor transport, mineral precipitation) in a simple parameterization setup.
The time-rate of change of porosity is proportional to the precipitation rate: where V s [m 3 /mol] is the molar volume of the considered mineral, here taken to be pure calcite for illustrative purposes. Permeability reduction can be represented as an exponential function of the fractional bulk volume (porosity difference) ε p [1,35]: where k 0 [m 2 ] is the initial effective permeability of porous medium and a [−] is an empirical coefficient that takes different values according to the soil type [1]. The reason to choose an exponential relationship is that the pore throat clogging can cause more permeability damage than pore surface deposition and decrease permeability to negligible values even in the presence of residual porosity [1,3].
Transport equation for inhibitor. Inhibitors act by establishing a chemical interaction with the rock surface, which can be modelled as a sorption-desorption process. Various kinetic and equilibrium models can be employed to address such chemical processes [25]. We consider here a nonlinear equilibrium Langmuir adsorption-desorption isotherm [18,25] to describe the inhibitor-rock affinity: where c i [kg/m 3 ] denotes inhibitor concentration in the injection fluid, b [m 3 /kg] is the inhibitor adsorption energy coefficient, F max [mg/g] is maximum inhibitor adsorption capacity and F is the mass of the adsorbed inhibitor per unit mass of the solid. Parameter F max indicates the maximum amount of solute that can be adsorbed on the solid rock per unit mass of the rock. The adsorption model in Eq. (9) is here selected as it is widely employed in the literature to handle sorptiondesorption processes in porous media, see e.g. [18,32,36]. A further motivation supporting our choice is that this model accounts for the maximum adsorption capacity, as opposed to a linear equilibrium sorption model, which, in principle, allows infinite mass to be adsorbed to a finite volume of rock mass. This feature allows testing this property of the system in our sensitivity analysis (see Section 4). Transport equation for the inhibitor concentration expressed in one-dimensional radial coordinates reads as [25].
where ρ b [kg/m 3 ] is the bulk density of the solid. The effect of the inhibitor on the precipitation is modelled as the reduction of the precipitation rate, R p defined in Eq. (6): where R pi denotes the precipitation rate under influence of the inhibitor. We define η as inhibitor efficiency coefficient and n as inhibitor efficiency exponent, both parameters being dimensionless.
Hence the time-rate of change of porosity Eq. (7) is modified as In summary, Eqs. (3-12) constitute the system of PDEs for the whole process.
Initial and boundary conditions. To close the system of PDEs initial and boundary conditions are required. In our model production and injection rate Q [m 3 /s] is fixed which implies a fixed radial gradient for pressure. At r = R Dirichlet boundary condition is imposed for the pressure, i.e., p = P R . The steady-state solution of the mass balance equation is taken as the initial condition.
Homogeneous Neumann condition is set for both c p (during production and injection) and c i (during production) at the well. A fixed inhibitor concentration c i = c well is set at the well during the injection phase. At r = R homogeneous Dirichlet condition c i (r = R) = 0 for the inhibitor and equilibrium concentration c p (r = R) = c eq for the dissolved minerals are imposed.
To ensure a smooth decay of the concentration of the precipitating species and to satisfy the homogenous Neumann boundary condition at the well, a parabolic behaviour is chosen as initial condition. The inhibitor concentration is set to zero for the whole domain at initial time.
Numerical model The numerical solution of the system is obtained through a finite difference (FD) method and the corresponding code is written in Python 3. The radial spatial domain is represented by a non-uniform grid with increasing grid size starting from the well. Discrete derivatives are approximated following [37]. The non-linear terms (arising because of the Langmuir sorption term) in Eq. (10) are dealt with by discretizing with an explicit time scheme. We verified our solver achieves satisfactory results in terms of mass conservation and performed a numerical grid convergence assessment to verify the convergence properties (see Supplementary Materials for details).

Model set-up
The simulation set-up is based on alternate production and inhibitor injection. The total simulation time T is selected upon simulating continuous production and recording the time at which the pressure at the well p well reaches a value equal to 30% of the value set at the initial time, p well, 0 . This simulation is performed with reference values of model parameters and leads to identifying a time window T = 3522 days, which is an estimation of the lifetime of the well in the absence of inhibitor injection (also termed as normal production scenario, in the following).
In the interval t ∈ (0, T] we simulate multiple realizations of production with inhibitor corresponding to the following cycles simulated periodically: 1) Injection of the inhibitor: The inhibitor is injected in the domain with a predefined concentration c well for a fixed injection time T i . 2) Production: When the injection finishes, the well is set back on production. Injection is applied again when the pressure at the well-bottom attains a value less than 95% of the initial pressure of a given production phase, p well, phase .
We perform periodical injection-extraction (production) processes until time T which is the lifetime of the well that is obtained from the simulation of normal production. Through this workflow we can easily assess the impact of the inhibitor treatment on the well pressure p well as compared to the normal production case (Fig. 2).
Production with inhibitor is here simulated considering a set of model parameters as stochastic variables. In our analysis the parameters are all considered to be mutually independent and uniformly distributed within a given interval. The extreme values associated with the probability distribution assigned to these parameters is shown in Table 1. Our choice of parameters allows us considering the joint effects of the treatment design (c well and T i ) together with the chemical characterization of the inhibitor-rock affinity (F max , b) and the inhibitor efficiency (η and n). Our target is then to understand the impact of such parameters on the various problem outputs, introduced in Sec. 2.4.

Model outputs
In this section we identify the key outputs of the numerical simulations performed with the model developed in the previous sections.
In the following graphs we compare the results of the normal production (NP) and production with inhibitor (PI). For an illustrative purpose we rely here on a reference realization of the model parameters provided by CHIMEC S.p.A. as reported in Appendix 3.
A first key simulation output is related to the value attained by the well (well-bottom) pressure p well at time T. Figure 3 depicts the temporal evolution of the pressure (a) and permeability (b) at the well-bottom k well during NP and PI. In the following we refer to a dimensionless pressure output In the example shown in Fig. 3, pressure at the well-bottom in PI simulation at time T (day 3522) is 92 bars, thus P well = 0.67. This factor can be directly linked to the observed permeability decrease (or damage) at the well, i.e., higher value of pressure obtained under production with inhibitor correspond to higher values of permeability as compared to a normal production case. In the specific example shown in Fig.  3b simulation permeability is predicted to decrease by more than 50% during NP, whereas under PI the permeability damage at the well is 33%. Another important factor in the inhibitor injection treatment is the number of injections n inj performed in the simulation time window. In practice, the biggest impediment to periodic injection-production process is to remove the pump, set up the injection devices and put back the pump to its place for the next extraction (production) phase [5]. Depending on the type of the well and the reservoir this procedure sometimes can be very time-consuming and costly. A second set of key model outputs is obtained from the analysis of the space-time profiles of inhibitor concentration. Figure 4a shows the inhibitor concentration profile in the near well region after the first injection phase for the reference case scenario (T i = 0.4 days). In reservoir treatment techniques, the minimum inhibitory concentration (MIC) is the lowest concentration of the chemical that prevents the precipitation [23]. In previous studies the recommended MIC is 10 mg/l (or 10 ppm) for the inhibitor in this reference realization. We observe that (Fig. 4a) the inhibitor tends to attain a value smaller than the recommended MIC for to r > 3 m and we identify the inhibitor influence radius, r inf . The influence radius has a direct influence on the pressure profile obtained, from Fig. 4b we clearly see the inhibitor is more effective very close to the well than at faraway locations and this can be explained by the fact that concentration of the inhibitor is higher close to the well. In fact, the two pressure profiles of NP and PI collapse right after the influence area of the inhibitor (r inf ≈ 3 m, in this realization of the model parameters).
Another important output of the inhibitor treatment design is the inhibitor injection lifetime, inh life . Inhibitor lifetime is characterized as the number of days that the inhibitor concentration is measured above the recommended MIC after the first injection. In Fig. 5 we observe that for the reference case scenario the inhibitor lifetime is around 84 days. The inhibitor lifetime is an important parameter often adopted to constrain optimization approaches, i.e., used as a predictive variable to understand the treatment efficiency. In practice, inh life can be determined from preliminary field tests and therefore can be used to gain insights on the efficacy of the treatment process [18].

Sensitivity analysis
We perform a global sensitivity analysis in the production with inhibitor scenario to quantify the effect of the input uncertainty for the inhibitor.
Sensitivity analysis is based on model inputs assigned according to quasi-Monte Carlo selection [29]. The influence of each uncertain parameter x k on the statistical moments of the pdf of the output f can be quantified by the AMA indices [29,30]. In particular, we rely on the indices AMAE and AMAV: where E and Var indicate that the index measures the importance of x k to the mean (expected value) or variance of a model output, respectively [29]. The AMAP sensitivity index can be used to rank parameters based on their impact on the probability to exceed a defined threshold [30].
where P thr Pr[f(x) > thr] is the unconditional probability that the quantity f(x) exceeds the threshold thr and Pr[f(x)| x k ] indicates the same probability conditional to parameter x k . AMAP complements the available AMA moment-based indices by targeting sensitivity with respect to the exceedance probability rather that the statistical moments of the output pdf and is therefore relevant in the context of risk or performance assessment [30].

Surrogate model
To explore sensitivity and uncertainty propagation we rely on extensive sampling of the parameters space. To this end we leverage a surrogate model of the system mimicking the inputoutput mapping at reduced computational effort. The surrogate model is expressed through a polynomial approximation formulated in terms of model input parameters where the generalized Polynomial Chaos Expansion (gPCE) technique [38,39] is employed. A given model output f(x) is approximated as where x is the vector of N uniformly distributed random input parameters, ψ j are orthonormal multivariate Legendre polynomials, and the number of polynomials terms O p is defined as where D indicates the maximum degree of polynomial approximation with respect to a single parameter. Evaluation of the gPCE coefficients, α j , entails solving the complete model to compute f(x) for several combinations of the uncertain parameters. Coefficients α j are computed through a least square minimization of the truncation error: where N coll denotes the number of collocation points in the parameter space and the entries of vector x i correspond to the combination of parameters used for the i-th simulation.  3 Numerical results and analysis 3.1 Global sensitivity analysis Figure 6 shows the probability distribution function (pdf) histogram of the outputs obtained by sampling the considered parameter space. As a result of the assumed parameter variability, we obtain values of the normalized pressure changing between 0.3 and 1.0. Because lower values indicate a loss of performance of the treatment it is important to understand which parameter drives the occurrence of such low values. Influence radius of the inhibitor mostly occurs around 4.5-7 m. The number of injections over 3522 days ranges mostly between 3 and 7. Finally, the inhibitor lifetime has its peak between 0 and 90 days and exhibits a right tail, with largest values approaching 3000 days (see Fig. 6d). AMAE and AMAV indices are shown in Figs. 7 and 8, respectively. The two indices show consistent parameters ranking for given output. As for what concerns the treatment design, the inhibitor injection time has generally a more marked influence than concentration on mean and variance of all the target outputs, i.e. both AMAE and AMAV indices associated with T i are larger than the ones associated with c well . AMAE indices show that all the input parameters have a similar weight on the expected value of P well (Fig. 7a). The efficiency coefficient η and the efficiency exponent n have negligible influence on the mean and variance of the inhibitor lifetime and influence radius (Figs. 7, 8b, d), the statistics of these two outputs being chiefly influenced by treatment design variables (T i , c well ) and the maximum adsorption capacity F max . The mean and variance of the required number of injections n inj display a very similar sensitivity pattern to those of P well (compare panels a and c of Figs. 7,8), which is consistent with the fact that injection of the inhibitor is driven by the pressure decay, as detailed in Sec. 2.3. Figure 9 displays the AMAP indices for the normalized wellbottom pressure and for the inhibitor lifetime. We apply here the definition provided in Eq. (15) with threshold set equal to 365 days and 0.8 for inh life and P well , respectively. Results are in line with those given by AMAE and AMAV thus rendering two different rankings and showing that each model parameter influences the two outputs to a different degree. Notably the parameters η and n have a relevant influence on the probability that the condition P well > 0.8 is satisfied, having AMAP values of approximately 0.15-0.2 (see Fig. 9a). However, the value assumed by the same parameters η and n is irrelevant for the probability to observe inh life > 365 days, as demonstrated by the corresponding AMAP value approaching 0 (see Fig. 9b). Overall, the results in Fig. 9 suggests that the two thresholds may not be necessarily satisfied under the same parameter combinations, as demonstrated in more detail in Sec. 3.2.

Implications for treatment design
We analyse here the implications of the input-output relationship reported in Sec. 3.1 on constrained treatment design. To this end we consider three key outputs, i.e. the dimensionless well pressure P well , the inhibitor lifetime inh life , and quantity M inh = c well T i n inj Q [kg], corresponding to the total mass of injected inhibitor along the whole simulation. Figure 10 shows the bivariate sample probability distributions between these three outputs, as rendered by a Monte Carlo simulation with 10 5 realizations. These results show that the knowledge of one of these outputs is not necessarily informative on the other two, i.e., the level of mutual correlation between the various outputs is generally poor. This behaviour is quantified in Table 2 where we list two indicators that are commonly employed quantify the mutual dependence between variables, i.e., linear correlation coefficient ρ and the uncertainty coefficient UC. This latter is based on mutual information and can also quantify a nonlinear relationship between variables [40]. Note that the values of UC are smaller than those attained by the linear correlation coefficient consistent with results found in [41]. The spread exhibited by the bivariate distributions shown in Fig. 10 and the results reported in Table 2 indicate that the three outputs are providing a markedly different response to the uncertain inputs and this is consistent with the diverse sensitivity patterns exhibited by the simulation outputs, as discussed in Sec. 3.1, and further demonstrated in the following.
Our aim is now to characterize the parameter spaces resulting from output variables satisfying three constraints which are used as performance indicators, leading to discern treatment success or failure. These constraints are set in terms of the three following criteria: The selected constraints are related to the overall efficiency achieved by the treatment as well as to the need to limit the cost associated with the treatment itself, which is proportional to the mass of injected inhibitor and the number of injections. Note that the thresholds imposed here are arbitrarily chosen for illustrative purposes. In the following we analyse the probability to fulfil these three constraints in three scenarios where we fix the values of the parameters b and F max to different combinations, as reported in Table 3. This choice is motivated by the fact that the parameters describing the affinity between the inhibitor and the rock matrix could be determined independently through laboratory experiments. While these estimates will be in principle affected by uncertainty, for the following analysis we neglect such uncertainty and fix them to deterministic values. Thus, the idea is to assess the impact of assuming complete knowledge of the inhibitor-rock affinity properties on the combinations of the remaining parameters leading to treatment success. The three scenarios analysed feature an increasing affinity between the inhibitor and the host rock from Scenario 1 to 3 (see Table 3). Table 3 presents the probability that the treatment is successful in fulfilling each of the three constraints A, B, C independently and the combination of the three criteria, for the three assumed scenarios. These results are obtained upon generating three independent Monte Carlo simulations of 10 5 realizations each through evaluating the surrogate model described in Sec. 2.4 while keeping fixed b and F max and sampling the remaining four uncertain parameters within the intervals reported in Table 1. An increase of b and F max leads to an increase of the probability to satisfy the three criteria independently as well as jointly. In general, we observe that the constraint applied to the inhibitor lifetime (constraint A) is satisfied with larger probability than the other two. For scenario 1 the most stringent criterion is related to the pressure (constraint C), while for increasing inhibitor-rock affinity the injected mass becomes the most limiting factor, i.e., Pr(B) < Pr(C). Table 3 also reports the conditional probability of failing to satisfy B or C given that A is satisfied, indicated as Pr(~B| A) and Pr(~C| A), respectively. Note that if only lifetime was used to parameterize the system (i.e., only constraint A is satisfied) there would be a large probability of failing to control the predicted mass of injected inhibitor as well as the pressure, i.e., Pr(~B| A) and Pr(~C| A) are both larger than 75% for low affinity (Scenario 1, see Table 3). This result is consistent with the poor level of correlation existing between the inhibitor lifetime and injected mass (see Table 2). Our results suggest that simulation assisted multi-objective optimization is essential to control the three constraints at the same time and that caution should be used in employing a single variable as a proxy of the overall efficacy of the treatment. Results obtained in Table 3 are specific to the selected thresholds, but they indicate a trend which is general for the investigated mathematical formulation of the problem.
The results shown in Table 3 demonstrate that even if the inhibitor-rock affinity was known without uncertainty failure probability would remain high, particularly if the three criteria are jointly considered. Therefore, we now investigate how the remaining uncertain parameters are influencing failure. The uni-and bi-variate marginal probability distributions of the parameters corresponding to the realizations satisfying the three constraints together (i.e., A ∩ B ∩ C) are reported in Figs. 11 and 12, respectively. We observe that for an inhibitor exhibiting relatively low affinity with the rock (Scenario 1, Fig. 11a-d), successful treatment is obtained only for parameter values lying close to the upper or lower bounds of the investigated intervals. In particular, joint occurrence of high values of injection time and inhibitor efficiency coefficient η appears to be essential for the treatment success (Fig. 12a). Note that while injection time is a variable subject to an engineering choice, the efficiency parameters η and n may be hard to estimate a priori in field settings. Therefore, in these conditions the treatment is prone to a considerable failure risk, i.e., any uncertainty related to the estimation of η and n may lead to a failure to satisfy one or more of the target optimization constraints. The spread associated with marginal distributions of model parameters increases with the inhibitor-rock affinity (compare Fig. 11i-l with ad). This result has two consequences: i) with increasing affinity more freedom is allowed in selecting the treatment design, given by the combination of injection time and concentration (see Fig. 12b, c), ii) the selected constraints are satisfied in a wider range of the parameters η and n. As a consequence, failure of the treatment due to uncertain estimation of the inhibitor efficiency parameters (η and n) becomes less likely. Finally, we observe that bivariate distributions associated with medium and high affinity scenarios feature well identified regions of the parameter space with positive success probability (Fig. 12b, c). These latter suggest the existence of highdimensional Pareto-like fronts, which could be exploited in an optimization framework. The shape assumed by these fronts is markedly influenced by the level of affinity assumed between the inhibitor and the rock matrix.

Conclusions
We model injection of inhibitor to control permeability formation damage due to mineral precipitation in the near-well region of subsurface reservoirs. Inhibitors hamper the precipitation process by influencing its kinetics. In this work we first develop a numerical simulation tool that considers the numerical solution of the problem in a homogeneous radial domain starting from first principles, upon assuming a carbonate rock domain. Detailed analysis of the fluid-rock geochemical interactions, and thus generalization to different mineral compositions, may be considered future investigations. This forward solver is then used to build a reduced complexity model, which is employed to fully explore propagation uncertainty from input to output variables.
We find that the inhibitor lifetime is largely affected by the rock-inhibitor affinity but is poorly sensitive to the parameters describing the inhibitor efficiency. In contrast well-bottom pressure is chiefly governed by inhibitor efficiency and only mildly sensitive to affinity parameters. Treatment design appears to affect all the considered outputs, the most relevant  factor being the injection time rather than concentration. Diverse sensitivity exhibited by the various quantities is likely related to the nonlinear nature of the investigated mathematical problem. Results of sensitivity analyses like the one we performed can assist practical investigations as they indicate which parameters should be further investigated to improve the control on the process. The nonlinearity of the input-output mapping also affects the impact of parameters on the probability to exceed a selected threshold, which becomes relevant with a view to a treatment optimization framework. Here, the predicted performance is evaluated in terms of three constraints; final well-bottom pressure, inhibitor lifetime and the total mass of injected inhibitor. Despite all variables being coupled in the adopted mathematical formulation, these outputs display a poor level of mutual dependence as quantified by linear and nonlinear indicators. This result suggests that multi-objective optimization frameworks are needed to constrain the treatment  Table 3 Values assigned to parameters F max , b and related sample probabilities of satisfying the three optimization constraint A,B,C individually and jointly, upon considering the uncertainty bounds assumed in Table 1 for the other parameters. The last two colums indicate conditional probability of failure to satisfy B and C when A is satisfied  properties under multiple sources of uncertainty. We project these results to the quantification of treatment success under various levels of rock-inhibitor affinity. As expected, inhibitors displaying a high affinity with the rock allow more freedom in the choice of treatment design and reduce the level of failure probability. For high affinity we find that treatment success is generally maximal under injection times larger than one day. Success probabilities are confined in well-defined regions of the parameter space thus suggesting the existence of high-dimensional Pareto fronts which could be exploited in an optimization design. Finally, we investigate the impact of rock-inhibitor affinity on the probability of failure under uncertainty characterizing the efficiency parameters, which may be hard to constrain in a field setting. The spread associated with marginal distributions of model parameters leading to a successful treatment is chiefly dependent on the inhibitor-rock affinity. Therefore, accurate estimation of the properties of the inhibitor are essential to infer failure probability under uncertain efficiency. Figure 13 shows the prediction through the surrogate model versus the original 10 3 simulation data. We observe that the prediction of the surrogate model accurately reproduces the data. Figure 13c displays the prediction obtained for the number of injections, this latter being a discrete sample the prediction of the injection number is rounded to the nearest integer. The inaccuracy for the influence radius (Fig. 13b) is due to the fact that the simulation domain is discrete, while the prediction is continuous.    Funding Open access funding provided by Politecnico di Milano within the CRUI-CARE Agreement.

Appendix 1
Data availability Numerical data are available upon reasonable request to the corresponding author.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.