© E. Tohidi & S. Sadeghnejad, published by IFP Energies nouvelles, 2021

Nomenclature

: Molality of injected H⁺, mol/g

a_i : The activity of species i, mole/vol. of solution

a : The hydrated ion size

A_l : The reactive surface area of the lth mineral, m²/m³

A_l0 : Initial the reactive surface area of the lth mineral, m²/m³

A_x : Cross-sectional area, cm²

A_j and A_i: Secondary and primary species

A : Solvent based constant

B : Solvent based constant

b : The decrease in solvent concentration

B_w@Res.C. : Water FVF at reservoir condition of P and T, bbl/bbl

T : Temperature, °F

T_ref : Reference Temperature, °K

p : Pressure, psi

p_out : Pressure at the outlet of the core, psi

C_i : Mass concentration of component i, mole/vol. of solution

C_inj : Injection mass concentration of component i, mole/vol. of solution

C₀ : Initial mass concentration, mole/vol. of solution

D : Diffusion coefficient, cm²/s

E_l : Activity energy of the lth mineral, J/mol

I : The ionic strength, mole/vol. of solution

K : Permeability, md

K_eql : Equilibrium constant of the lth mineral reaction

k_l : Rate constant (mol/cm²/s)

k_l ref : Reaction rate constant of the lth mineral at reference temperature of T_ref

L : Core length, cm

Mw_l: Molar mass, g/mol

N : Total number of aqueous species

N_Da : Damköhler number

N_c/N_x: Number of independent/dependent components

N_rs : Number of dissolution/precipitation kinetic reactions

q_m : Production or injection rate, cm³/s

Q_s : Ion activity product

R : Universal gas constants, J/mol/K

R_i : Rate of generation or depletion of species i in solution per unit volume, mol/cm³/s

: Heterogeneous reaction rate, mol/cm³/s

: Maximum H⁺ reaction rate, mol/cm³/s

S : Salinity, percent

t : Elapsed time, s

u_x : Velocity of the bulk liquid solution in the x-direction, cm/s

u_inj : Injection velocity, cm/s

U_i : Total concentration of component I, mole/vol. of solution

V_b : Bulk volume of the grid, cm³

VF_l: Fraction of the lth mineral volume per grain volume

VF_l0: Fraction of the lth mineral volume per grain volume

X_j : Concentration of primary species, mole/vol. of solution

∆x: Difference along x direction, cm

z : Ionic charge

1 Introduction

Waterflooding, as a secondary oil recovery method, is used to improve sweep efficiency in oil reservoirs. However, complicated geochemical reactions can occur due to the incompatibility of injected water with the indigenous pore fluid. Various key parameters influence the resulted non-equilibrium state, such as the type and number of the existing components within the system. Therefore, injected fluid composition and overall salinity are of prime importance to determine scale precipitations as an undesired phenomenon during waterflooding.

Various side effects are expected due to scale formations, including production decline, excess work over, and facilities expenses (Abouie et al., 2017). The formation damage caused by waterflooding operation was reported for oil fields worldwide (Elmorsey, 2013; Mitchell et al., 1980; Moghadasi et al., 2004; Qing et al., 2002; Shuler and Jenkins, 1991). Fu et al. (2012) developed a geochemical model to quantify the reactions occurred among rock minerals, indigenous pore fluid, and injected sea water into a sandstone reservoir. They detected considerable differences between the composition of produced water and simple mixing of seawater and formation water. The observed discrepancies were assumed as a result of geochemical reactions influenced by the fraction of seawater admixed. Lu et al. (2017) conducted an integrated experimental and simulation study to investigate the geochemical characteristic of a sandstone formation rock and four injected brine samples with salinity ranging from 10 000 up to 1 000 000 ppm. The results showed negligible rock porosity and permeability changes. Engineered water injection was also a promising technique to enhance oil recovery with logical cost and limited environmental effects. Claes et al. (2018) conducted experimental and numerical study to characterize the pore structure of sandstone reservoirs. The alteration of flow properties including porosity and permeability caused by diagenetic phases was evaluated in samples with different calcite cementation grad, indicating specific dissolution patterns. Ahusborde et al. (2018) developed a comprehensive model coupling geochemical reactions with multiphase multicomponent fluid flow. This model was proved to be applicable in 2D and 3D simulation of highly heterogeneous fields. The validity of modeling results was checked via existing simulations in the literature. Rostami et al. (2019) proposed a new correlation to describe permeability changes caused by formation damage using a large databank. System parameters including injection rate, pressure, temperature, and fluid composition were employed to predict the permeability damaged. Al-Saedi et al. (2020) implemented the CrunchFlow software to investigate the impact of injected water overall salinity on oil recovery in sandstone reservoirs. They conducted a number of core flooding experiments including the injection of formation water and several diluted water samples. The experimental results were used to calibrate the developed simulation model. Bourbiaux (2020) developed a model capable of optimizing the injection of low salinity water into porous media. The simulation results showed that the effect of injected fluid composition and concentration on process efficiency and residual oil saturation depends on the crude oil-brine-rock system.

The geochemical modelling of reactive flow was the subject of many studies (Amirian et al., 2019; Bhuyan et al., 1990; Safari et al., 2016; Sharma and Mohanty, 2018; Taheri et al., 2019; Younesian-Farid and Sadeghnejad, 2020). Economides and Nolte (1989) used two mineral models to describe acidic reactions during sandstone acidizing treatment. Nghiem (2003b) and Lake et al. (2002) stated the impact of proper selection of necessary reactions on the final history match results. The reaction rate of minerals, as a critical parameter, was investigated in terms of their relation with pH and temperature (Allan et al., 2011). Multiple models were developed considering the impact of injection rate, type, and concentration of injected acid (Benson, 2007; Chen et al., 2010; Choi, 2008; Kazempour and Alvarado, 2011; Kazempour et al., 2012; Mahrous et al., 2017; Patterson, 2014; Younesian-Farid and Sadeghnejad, 2019). The results show that the characterization of dominant reactions is crucial to achieving a trustable simulation geochemical model.

Much hydrocarbon remains upswept at the end of a waterflooding operation because of the injected water channeled through high-permeable zones. Several attempts have been made to increase the sweep efficiency and correct the conformance problem in heterogeneous reservoirs (de Aguiar et al., 2020, Ashrafizadeh et al., 2017; Farasat et al., 2017; Ghriga et al., 2019; Hashmet et al., 2017; Koochakzadeh et al., 2021a; Sadeghnejad and Masihi, 2017). Among them, pH-sensitive microgel injection recently received attention. This method was first introduced by Al-Anazi and Sharma (2002). As the term implies, environmental pH is an essential factor in the final efficiency of these microgels. The reservoir pH reduction by acid pre-flushing increases the performance of this method (Choi, 2008; Koochakzadeh et al., 2021b; Lalehrokh et al., 2008; Younesian-Farid and Sadeghnejad, 2019). Low pH environments decrease the polymer viscosity and enable it to penetrate deep into a reservoir. The geochemical reactions of acid and rock minerals eventually increase the system pH. Microgels swell at a critical pH and block high-permeable strata and divert the post-water flow toward low-permeable upswept zones. Some new crossed-linked pH-sensitive nanohydrogels with controllable pH changes were also developed recently (Ashrafizadeh et al., 2019; Ashrafizadeh et al., 2020). Thus, correct geochemical evaluations of acid-reservoir fluid–rock minerals are essential for a successful conformance control project (Sadeghnejad et al., 2021).

Unlike carbonate reservoirs with only one dominant reaction, modelling acid flooding in sandstone reservoirs is most complicated. Several reactions are to be considered to couple the transport with chemical reactions between the aqueous phase and sandstone minerals. The geochemical investigation of acid pre-flushing operation before pH-sensitive microgel gel flooding was extensively investigated. Choi (2005, 2008) evaluated the impact of different variables, such as the injected fluid acidity, core length, and flow rate on the final equilibrium pH. Lalehrokh et al. (2008) investigated hydrochloric acid, citric acid, and acetic acid on crushed sandstone and limestone rocks. A set of experiments to evaluate HCl injections at reservoir temperature was also analyzed (Kazempour and Alvarado, 2011). Kazempour et al. (2012) conducted multiple high-pH flooding on Berea sandstone samples and characterized mineralogy, flow rate, and composition. The dissolution of cement fractures resulting from injecting strong acids before microgel flooding was also investigated (Patterson, 2014). Experimental results achieved by Teimouri et al. (2020) showed that the rate of fracture aperture changes is a function of acid injection flow rate. Koochakzadeh et al. (2021b) conducted both batch and flooding experiments to characterize the geochemical performance of acid pre-flushing into carbonate rocks showing acetic acid as the prefered preflushing fluid to maintain a low-pH environment.

On top of discussed experimental approaches, various geochemical models were also implemented to predict the geochemical interactions during acid preflushing. The Darcy-scale models developed by Benson (2007) and Choi (2008), history matched on the experimental results of Choi (2005), failed to simulate the outlet concentration of all participated ions in reactions. Choi et al. (2006) developed a single-phase 1-D geochemical model using KGEOFLOW simulator to predict the effluent ion concentration profile during acid injection in sandstones. They evaluated the alteration of quartz, kaolinite, and anorthite as the representative of Berea sample mineralogy while considering constant rock porosity and permeability. Benson (2007) implemented GEM-GHG package to modify the developed model by Choi et al. (2006). They considered the dispersion term and porosity alteration in their study. The Berea sandstone mineralogy was simplified by considering a four-mineral model including quartz, kaolinite, potassium feldspar, and calcite. Several significant simplifications were assumed, and the impact of heterogeneity and permeability alteration on acid flooding results was disregarded.

A comparison was conducted to evaluate the performance of strong and weak acids during pH-sensitive microgel flooding (Kazempour and Alvarado, 2011; Sharma et al., 2008; Younesian-Farid and Sadeghnejad, 2019). Sharma et al. (2008) claimed that the observed mismatches of the geochemical models were due to missing some critical reactions or failure of the partial local equilibrium assumption during geochemical modelling. Kazempour and Alvarado (2011) developed a 1D radial model and 2D Cartesian model considering the proposed mineralogy by Benson (2007). They also defined the reaction rate constants of kaolinite and dolomite minerals as a function of the H⁺ concentration. Moreover, the developed geochemical model by Younesian-Farid and Sadeghnejad (2019) was calibrated on Berea sandstone rocks and then extended to carbonate rocks. Even though proper matches of the effluent pH have been achieved, the prediction of the solution chemistry changes was not accurate enough. Koochakzadeh et al. (2021b) considered permeability changes as a function of porosity and found acetic acid as a proper acid to maintain a low-pH environment prior pH-sensitive microgel injection. All discussed geochemical models were simplified through the injection of pure and additive-free acid solutions into a homogenous core sample under laboratory conditions. Disregarding the dependency of fluid properties on the fluid salinity at reservoir conditions limits the applicability of these models for prediction in high salinity reservoirs. Simplifications of the aforementioned models seem to be compensated by altering matching parameters and can cause significant errors in real case studies.

In this study, we develop a geochemical model tuned by experimental results (Choi, 2005) to predict acid injection behaviour in sandstones at reservoir condition. The rock properties alterations by the acid injection can be predicted accurately by considering influential parameters like reservoir temperature and salinity. Our model incorporates the composition and salinity of the injected fluid on the process efficiency by considering the viscosity and density of the injected fluid as a function of ionic concentrations. This can be more crucial in high salinity reservoirs such as the Middle East fields with in-situ brine salinities higher than 200 000 ppm (Al-Keebali et al., 2019; Al-Murayri et al., 2019; Jouenne, 2020; Kornilov et al., 2019; Masalmeh et al., 2019). The single-phase geochemical model in a one dimensional configuration is coded in MATLAB. The developed model is solved by a semi-implicit numerical scheme for flow transport and geochemical reaction governing equations. The heterogeneity of rock mineral distribution across the sample, absent in the previous studies, is also investigated.

2 Model development

2.1 Fluid continuity equation

The mole balance equation represents the material flowing into and out of a porous system and mass changes inside the system. The corresponding continuity equation can be expressed via the following equation for a 1D slightly compressible single-phase fluid,(1)here, α_c is the conversion factor; ϕ and ρ are the porosity and density, respectively; V_b stands for the bulk volume of the grid; A_x is the cross-sectional area; q_m is the production or injection rate, and t is time. u_x is the velocity of the bulk liquid solution in the x-direction and can be obtained using Darcy’s law. In the case of a constant injection flow rate, initial and boundary conditions can be considered as follows,(2)(3)(4)where p is the pressure; p_out stands for the pressure at the outlet of the core sample; u_inj is the injection velocity of the bulk liquid solution, and L represents the core length.

2.2 Geochemical reactions

The driving force through the bulk solution is supplied by the solution bulk flow and molecular concentration gradient. The equation of coupling multi-component chemical reactions and transport in the saturated porous media is defined as,(5)where C_i is the mass concentration of component i (mole/vol. of solution); N_t is the total number of aqueous species; R_i stands for the rate of generation or depletion of species i in the solution per unit volume, and D is the diffusion coefficient. A relation between the effective diffusion coefficient and porosity changes (e.g., similar to Archie’s law) was also introduced in previous studies (Boving and Grathwohl, 2001),(6)

Adding equation (6) into the system of equations brings another matching parameter (i.e., n). We considered a constant diffusion coefficient in this study as the convection term is dominant on diffusion in our simulations. However, such dependencies can be implemented in future studies. The considered initial and boundary conditions are then calculated by,(7)(8)(9)

Here, C_inj and C₀ are the injection and initial concentration of ions, respectively. In the case of rock–fluid interactions, the consideration of chemical equilibrium brings errors to the calculations. Therefore, a partial local equilibrium assumption was implemented to simplify the geochemical model of the present study. Aside from the instant equilibrium among aqueous phase species, kinetic reactions were considered to control the interactions between minerals and the aqueous phase. The existing components are categorized into dependent and independent variables, designated as X_j and C_i, respectively. The consideration of the equilibrium state between the two groups of components led to reducing the numbers of model unknowns as follows (Mahrous et al., 2017; Steefel et al., 2015),(10)where N_c and N_x represent the number of independent and dependent components, respectively; while ν_ij is the number of independent component i per mole of dependent component j. A_j and A_i stand for secondary and primary species, respectively. The final form of governing equation qualifying the conservation of solute mass can be written as:(11)(12)where U_i is defined as the total concentration and represents the heterogeneous reaction rate of independent component i.

It is worth mentioning that only the dominant geochemical species reactions are considered in our model. The reaction rate can be obtained using the following relation by assuming the surface reaction rate as a key controlling factor in mineral dissolution and precipitation processes,(13)where k_l is the rate constant; v_il is the stoichiometric coefficient of species i in the lth reaction; K_eql and A_l are the equilibrium constant and reactive surface area of the lth mineral, and N_rs is the number of dissolution/precipitation kinetic reactions. Q_s represents the ion activity product, which is defined using the following equation,(14)where a_i is the activity of species i. The discrepancy of ion activity and concentration is highlighted in modeling of high salinity solutions. The geochemical modeling of acid pre-flushing at reservoir conditions requires the implementation of an activity coefficient to correct this deviation. The extended Debye–Hückel model was used to calculate the activity coefficient (γ) (Bates, 1955; Kielland, 1937),(15)where I is the ionic strength; z stands for ionic charge; a corresponds to the hydrated ion size and b represents the decrease in solvent concentration. A and B are solvent-based constants and the values of 0.5125 and 0.3294 at 25 °C (0.5558 and 0.3384 at 65 °C) were considered at laboratory and reservoir conditions, respectively (Punase, 2015). The ionic strength can be computed as the following,(16)

Multiple experimental studies displayed the dependency of reaction rates constant on temperature (Abdelgawad et al., 2019; Benson, 2007). The reaction rate constant at the intended temperature can be quantified using the following equation,(17)where k_l ref is the reaction rate constant of the lth mineral at the reference temperature of T_ref; R represents the universal gas constant, and T stands for the temperature. The activity energy (E_l) of the final selected minerals was acquired from Palandri and Kharaka (2004).

The Damköhler number is defined as follow (Nghiem, 2003a),(18)where is the injected H⁺ molality and represents the maximum H⁺ reaction rate. This dimensionless controlling parameter was implemented to make a comparison between the convective flow and chemical equilibration.

2.3 Fluid properties

The dependency of the aqueous phase density on salinity was more significant than pressure (Alexeev et al., 2015; Schmidt et al., 2018). A modification was made to improve the estimation of fluid viscosity and density considering the effect of ionic concentration, which can be more crucial in high salinity sandstone reservoirs. It should be noted that these dependencies were ignored in all previous studies. In addition, the intense effects of reservoir pressure and temperature on fluid properties were considered (El-Dessouky and Ettouney, 2002). Fluid density in the continuity equation was implicitly updated via the following equations at each time step,(19)(20)(21)(22)(23)where ρ_w@S.C. is the water density at standard conditions calculated as a function of salinity (S) while B_w@Res.C. is the average water formation volume factor used to convert the standard condition into the reservoir condition.

The following equations were used to calculate fluid viscosity as a function of pressure, temperature, and salinity at each time step,(24)(25)(26)(27)where μ_w and S are the water viscosity and salinity, respectively. As illustrated in Figure 1, the temperature was the most influential parameter influencing water viscosity. Therefore, it can be inferred that fluid salinity variation cannot cause considerable alteration in water viscosity, while its effect on density is crucial and should be considered.

Fig. 1 The impact of salinity on (a) water density variation with pressure, (b) water viscosity variation with pressure, (c) water density variation with temperature, and (d) water viscosity variation with temperature.

2.4 Updating rock properties during geochemical reactions

Rock property alteration during a geochemical reaction can be predicted accurately by considering the most influential factors. It is crucial to quantify the effects of porosity and permeability alterations on transport properties. In this study, the porosity alterations caused by minerals precipitations and dissolutions were quantified by determining the variation of mineral volumes per grain Volume Fractions (VF_l). In addition, a highly correlated relation was considered between porosity and permeability, especially in sandstones (Adeoti et al., 2017). Consequently, the permeability of each grid was updated by a log-derived porosity-permeability relation to match the experimental average core permeability,(28)(29)(30)(31)where K stands for the permeability; Mw_l is the molecular weight of the lth mineral; R_l is the heterogeneous reaction rate of the lth mineral; A_l0 and VF_l0 are the initial surface area and mineral volumes per grain volume fractions of the lth mineral.

In this study, the iterative coupling approach is applied to solve the system of equations. The whole equation systems are decoupled into separate concentration and pressure equation systems, which are sequentially solved in an implicit manner by implementing the benefit of the classical Newton–Raphson approach. The advantages of this operator-splitting technique was shown by Lu and Wheeler (2009) over different numerical examples. Considering the fact that fluid continuity equation should be solved for all grids, a Jacobian matrix for N_x (grid numbers) equation was established to determine the pressure distribution along the core. In case of determining the concentration profile, the related equation was applied for all primary species (N_c) over all grids. A numerical solution of cell-centered finite difference scheme was applied for solving the system of equations.

Figure 2 provides the flow chart for our geochemical model, which implicitly calculates pressure and velocity profiles. Afterwards, fluid properties are updated as a function of pressure, temperature, and salinity before being used in the continuity equation. Next, by calculating the reaction rate of minerals at the intended temperature, an implicit method is applied to determine the concentration of particles across the core. Finally, rock properties (e.g., porosity, permeability, and mineral volume fractions) are updated. The repetition of the summarized modelling steps stops at the end of simulation time.

Fig. 2 The flowchart used for geochemical reaction modeling using the semi-implicit numerical scheme.

2.5 Model validation by experimental data

The core flooding experiments conducted by Choi (2005) are used to validate the introduced geochemical model. Choi (2005) conducted a set of experiments to evaluate the injection of 0.1 molar HCl into a Berea sandstone rock. It should be noted that the core permeability variation was reported for a core flood experiment with a constant injection rate of 2 mL/min. In this regard, the effluent pH and components concentration were applied to determine matching parameters of the present geochemical model. The overall conditions of the experiment are provided in Table 1.

Ideally, considering the entire chemical reactions that taken place in a system increases the accuracy of predictions. However, this is infeasible due to insufficient knowledge about the reactions and also considering computational expenses. Therefore, the rock mineralogy, selection of dominant reactions, and reaction rate constants are considered as matching parameters in this study. An acceptable match for effluent species concentration requires a good understanding of a complicated set of heterogeneous and homogeneous reactions. Since no exact determination of rock mineralogy was available, the typical Berea mineralogy, investigated by Benson (2007), was considered. As presented in Table 3, the mineral percentages of quartz, feldspar, kaolinite, and calcite reported by Choi (2005) were used as the initial guess.

Since the reactions between minerals and species present in the aqueous phase take place much slower, Sevougian et al. (1995) suggested the kinetic-type reactions for dissolution/precipitation and the local equilibrium assumption for inter-aqueous reactions. In addition, the equilibrium equations presented by Morgan and Stumm (1996) were implemented to determine the dominant intra-aqueous reactions. The final selected heterogeneous and homogeneous reactions and their equilibrium constants are given in Tables 2 and 3, respectively. The initial estimates for reaction constants were obtained from Palandri and Kharaka (2004).

2.6 High-salinity sandstones

The validated model from the previous section was used as a based model to evaluate the impact of multiple parameters on core flooding results. First, the base model was saturated with formation water and was flooded with a 0.1 M HCl solution. In addition, the pressure and temperature of the model were set to 170 atm and 150 °F, respectively. Finally, a set of cases was simulated to study the effect of fluid injection rate, composition, and the dependency of fluid properties on fluid salinity. In all cases, the dependency of fluid properties and the reaction rate constants on temperature were also considered.

In order to evaluate acid pre-flushing in high-salinity reservoir conditions, the composition of formation water sampled from a formation (Y-Field) in the south of Iran was implemented. This downhole water sample was taken during a Modular Dynamic Test (MDT) operation and had a pH of 6.3 and the Total Dissolved Solids (TDS) of 293 000 mg/L and considered as a formation brine in our cases. The composition and corresponding stiff diagram are presented in Table 4 and Figure 3, respectively.

Fig. 3 Stiff diagram of formation water bottom-hole sample taken during MDT operation from Y-Field implemented as the pore space fluid.

Three HCl solutions with the same injection pH of 1 but prepared with brine solutions with various salinities were considered to analyze the effect of injected fluid composition on modelling results. Table 5 presents these three fluids. The composition of low salinity and the seawater sample was selected from Rahimi et al. (2020).

3 Results and discussion

The developed model capable of simulating a wide variety of homogenous and heterogeneous reactions is verified using experimental data reported by Choi (2005). This section is divided into five sections starting with grid sensitivity analysis and continuing with model validity analysis. The effect of multiple parameters, including injection fluid composition, injection rate, and heterogeneous minerals distribution along the core was investigated. In this regard, a set of core flooding cases were modelled and thoroughly discussed in this section.

3.1 Grid independence analysis

A sensitivity analysis was performed on the number of grids and time steps to determine the optimum values considering model run-time and simulation accuracy. The experimental results were compared with models composed of 20 through 70 grid blocks to find the model sensitivity to grids. The result of this comparison for pH (Fig. 4a) and effluent K composition (Fig. 4b) is presented as an example. As indicated, an insignificant difference was found for the models with more than 40 grid blocks. So, this value was selected for our further simulations.

Fig. 4 Sensitivity analysis on the impact of grid number on (a) pH and (b) K concentration at the outlet of the Berea sandstone core.

As mentioned earlier, the fluid continuity equation and geochemical transport reactions are separately solved in an implicit manner. Since the modeling steps were treated sequentially, the solving procedure is not unconditionally stable and requires some restrictive precautions on the time step to avoid convergence failure. Taking into the account that most dramatic changes in unknowns (i.e., pressure and concentrations) occur during the initial simulation stages, small time steps were considered. As the simulations progress, time steps were multiplied by a constant coefficient (time step multiplier coefficient, δ) to speed up the simulations. To find the proper value for δ, a sensitivity analysis was performed. A simulation with a fixed small time step (t_f = 0.01 s) was run, which has the benefit of the excellent solution accuracy. The experimental results were compared to the models with various time step multiplier coefficients to find the model sensitivity to time steps. Finally, the time step multiplier coefficient at later time steps was determined based on the deviation from the constant time step method. The result of this comparison for pH (Fig. 5a) and effluent K composition (Fig. 5b) is presented as an example. As it is indicated, the models with δ of 1.2 indicated similar accuracy to the constant time step method. Thus, this value was selected for the rest of our simulations.

Fig. 5 Sensitivity analysis on the impact of time step size on (a) pH and (b) K concentration at the outlet of the Berea sandstone core.

3.2 Model validation

The final matching results of the mineral reaction rate constant and the grain composition as a fraction of bulk volume and initial specific reaction area are provided in Table 6. The final match was more sensitive to the calcite reaction rate constant than the other reactions. Moreover, the Kaolinite and potassium feldspar reaction constants were considered a function of hydrogen concentration in the system.

A comparison between the results of our model and previous studies along with the experimental observations is provided in Figure 6. Because a tracer dose does not participate in reactions, its injection was monitored to evaluate the convective mass transport across the core. According to the simulation results, it only takes 400 s for the tracer to cross the core length. Therefore, it can be concluded that the first significant drop in pH trend can be assumed due to H⁺ crossing the core length by flow (Fig. 6a). It should be noted that H⁺ ions participate in multiple reactions, and their transport is kinetically controlled (see Tables 2 and 3). The second sharp decline observed in the pH trend (Fig. 6a) was primarily because of the consumption of the main content of calcite minerals in the rock sample. As a result, the calcite reaction rate constant was tuned while matching the second observed pH drop after about 146 injection pore volume (elapsed-time of 1300 min). This phenomenon is associated with a sharp decrement in the effluent Ca concentration.

Fig. 6 Comparison between the model predictions and experimental results of Choi (2005) for (a) the effluent pH; (b) average core permeability; (c) aluminium; (d) calcium; (e) potassium, and (f) silicon concentrations at the outlet of Berea sandstone core.

The increasing H⁺ concentration in the system is expected to result in a higher dissolution rate of minerals at the nearby inlet regions by the start of HCl injection. The dependency of kaolinite mineral reaction on pH, as the only source for Al in effluents, was assumed as the reason for the upward trend observed in the effluent Aluminium concentration (Fig. 6c). The lower dependency of potassium (Fig. 6e) and silicon (Fig. 6f) on pH resulted in negligible variations observed for their concentration profiles. The reasonable assumptions and characterization of the reactions provided a higher consistency between model results and the experimental observations than previous studies (Benson, 2007; Choi, 2008; Kazempour and Alvarado, 2011; Sharma et al., 2008; Younesian-Farid and Sadeghnejad, 2019). As shown in Figure 6b, the average permeability at the end of the experiment exceeded almost twice its initial value. A higher permeability increment was observed after 1300 min of injection when the entire core was almost under minimum pH. It is worth mentioning that permeability alteration was not reported in the previous studies.

3.3 Reaction rate versus convective mass transport

The Damköhler number (Da) was implemented to evaluate the core flooding model behaviour at reservoir conditions. When this dimensionless number has a low value, the convective flow term is more dominant than the reaction term in the system, which means that the reaction process is controlled by the transport of active ions. In addition to the modelling of the experiment by Choi (2005) with a Da of 5, two sets of simulations with a Da of 1 and 10 were considered by changing the injection rate (2.12, 0.21 cm/min).

As illustrated in Figure 7a, a constant pH of 7 can be observed for the model with a zero-injection rate. This indicates an equilibrium condition across the core. This model was used as an initial condition for the rest of the analyses. It should be noted that the final pH value depends on the model properties such as mineralogy and dominant reactions. The final pH of 1, equal to the injection pH, was achieved for all three cases (Da = 1, 5, 10). However, the time reaching this pH value increases by lowering the injection rate and increasing residence time. The time elapsed to achieve the final pH value for the injection case with the Da of 10 and 5 is about 1000 and 500 min, respectively. Nevertheless, it takes only a few seconds for the highest injection rate case (Da = 1) to reach this pH value. As expected, higher injection rates (e.g., 2.12 cc/min in Da = 1) highlighted the effect of convective mass transport over the reaction term. In other words, increasing the residence time for active species in the system allows heterogeneous reactions to approaching their equilibrium condition before leaving the system.

Fig. 7 (a) pH, (b) K-Feldspar composition, and (c) porosity changes at the outlet of the Berea sandstone core for injections with various Das.

As mentioned earlier, the reaction rate constant of some minerals (Kaolinite and K-Feldspar) was found dependent on the H⁺ concentration (see Table 6). Achieving the final pH value with a lower elapsed time for the injection with the Da of 1 resulted in a higher reaction rate constant for those minerals. For example, according to Figure 7b, approximately 70% of the K-Feldspar content was reduced at the end of injection for the case with the Da of 1, while this value was 50% for the case with the Da of 10. It seems that the impact of the H⁺ concentration dominated the residence time effect. The average porosity changes over the injection period are illustrated in Figure 7c. The higher H⁺ concentration in the case with a higher reaction rate (Da = 1) resulted in a sharp increment in porosity, especially at earlier times (T < 400 min). Afterwards, the slope of the porosity curve was reduced by getting closer to the equilibrium condition. Meanwhile, reaching an equilibrium takes more time in cases with high Da. For example, in the case with the Da of 10, porosity increases with an approximately constant slope during the injection phase. Thus, an optimum injection rate can be found considering existing practical limitations for injection rate and time elapsed to achieve the desired pH value.

3.4 Evaluation of fluid composition

As illustrated in Figure 8a, three injection cases with different water compositions (see Table 5) but identical pH was modelled and compared with the case of 0.1 M HCl acid injection with zero salinity. Three water samples were used for saline water injection cases as the injection fluids with total dissolved salinities of 35 000, 28 000, and 21 000 ppm. Since underlying geochemical reactions depend on the type and the number of existing components, the injection of saline waters creates a new chemical state in the porous medium. Thus, the rate of both homogeneous and heterogeneous reactions and even their direction might change. The results of this evaluation prove that in addition to the impact of injection pH investigated in the previous studies, the process efficiency can be altered by the composition of injected fluid (Al-Otaibi, 2020; Ghasemian et al., 2019; Hu et al., 2020; Olayiwola and Dejam, 2020). Comparing the injecting fluids with different salinity in Figure 8a indicates that the available active ions in the system created a new chemical state with a lower tendency toward minerals dissolution. As salinity increases, more ions are coming into the porous system with the injected fluid. For example, in the case of seawater injection, 419 ppm Ca was injected with a constant rate of 0.4 cc/min and reduced calcite dissolution rate. Therefore, more remarkable porosity changes are expected for the case with an injection rate of 0.1 M HCl solution with zero salinity. As illustrated in Figure 8a, a lower final pH value (0.8) at the end of the injection period was observed for the case with zero salinity. Figure 8b illustrates different final average porosities for the four saline cases because of each system’s various geochemical conditions. Considering a direct relation between the mineral reaction rate and the H⁺ concentration, the low pH observed for the zero-salinity case resulted in high minerals reaction rates and core porosity changes. It can be inferred that the composition of the injected fluid can alter the system equilibrium, which accounts for an essential parameter in modelling reactive flow in porous media.

Fig. 8 Comparison of (a) effluent pH and (b) average porosity for cases with different injection salinities. (c) Comparison of the effluent pH for cases with and without considering the variation of fluid properties with salinity.

A higher density and viscosity are expected for injection cases with higher salinity by considering the dependency of fluid properties on salinity (especially density). Therefore, it can be assumed that apart from the new chemical state, the observed discrepancies in the outlet pH can be due to different fluid properties of the developed cases. A comparison was made for the seawater injection model by considering and not considering the effect of salinity on fluid properties (density and viscosity changes) (Fig. 8c). As indicated, ignoring the variation of fluid property by salinity has a more dominant effect on model predictions at early simulation times. An effluent pH of 3.2 was observed for the case with the constant fluid property after 60 min from the start of the seawater injection, while the corresponding value for the case with the variable fluid property case was around 1.4. This deviation can be more crucial in cases with higher salinities. As expected, the final equilibrium pH of both cases was identical since it only depends on the geochemical specifications of the model (e.g., mineralogy and present reactions).

3.5 Heterogeneous mineral distribution

The impact of non-uniformity of minerals distribution across the core was evaluated. To do so, the four rock minerals (see Table 6) were distributed across the core by considering a standard deviation. The mean content for each mineral across the sample was maintained fixed to keep the geological properties constant. It should be noted that the model was initialized with a constant porosity of 0.177. Figure 9 illustrates the final volumetric percentage of four minerals of quartz, feldspar, calcite, and kaolinite, randomly distributed across the core. Since chemical reactions occur concerning the content of minerals, different chemical states were expected for each grid block. Figure 10a illustrates the porosity alteration across the core at different elapsed times from the start of injection. Since the mineral content in the intended grid directly affects mineral reaction rate, heterogeneous mineral distribution led to oscillations in the mineral dissolution and precipitation across the core and resulted in variable porosity changes. As shown in Figures 10b and 10c, the effluent pH and Al concentration were mainly identical for both heterogeneous and homogeneous models. Therefore, it can be inferred that the outlet results are only affected by the average mineral contents within the rock system, which means that the equivalent condition of a system depends on its average properties and not its local variations.

Fig. 9 The volumetric percentages of four main minerals present in Berea sandstone along the core in the heterogeneous model.

Fig. 10 (a) Porosity changes across the core for the heterogeneous model at different elapsed times after the start of injection; Comparison of (b) the effluent pH and (c) Al concentration of heterogeneous and homogeneous models.

4 Conclusion

Incompatibility of the injected fluid with the pore fluid during acid pre-flushing can result in a chemical non-equilibrium state within the porous medium. Therefore, different chemical reactions can occur depending on the type and the number of existing components in the system. In this study, a generalized model capable of simulating a wide variety of homogenous and heterogeneous reactions was developed to predict the pH and ions concentration while considering alterations of rock and fluid properties. The material balance, Darcy law, and equilibrium and rate-dependent reactions balance were coupled to quantify the governing mechanisms during the reactive geochemical transport. The model has been coded in MATLAB and validated using the experimental data presented by Choi (2005). It can be inferred that an accurate model was developed by the variation of rock and fluid properties into account. The initialized model under reservoir conditions was used to evaluate the core flooding performance by altering multiple items, including heterogeneity, injected fluid composition, Damköhler number (Da), and dependency of fluid physical properties on salinity.

Despite the observed discrepancy in pH and average porosity of the investigated models with various Da numbers, the final pH values were found identical for all the cases because of their similar geochemical characterizations. In addition, modelling the injection with various salinities provided a secondary source for ion supplement into the core system. This supply of active ions into the system provided a new chemical state with a lower tendency toward mineral dissolution. This study proves that in addition to injection pH, the composition of the injected fluid can alter the efficiency of the process. Even though the final pH value does not depend on the variation of fluid properties with salinity. However, considering constant fluid properties was associated with numerous errors, especially for the early-time period. In other words, fluid properties dependency on salinity and the composition of injected solution are essential enough to be considered for mathematical simulation of such a geochemical fluid flow process. Finally, modelling the non-uniform distribution of minerals across the core indicated that the outlet results are only affected by the average mineral contents in the system.

References

All Tables

All Figures

	Fig. 1 The impact of salinity on (a) water density variation with pressure, (b) water viscosity variation with pressure, (c) water density variation with temperature, and (d) water viscosity variation with temperature.
In the text

	Fig. 2 The flowchart used for geochemical reaction modeling using the semi-implicit numerical scheme.
In the text

	Fig. 3 Stiff diagram of formation water bottom-hole sample taken during MDT operation from Y-Field implemented as the pore space fluid.
In the text

	Fig. 4 Sensitivity analysis on the impact of grid number on (a) pH and (b) K concentration at the outlet of the Berea sandstone core.
In the text

	Fig. 5 Sensitivity analysis on the impact of time step size on (a) pH and (b) K concentration at the outlet of the Berea sandstone core.
In the text

	Fig. 6 Comparison between the model predictions and experimental results of Choi (2005) for (a) the effluent pH; (b) average core permeability; (c) aluminium; (d) calcium; (e) potassium, and (f) silicon concentrations at the outlet of Berea sandstone core.
In the text

	Fig. 7 (a) pH, (b) K-Feldspar composition, and (c) porosity changes at the outlet of the Berea sandstone core for injections with various Das.
In the text

	Fig. 8 Comparison of (a) effluent pH and (b) average porosity for cases with different injection salinities. (c) Comparison of the effluent pH for cases with and without considering the variation of fluid properties with salinity.
In the text

	Fig. 9 The volumetric percentages of four main minerals present in Berea sandstone along the core in the heterogeneous model.
In the text

	Fig. 10 (a) Porosity changes across the core for the heterogeneous model at different elapsed times after the start of injection; Comparison of (b) the effluent pH and (c) Al concentration of heterogeneous and homogeneous models.
In the text