Jan Vinogradov

Miftah Hidayat

Jos Derksen

David Vega-Maza

Stefan Iglauer

Mohammad Sarmadivaleh

Damien Jougnot

Mohamed Azaroual

Philippe Leroy

Predictive surface complexation model of the calcite-aqueous solution interface: The impact of high concentration and complex composition of brines

Introduction

Calcite is a very common mineral on the Earth's surface comprising approximately 4% by weight of Earth's crust. Physico-chemical processes that take place in porous calcite-aqueous solution (water for simplicity) systems are of great importance for a variety of scientific and technological areas. Properties of the calcite mineral surface and the interface between calcite and aqueous solution are important for a broad range of applications including geological storage of nuclear waste and CO2 [START_REF] Hester | Carbon Capture: Sequestration and Storage[END_REF], freshwater aquifer management [START_REF] Graham | Self-Potential as a Predictor of Seawater Intrusion in Coastal Groundwater Boreholes[END_REF][START_REF] Macallister | Remote Detection of Saline Intrusion in a Coastal Aquifer Using Borehole Measurements of Self-Potential[END_REF], hydrocarbon production [START_REF] Saunders | Streaming potentials at hydrocarbon reservoir conditions[END_REF], paper production [START_REF] Carrasco | Control of retention in paper-making by colloid titration and zeta potential techniques[END_REF]. However, electrochemical properties of the interfaces between calcite and aqueous solutions (specifically, electrical surface charge and zeta potential), especially under conditions relevant to natural subsurface systems (i.e., water chemical composition, temperature, pressure) remain poorly understood because of their microscopic nature, the high reactivity of the calcite mineral, and the lack of relevant measurements and models [START_REF] Heberling | Reactivity of the calcite-water-interface, from molecular scale processes to geochemical engineering[END_REF].

The zeta potential is defined as the electrical potential at the shear plane when water flow occurs parallel to the pore surface [START_REF] Hunter | Zeta potential in colloid science: principles and applications[END_REF]. It is most commonly used to estimate the electrochemical properties of calcite because acid-base potentiometric titration, usually employed to measure the surface charge of minerals, is extremely difficult to carry out for highly reactive calcite [START_REF] Eriksson | The calcite/water interface: I. Surface charge in indifferent electrolyte media and the influence of low-molecular-weight polyelectrolyte[END_REF]. There are many studies that have reported measurement of the zeta potential of various calcite-water systems. However, many reported results are contradictory (see for example Fig. 3 in Jackson et al. [START_REF] Jackson | Measurements of spontaneous potential in chalk with application to aquifer characterization in the southern UK, Q[END_REF]). Most of the experimental studies have reported zeta potential measured with low salinity single-salt solutions (e.g., [START_REF] Heberling | Structure and reactivity of the calcite-water interface[END_REF][START_REF] Li | Influence of surface conductivity on the apparent zeta potential of calcite[END_REF]), and there is inconsistency in the experimental conditions when some experiments were conducted at equilibrium between the mineral and water while others reported the zeta potential with non-equilibrium solutions (see for example discussion in Section 3 and 4 Al Mahrouqi et al. [START_REF] Al Mahrouqi | Zeta potential of artificial and natural calcite in aqueous solution[END_REF]). The experimental studies have also employed a variety of electrokinetic methods including the electrophoretic mobility (EPM [START_REF] Foxall | Charge determination at calcium salt/aqueous solution interface[END_REF]) and streaming potential method (SPM; [START_REF] Thompson | Surface electrical properties of calcite[END_REF]). Each of these experimental techniques suffers from several limitations [15]. The EPM applies an electrical field (of magnitude hundred Volts per meter) that mobilizes charged particles relative to aqueous solution and measuring the resulting mobility. The EPM is conducted on powdered calcite samples suspended in aqueous solution thus potentially exposing fresh mineral surfaces. Moreover, the method is usually limited to solution concentration below 1 M (M = mol/L) [START_REF] Kosmulski | High ionic strength electrokinetics of clay minerals, Colloids Surfaces A Physicochem[END_REF]. The ratio of the mineral to the solution is also very low (typically few grams per liter) and not representative of real rocks. On the other hand, the SPM applies a pressure difference across the sample and measures a resulting voltage, it can be carried out under equilibrium conditions, on intact rock samples at salinities up to full saturation (e.g., [START_REF] Cherubini | Influence of CO2 on the Electrical Conductivity and Streaming Potential of Carbonate Rocks[END_REF][START_REF] Cherubini | Streaming Potential Coupling Coefficient and Transport Properties of Unsaturated Carbonate Rocks[END_REF]). However, unlike the EPM method, SPM can be time consuming and very challenging for highly reactive minerals such as calcite in highly saline conditions thus limiting its ability to cover a wide range of rock-solution permutations (see for example typical duration of both the equilibration process and individual SPM experiment, as well as the voltage stability issues at high salinity reported in Alroudhan et al. [START_REF] Alroudhan | Zeta potential of intact natural limestone: Impact of potential-determining ions Ca, Mg and SO4, Colloids Surfaces A Physicochem[END_REF]). Moreover, the zeta potential values obtained with SPM can be specific to rock type as different carbonate samples may contain different impurities such as anhydrite and dolomite (e.g., [START_REF] Al Mahrouqi | Zeta potential of artificial and natural calcite in aqueous solution[END_REF]).

A model capable of predicting the zeta potential of calcite-water systems at equilibrium as well as non-equilibrium conditions, low to high salinity, and complex solution composition is crucially important. Such model will not only be capable of predicting the zeta potential, but also improve our understanding of the complex electrochemical processes that take place at calcite-water interface. To date, there have been numerous attempts to develop such model with most of published studies using the so-called surface complexation modelling (SCM) approach (e.g., [START_REF] Heberling | Reactivity of the calcite-water-interface, from molecular scale processes to geochemical engineering[END_REF][START_REF] Heberling | Structure and reactivity of the calcite-water interface[END_REF][START_REF] Li | Influence of surface conductivity on the apparent zeta potential of calcite[END_REF]). In this method the calcite surface is presented as an ensemble of chemically active surface sites that interact with ions from the adjacent bulk electrolyte. The strength of these interactions is described by the equilibrium constants (similar to chemical reactions), so that the resulting equilibrium concentration of positive and negative surface complexes establishes the net surface charge and the corresponding surface potential [START_REF] Davis | Surface ionization and complexation at the oxide/water interface: I. Computation of electrical double layer properties in simple electrolytes[END_REF]. Depending on the complexity of the aqueous solution, such models can be realized analytically (e.g., [START_REF] Revil | Streaming potential in porous media: 1. Theory of the zeta potential[END_REF]) or by using numerical methods already available in software packages such as PHREEQC (e.g., [START_REF] Parkhurst | Description of Input and Examples for PHREEQC Version 3-A Computer Program for Speciation, Batch-Reaction, One-Dimensional Transport, and Inverse Geochemical Calculations[END_REF]).

Studies published to date on SCM have used different approaches to model the zeta potential (Figure 1): Basic Stern Model (BSM, e.g., [START_REF] Heberling | Structure and reactivity of the calcite-water interface[END_REF]), Triple Layer Model (TLM, e.g., [START_REF] Wolthers | The surface chemistry of divalent metal carbonate minerals; a critical assessment of surface charge and potential data using the charge distribution multi-site ion complexation model[END_REF]), Quad-Layer Model (QLM, e.g. [START_REF] Alizadeh | Flexibility of inactive electrokinetic layer at charged solid-liquid interface in response to bulk ion concentration[END_REF]) and Diffuse Layer Model (DLM, e.g. [START_REF] Ding | Investigation of the Impact of Potential Determining Ions from Surface Complexation Modeling[END_REF]). BSM, TLM and QLM explicitly describe the Stern layer and a layered arrangement of ionic species adsorbed on the mineral surface in the diffuse layer. Both the TLM and QLM distinguish between the Inner Helmholtz plane (IHP) and the Outer Helmholtz plane (OHP), at which adsorption of de-hydrated and hydrated ions, respectively, takes place [START_REF] Wolthers | The surface chemistry of divalent metal carbonate minerals; a critical assessment of surface charge and potential data using the charge distribution multi-site ion complexation model[END_REF]. On the other hand, BSM for calcite assumes that adsorption of all ions in compact layer takes place at the OHP (e.g., [START_REF] Heberling | Reactivity of the calcite-water-interface, from molecular scale processes to geochemical engineering[END_REF][START_REF] Heberling | Structure and reactivity of the calcite-water interface[END_REF][START_REF] Li | Influence of surface conductivity on the apparent zeta potential of calcite[END_REF]). All three models that explicitly describe the adsorption layers (i.e., BSM, TLM, QLM) also assume that only protonation/deprotonation of hydrated calcite lattice ions take place at the mineral surface. In these models, the surface charge densities are computed at each plane, from which linear, capacitor-like variation of the electric potential between the planes is obtained. Consequently, TLM and QLM simulate separately protonation/deprotonation and salt ion adsorption reactions at three different planes (mineral surface, IHP, OHP), and ascribe two capacitance values to layers confined between the mineral surface and IHP, and between IHP and OHP. In addition, the QLM considers a stagnant diffuse layer, implying that the shear plane is located further away from the OHP. This makes these two models more computationally expensive compared with BSM, which uses only one capacitance and all adsorption reactions are considered to take place at OHP for calcite.

On the other hand, DLM assumes that both the protonation/deprotonation and salt ion adsorption reactions take place at the mineral surface thus ignoring the complexity of the Stern layer that comprises the IHP and OHP. Hence, compared to the BSM, TLM, and QLM, the DLM is not sensitive to the capacitance values and is only sensitive to the location of the so-called shear (or slip) plane at which ions can be mobilized by the flow and at which the electric potential is defined as the zeta potential. According to the generally agreed theory of the electrical double layer, concentrations of counter-ions that populate the region beyond the OHP (in the diffuse layer) obey the Boltzmann distribution law where they decrease exponentially towards concentrations of ions in electroneutral bulk electrolyte (outside the diffuse layer). The location of the slip plane is uncertain and is usually used as a fitting parameter to match the model predictions to experimental data [START_REF] Heberling | Structure and reactivity of the calcite-water interface[END_REF]. The assumed location of the slip plane can also be used in BSM, TLM and QLM as additional fitting parameter, but in contrast with DLM the modelled zeta potential using these models depends on adjustable capacitance values and therefore, the location of the slip plane can be kept constant for matching different experimental datasets. The QLM published by Alizadeh and Wang [START_REF] Alizadeh | Flexibility of inactive electrokinetic layer at charged solid-liquid interface in response to bulk ion concentration[END_REF] assumed that the slip plane was dynamic and moved further away from the mineral surface as the ionic strength of the bulk electrolyte decreased. Although the dynamic slip plane model successfully reproduced experimental results, the hypothesis was not clearly justified by the authors or confirmed experimentally. at each plane, the slip plane does not coincide with OHP and its distance from OHP can be constant or varying with salinity); d) -DLM (IHP and OHP are not considered at all as if there is no separation between the mineral surface and OHP, the slip plane may or may not coincide with mineral/OHP).

However, despite the plethora of published SCM studies, all of the simulated results appear to be specific to an experimental dataset. Heberling and co-authors [START_REF] Heberling | Reactivity of the calcite-water-interface, from molecular scale processes to geochemical engineering[END_REF][START_REF] Heberling | Structure and reactivity of the calcite-water interface[END_REF] used fractional charges of individual surface site, the approach that was consistent with that taken by Wolthers et al. [START_REF] Wolthers | The surface chemistry of divalent metal carbonate minerals; a critical assessment of surface charge and potential data using the charge distribution multi-site ion complexation model[END_REF], who used the Charge Distribution MUltiSite Ion Complexation or CD-MUSIC modeling approach. However, the magnitude of the surface site charges was different. On the other hand, Van Cappellen et al. [START_REF] Van Cappellen | A surface complexation model of the carbonate mineral-aqueous solution interface[END_REF] used integer surface site charges in their BSM. Furthermore, Heberling et al. [START_REF] Heberling | Reactivity of the calcite-water-interface, from molecular scale processes to geochemical engineering[END_REF][START_REF] Heberling | Structure and reactivity of the calcite-water interface[END_REF] used equilibrium constants for adsorption reactions at OHP limited by ion-binding ones that corresponded to pairing of ions in the bulk electrolyte. In contrast, Pokrovsky and Schott [START_REF] Pokrovsky | Surface Chemistry and Dissolution Kinetics of Divalent Metal Carbonates[END_REF] estimated their equilibrium constants for adsorption reactions from the correlation between stability of surface sites and stability of the same molecules in the solution. As a result, the equilibrium constants used by Pokrovsky and Schott [START_REF] Pokrovsky | Surface Chemistry and Dissolution Kinetics of Divalent Metal Carbonates[END_REF] in their BSM exceeded by several orders of magnitude those used by Heberling et al. [START_REF] Heberling | Reactivity of the calcite-water-interface, from molecular scale processes to geochemical engineering[END_REF][START_REF] Heberling | Structure and reactivity of the calcite-water interface[END_REF] in their BSM. Considering the equilibrium constants of protonation/deprotonation reactions that take place at the mineral surface, a significant variation in their values was also evident across different studies. For instance, Heberling et al. [START_REF] Heberling | Structure and reactivity of the calcite-water interface[END_REF] used fivefold larger equilibrium constant for deprotonation of hydrated Ca site compared with their later study [START_REF] Heberling | Reactivity of the calcite-water-interface, from molecular scale processes to geochemical engineering[END_REF]. At the same time, Wolthers et al. [START_REF] Wolthers | The surface chemistry of divalent metal carbonate minerals; a critical assessment of surface charge and potential data using the charge distribution multi-site ion complexation model[END_REF] used twelve orders of magnitude larger equilibrium constant for the same reaction compared with Heberling et al. [START_REF] Heberling | Structure and reactivity of the calcite-water interface[END_REF].

Previously published studies that used either BSM or TLM relied on constant Stern layer capacitance, but the reported values used in these studies varied between 0.45 F/m 2 , 1.24 F/m 2 [START_REF] Li | Influence of surface conductivity on the apparent zeta potential of calcite[END_REF] , 17 F/m 2 [START_REF] Pokrovsky | Surface Chemistry and Dissolution Kinetics of Divalent Metal Carbonates[END_REF], 52 F/m 2 [START_REF] Pokrovsky | Surface Speciation Models of Calcite and Dolomite/Aqueous Solution Interfaces and Their Spectroscopic Evaluation[END_REF] and 100 F/m 2 [START_REF] Wolthers | The surface chemistry of divalent metal carbonate minerals; a critical assessment of surface charge and potential data using the charge distribution multi-site ion complexation model[END_REF]. Therefore, the inconsistency in these values suggests a high degree of freedom and unconstrained variation of simulated zeta potential. Moreover, most of the published SCM simulated the zeta potential of calcite in contact with low salinity (typically less than 0.1 M) single salt solutions (e.g., NaCl) therefore ignoring the effect of divalent ions such as Mg 2+ and SO4 2-and high ionic strength. Studies that included such reactions reported inconsistent equilibrium constants, with Heberling et al. [START_REF] Heberling | Reactivity of the calcite-water-interface, from molecular scale processes to geochemical engineering[END_REF] predicting adsorption of Ca 2+ on negative surface sites to be controlled by an equilibrium constant ten orders of magnitude smaller compared with the equilibrium constant calculated from two consecutive reactions described by Pokrovsky et al. [START_REF] Pokrovsky | Surface Speciation Models of Calcite and Dolomite/Aqueous Solution Interfaces and Their Spectroscopic Evaluation[END_REF] and Pokrovsky and Schott [START_REF] Pokrovsky | Kinetics and Mechanism of Dolomite Dissolution in Neutral to Alkaline Solutions Revisited[END_REF]. All in all, a unique combination of adsorption reactions, their equilibrium constants, capacitance values and charges of the surface sites used in published studies appear to be very specific to the modelled experimental datasets, thus limiting the use of each model for a very narrow range of ionic strength, pH and composition of electrolytes in contact with the calcite mineral.

Therefore, the aim of this study is to develop and validate a new SCM that accurately computes the zeta potential of any calcite-water system as long as the experimental conditions are reproduced in the model. The developed model has proven to be accurate in simulating and predicting zeta potentials at both equilibrium and non-equilibrium conditions, high to low salinity, for any carbonate rock type as long as the dominating mineral is calcite, and for any solution composition as long as the dominating salt is NaCl and concentration of SO4 2-does not exceed that of Ca 2+ and Mg 2+ .

Methodology

To develop a robust predictive model requires detailed description of experimental parameters and conditions. Therefore, we have chosen the data reported by Al Mahrouqi et al. [START_REF] Al Mahrouqi | Zeta potential of artificial and natural calcite in aqueous solution[END_REF] as the most comprehensive experimental work and will briefly describe the experimental conditions and main conclusions.

2.1.

Description of the experimental data on zeta potential in intact limestone samples

Streaming potential measurements were performed by Al Mahrouqi et al. [START_REF] Al Mahrouqi | Zeta potential of artificial and natural calcite in aqueous solution[END_REF] and Alroudhan et al. [START_REF] Alroudhan | Zeta potential of intact natural limestone: Impact of potential-determining ions Ca, Mg and SO4, Colloids Surfaces A Physicochem[END_REF] on three different carbonate rock samples: Estaillades [START_REF] Al Mahrouqi | Zeta potential of artificial and natural calcite in aqueous solution[END_REF], Ketton [START_REF] Al Mahrouqi | Zeta potential of artificial and natural calcite in aqueous solution[END_REF] and Portland [START_REF] Alroudhan | Zeta potential of intact natural limestone: Impact of potential-determining ions Ca, Mg and SO4, Colloids Surfaces A Physicochem[END_REF]. The petrophysical and mineralogical properties of the rock samples are shown in Table 1. amineralogy reported by Al Mahrouqi et al. [START_REF] Al Mahrouqi | Zeta potential of artificial and natural calcite in aqueous solution[END_REF] using X-ray diffraction (XRD). Magnesium was reported to be likely incorporated into the limestone as dolomite.

Property
brecent study by Udoh and Vinogradov [START_REF] Udoh | Controlled Salinity-Biosurfactant Enhanced Oil Recovery at Ambient and Reservoir Temperatures-An Experimental Study[END_REF], identified and confirmed that there was 4% dolomite (CaMg(CO3)2) and 1% anhydrite (CaSO4) in the Estaillades sample using XRD.

The zeta potential measurements were conducted at room temperature (23 o C) with NaCl solutions of varying concentration between 0.05 M (mol/L) and 5.0 M NaCl at equilibrium conditions. Effluent fluid samples were regularly collected for chemical analyses to determine the concentration of Na + , Ca 2+ , Mg 2+ , Cl -, and SO4 2-. The experimentally confirmed presence of Ca 2+ , Mg 2+ , and SO4 2-was explained by partial dissolution of rock samples during the equilibration with NaCl solutions and leaching of these ions.

These studies concluded that the zeta potential of carbonate rock samples was controlled by the concentration of the potential determining ions (PDIs), which were identified to be Ca 2+ , Mg 2+ , and SO4 2-.

Calcite surface complexation model development

Basic definitions and assumptions

A BSM was used in this study to model the zeta potential of pure calcite in contact with aqueous solutions. The modelling approach was consistent with the most comprehensive and empirically justified model developed by Heberling et al. [START_REF] Heberling | Reactivity of the calcite-water-interface, from molecular scale processes to geochemical engineering[END_REF][START_REF] Heberling | Structure and reactivity of the calcite-water interface[END_REF]. Since the content of calcite in all three tested carbonate rock samples (Estaillades, Ketton, Portland) was found to be greater than 95%, it was assumed that calcite was the most dominant mineral in the samples so that only calcite surface sites were considered in the model.

The measured zeta potential with all tested carbonate rocks reflected presence of minerals other than calcite (e.g., dolomite and anhydrite), which led to a non-zero concentration of Mg 2+ and SO4 2- found in the bulk solution as a result of mineral dissolution during the equilibration process. Therefore, we artificially added the equilibrium concentration of Mg 2+ and SO4 2-measured in the experiment to our geochemical model of pure calcite. The amount of Ca 2+ leached into NaCl solutions during the equilibration was either added artificially to the model and equilibration was switched off, or it was not added to the initial solution while the equilibration was activated (see Section 2.3 for details).

Consistent with Heberling et al. [START_REF] Heberling | Reactivity of the calcite-water-interface, from molecular scale processes to geochemical engineering[END_REF][START_REF] Heberling | Structure and reactivity of the calcite-water interface[END_REF], calcite (1 0 4) face was assumed to dominate the interface due to its high abundance and the total surface site density of calcium and carbonate surface sites on the calcite (1 0 4) face was set to 4.95 sites per nm 2 for each type of site. However, unlike Heberling et al. [START_REF] Heberling | Reactivity of the calcite-water-interface, from molecular scale processes to geochemical engineering[END_REF], this study assumed that the zeta potential could be calculated by considering the shear plane coinciding with the OHP (see Figure 2), i.e., 𝜙 𝛽 = 𝜙 𝑑 = 𝜁the assumption that was also used by Li et al. [START_REF] Li | Influence of surface conductivity on the apparent zeta potential of calcite[END_REF] (no stagnant diffuse layer). Li et al. [START_REF] Li | Influence of surface conductivity on the apparent zeta potential of calcite[END_REF] considered that the assumption of the presence of a stagnant diffuse layer in Heberling et al. [START_REF] Heberling | Structure and reactivity of the calcite-water interface[END_REF] was due to the nature of electrokinetic measurements that were not corrected for surface conductivity effects thus reducing the experimental zeta potentials and then necessitating smaller modelled zeta potentials to be further away from the mineral surface, to match the data.

The 0-plane was defined as the hydrolysis layer located 2.3 Å from the calcite surface defined by the surface of calcium ions, at which only protonation/deprotonation of surface sites took place [START_REF] Heberling | Reactivity of the calcite-water-interface, from molecular scale processes to geochemical engineering[END_REF].

The surface charge at 0-plane was found to be negative regardless of the concentration and composition of the adjacent solution and pH range of 5-11 tested in our study, similar to the results reported by Heberling et al. [START_REF] Heberling | Reactivity of the calcite-water-interface, from molecular scale processes to geochemical engineering[END_REF]. On the other hand, the β-plane was set at a distance of 4.6 Å from the surface [START_REF] Heberling | Reactivity of the calcite-water-interface, from molecular scale processes to geochemical engineering[END_REF], and we assumed that all the specific adsorption reactions between the surface sites and salt ions from the bulk electrolyte only occurred at the OHP where ions are adsorbed as outersphere complexes. 236

Ca CO 3 Ca CO 3 Ca CO 3 OH 2 OH OH 2 +x -x x -1 Na + -x Ca CO 3 Ca CO 3 Ca CO 3 x -1 OH +x OH 2 +x OH 2 +x -x -x -x -x Ca 2+ Mg 2+ CO 3 2- SO 4 2- Cl - Ca OH 2 CO 3 Ca 2+ -x +x Ca 2+ HCO 3 - Mg 2+ Mg 2+ Calcite Electric potential, f Distance f 0 z 0 - plane b - plane C SL Figure 2.
Illustration of the Basic Stern Model used in this study. The 0-plane corresponds to the mineral surface where only protonation/deprotonation reactions take place. We used fractional surface charge x = 0.5 for all surface sites. The Stern Layer is confined between the 0-plane and b-plane, with the latter corresponding to the OHP where all adsorption reactions take place and where the zeta potential (z) is defined. CSL is the Stern layer capacitance.

Adsorption reactions and equilibrium constants

All modelled adsorption reactions between calcite surface sites and ions are listed in Table 2. The original set of adsorption reactions proposed by Heberling et al. [START_REF] Heberling | Reactivity of the calcite-water-interface, from molecular scale processes to geochemical engineering[END_REF][START_REF] Heberling | Structure and reactivity of the calcite-water interface[END_REF] and Li et al. [START_REF] Li | Influence of surface conductivity on the apparent zeta potential of calcite[END_REF] was extended to include the reactions between the surface sites and Mg 2+ and SO4 2-(R9 -R11 in Table 2). The equilibrium constant of R1 that corresponds to deprotonation of hydrated calcium surface sites was fixed to the value of LogK1 = -0.5 to yield the relative abundance of positive and negative surface complexes at 0-plane according to Heberling et al. [START_REF] Heberling | Reactivity of the calcite-water-interface, from molecular scale processes to geochemical engineering[END_REF].

Among published studies on SCM of calcite, there is no consensus on whether protonation of > CO 3 -0.5 surface sites on 0-plane should be considered. Some of these studies included the reaction > CO 3 -0.5 + H + ⟺> CO 3 H +0.5 (e.g., [START_REF] Heberling | Structure and reactivity of the calcite-water interface[END_REF][START_REF] Wolthers | The surface chemistry of divalent metal carbonate minerals; a critical assessment of surface charge and potential data using the charge distribution multi-site ion complexation model[END_REF][START_REF] Pokrovsky | Surface Chemistry and Dissolution Kinetics of Divalent Metal Carbonates[END_REF]), while other studies (e.g., [START_REF] Heberling | Reactivity of the calcite-water-interface, from molecular scale processes to geochemical engineering[END_REF][START_REF] Li | Influence of surface conductivity on the apparent zeta potential of calcite[END_REF]) disregarded it. We could not find any published justification for inclusion of > CO 3 -0.5 protonation based on the experimental data or theoretical analysis. However, there have been several publications in support of excluding the protonation of the > CO 3 -0.5 surface sites. One theoretical study on hydration of (1 0 4) calcite surface [START_REF] Villegas-Jiménez | Theoretical Insights into the Hydrated (10.4) Calcite Surface: Structure, Energetics, and Bonding Relationships[END_REF] demonstrated that neutral H2O molecules preferentially 'chemisorbed' into > Ca +0.5 to form > CaOH 2 +0.5 site, and existence of > CO 3 H +0.5 surface sites was attributed to hydration of carbonate surface groups (> CO 3 OH 2 -0.5 ) but not direct protonation of the > CO 3 -0.5 sites. Therefore, in the absence of > CO 3 OH 2 -0.5 there could be no > CO 3 H +0.5 sites resulting from subsequent dissociation of water molecules > CO 3 OH 2 -0.5 ⟺> CO 3 H +0.5 + OH -. These conclusions were supported by a variety of experimental measurements reported in the paper. Another numerical study [START_REF] Kurganskaya | Carbonate Dissolution Mechanisms in the Presence of Electrolytes Revealed by Grand Canonical and Kinetic Monte Carlo Modeling[END_REF] also predicted low fraction of protonated > CO 3 H +0.5 sites (<20%) on a variety of calcite mineral geometries for pH between 5 and 11, which is consistent with pH values tested with our model. For pH greater than 7, the fraction of protonated carbonate sites decreased to below 5% according to this study. Moreover, inclusion of the > CO 3 -0.5 surface site protonation would result in increased pH under simulated equilibrium conditions, thus mismatching the experimentally measured values. Effect of > CO 3 -0.5 protonation on equilibrium pH will be discussed in more detail in Section 3. For these reasons and following Heberling et al., [START_REF] Heberling | Reactivity of the calcite-water-interface, from molecular scale processes to geochemical engineering[END_REF] we did not consider hydration of carbonate surface sites in our model. a -Song et al. [START_REF] Song | Characterizing the Influence of Organic Carboxylic Acids and Inorganic Silica Impurities on the Surface Charge of Natural Carbonates Using an Extended Surface Complexation Model[END_REF].

b -Heberling et al. [START_REF] Heberling | Reactivity of the calcite-water-interface, from molecular scale processes to geochemical engineering[END_REF].

conly one value for the equilibrium constant for R3 was found in the literature (0.45), hence we considered LogK min to be equal to that of R10 due to similarity.

d -Li et al. [START_REF] Li | Influence of surface conductivity on the apparent zeta potential of calcite[END_REF].

eonly one value for the equilibrium constant for R8 was found in the literature (1.68), hence we considered LogK max to be equal to that of R7 due to similarity.

fcombined from two reactions reported in Pokrovsky and Schott [START_REF] Pokrovsky | Kinetics and Mechanism of Dolomite Dissolution in Neutral to Alkaline Solutions Revisited[END_REF].

g -Qiao et al. [START_REF] Qiao | Modeling low salinity waterflooding in mineralogically different carbonates[END_REF].

hcombined from two reactions reported in Pokrovsky et al. [START_REF] Pokrovsky | Surface Speciation Models of Calcite and Dolomite/Aqueous Solution Interfaces and Their Spectroscopic Evaluation[END_REF].

ionly one value for the equilibrium constant for R8 was found in the literature (0.56), hence we considered LogK max to be equal to that of R7 due to similarity. jonly one value for the equilibrium constant for R9 was found in the literature (1.66), hence we considered LogK max to be equal to that of R4 due to similarity.

k -Qiao et al. [START_REF] Qiao | A Mechanistic Model for Wettability Alteration by Chemically Tuned Waterflooding in Carbonate Reservoirs[END_REF].

lthis is a new additional reaction introduced in our model with no values for the equilibrium constant found in the literature, hence we assumed both LogK min and LogK max to be equal to those of R7 due to similarity.

For reactions R2 -R11 we identified ranges of possible variation of equilibrium constants. The ranges were based on previously published values thus defining the minimum and maximum possible equilibrium constant for each reaction, LogK min and LogK max , respectively. Note, that due to lack of published information on some equilibrium constants, their maximum and minimum values were assumed to be similar to those of other adsorption reactions or calculated from published values of consecutive reactions (see caption of Table 2 for details).

We used the Pitzer theory (pitzer.dat database of PREEQC [START_REF] Pitzer | Thermodynamics of electrolytes. I. Theoretical basis and general equations[END_REF][START_REF] Harvie | The prediction of mineral solubilities in natural waters: The Na-K-Mg-Ca-H-Cl-SO4-OH-HCO3-CO3-CO2-H2O system to high ionic strengths at 25°C[END_REF]) to calculate the ion activity coefficients for all complexation reactions that take place in the bulk solution at high ionic strength systems (> 0.5 M). For ionic strengths below 0.5 M it was found that there was no noticeable difference in the modelled results obtained using either the Debye-Huckel (phreeqc.dat database of PHREEQC) or Pitzer theory to calculate the ion activity coefficients. Therefore, the Pitzer theory was used throughout the entire range of tested ionic strengths.

BSMcapacitance of the Stern layer and electric potential distribution within EDL

The electric potential distribution within the Stern layer, between 0-plane and β-plane, is described by a linear variation, similar to the concept of two parallel plates capacitor which it can be expressed as follows:

𝜎 0 = 𝐶 𝑆𝐿 (𝜙 0 -𝜙 𝛽 ), (1)
where 𝜎 0 is the surface charge density at 0-plane (C•m -2 ), 𝐶 𝑆𝐿 is the specific Stern layer capacitance (or Stern layer capacitance for simplicity) between 0-plane and β-plane (F•m -2 ), 𝜙 0 is the electric potential at 0-plane (V) and 𝜙 𝛽 is the electric potential at β-plane (V). The Stern layer capacitance can also be expressed as:

𝐶 𝑆𝐿 = 𝜀 o 𝜀 𝑟 𝑥 , (2)
where 𝜀 o is vacuum permittivity (F•m -1 ), 𝜀 𝑟 is relative permittivity of the Stern layer and x is the distance between 0-plane and β-plane or the thickness of the Stern layer (m). Beyond the β-plane, in the diffuse layer, the electrical potential distribution is described using the Gouy-Chapman theory based on the Poisson-Boltzmann equation. The electrical potential magnitude decreases exponentially with the distance from the mineral surface and ion concentrations follow a Boltzmann distribution.

Two-step optimization

Our model was developed through a two-step optimization process applied to experimental zeta potential results obtained with Estaillades and Ketton rock samples saturated with equilibrated NaCl solutions of ionic strengths between 0.05 M and 5 M.

The first-step optimization aimed at determining the equilibrium constants of R2-R11, while the log 𝐾 R1 and the Stern layer capacitance (CSL) were kept constant and equal to -0.5 and 1.24 F/m 2 , respectively, in accordance with the study of Li et al. [START_REF] Li | Influence of surface conductivity on the apparent zeta potential of calcite[END_REF]. The optimization of equilibrium constants of the surface complexation reactions and Stern layer capacitance was a necessary step, as previous studies considered either low salinity (e.g., [START_REF] Li | Influence of surface conductivity on the apparent zeta potential of calcite[END_REF]) or moderate salinity solutions (e.g., seawater model by Song et al. [START_REF] Song | Characterizing the Influence of Organic Carboxylic Acids and Inorganic Silica Impurities on the Surface Charge of Natural Carbonates Using an Extended Surface Complexation Model[END_REF]) using different sets of the parameters for each ionic strength, thus making their models suitable to very specific experimental conditions. To the best of our knowledge there has been no study that tried SCM for a range of salinities, and especially for ionic strengths above 0.5 M which are typical for many subsurface settings. The first-step optimization was conducted in the following manner:

 Concentration of all ions (Na + , Ca 2+ , Mg 2+ , Cl -and SO 4 2-) and pH reported by Al Mahrouqi et al. [START_REF] Al Mahrouqi | Zeta potential of artificial and natural calcite in aqueous solution[END_REF] were used as input parameters. The concentration of carbon related ions (HCO3 -, CO3 2-), termed C(4) ions in PHREEQC, was not reported in the corresponding paper, therefore it was calculated from the charge balance equation. It was also assumed that the calculated C(4) concentration could exceed that of Ca 2+ and C(4) ions that originate from the dissolved atmospheric CO2 in aqueous solution at equilibrium with air. This assumption is consistent with the experimentally confirmed dissolution of CaCO3 during the equilibration process that resulted in 10 -3 M of dissolved Ca 2+ , and hence the same concentration of C(4) [START_REF] Al Mahrouqi | Zeta potential of artificial and natural calcite in aqueous solution[END_REF][START_REF] Alroudhan | Zeta potential of intact natural limestone: Impact of potential-determining ions Ca, Mg and SO4, Colloids Surfaces A Physicochem[END_REF], which is higher than the equilibrium concentration of 1.5×10 -5 M of C( 4) in solution at equilibrium with air [START_REF] Li | Influence of surface conductivity on the apparent zeta potential of calcite[END_REF].



The dissolution and precipitation of calcite was switched off at this optimization step by implementing PREEQC code 'Calcite 0 0' in 'equilibrium_phases' section [START_REF] Heberling | Reactivity of the calcite-water-interface, from molecular scale processes to geochemical engineering[END_REF][START_REF] Parkhurst | Description of Input and Examples for PHREEQC Version 3-A Computer Program for Speciation, Batch-Reaction, One-Dimensional Transport, and Inverse Geochemical Calculations[END_REF].



The R2 -R11 equilibrium constants were optimized by minimizing the objective function, f, that defines the normalized difference between the observed and the simulated zeta potential and pH:

𝑓 = ∑ ( 𝜁 𝑜𝑏𝑠 -𝜁 𝑠𝑖𝑚 𝛿 𝜁 ) 2 + ∑ ( pH 𝑜𝑏𝑠 -pH 𝑠𝑖𝑚 𝛿 pH ) 2 , ( 3
)
where 𝜁 𝑜𝑏𝑠 is the observed zeta potential (mV), 𝜁 𝑠𝑖𝑚 is the simulated zeta potential (mV), 𝛿 𝜁 is the experimental uncertainty of the zeta potential (mV), pH 𝑜𝑏𝑠 is the observed pH, pH 𝑠𝑖𝑚 is the simulated pH and 𝛿 pH is the experimental uncertainty of pH. Note that pH was included in the objective function since the computed value for a given set of equilibrium constants could be different from the input (experimental) value. Therefore, both directly measured properties (zeta potential and pH) were included in Eq. 3. The optimization was conducted by combining the PHREEQC geochemical code with an optimization software, PEST, in which the Gauss Marquardt Levenberg method was implemented.



The optimization was repeated three times, with the initial values of all equilibrium constants being either i) the minimum found in literature (LogK min ), ii) the maximum reported in the literature (LogK max ), or iii) the median within the tested range (see more details in Table 2 caption). This repetitive optimization procedure was required to find the global minimum of 𝑓, which corresponded to the best match between the experimentally measured and computed zeta potential and pH, so that the equilibrium constants that corresponded to the smallest 𝑓 were taken forward.



The results of the first step optimization procedure are shown in Table 2 (LogK opt column).

It was found that the calculated zeta potentials of Estaillades and Ketton were positive except for Estaillades at ionic strength below 1 M. In contrast, the experimental results showed that the zeta potential of both rock samples remained negative at low ionic strength to become less negative/positive at high ionic strength (Figure 3a). Moreover, the calculated Ketton zeta potential was more positive compared with that of Estaillades, while the experimental data showed an opposite trend. The more positive zeta potential measured in Estaillades sample could not be explained by different rates of calcite dissolution as the equilibrium Ca 2+ concentration was found to be higher in Ketton sample (Figure 9a in Al Mahrouqi et al. [START_REF] Al Mahrouqi | Zeta potential of artificial and natural calcite in aqueous solution[END_REF]), which would imply a more positive zeta potential, consistent with our model.

However, a considerably higher concentration of SO4 2-was also reported in experiments with Ketton (Figure 9b in Al Mahrouqi et al. [START_REF] Al Mahrouqi | Zeta potential of artificial and natural calcite in aqueous solution[END_REF]), and more negative zeta potential compared with that of Estaillades sample was attributed to it. The authors demonstrated that the zeta potential was not very sensitive to SO4 2-concentration (Figure 8b in Al-Mahrouqi et al. [START_REF] Al Mahrouqi | Zeta potential of artificial and natural calcite in aqueous solution[END_REF]), as long as concentration of all other ions remained unchanged thus implying that high sensitivity to SO4 2-only took place when sulfate content was less or equal to that of the divalent cations. However, the paper did not consider a possibility of Ca 2+ ions acting as bridges to SO4 2-ions [START_REF] Nachbaur | Electrokinetic Properties which Control the Coagulation of Silicate Cement Suspensions during Early Age Hydration[END_REF][START_REF] Austad | Is wettability alteration of carbonates by seawater caused by rock dissolution?[END_REF][START_REF] Elakneswaran | Electrokinetic potential of hydrated cement in relation to adsorption of chlorides[END_REF] thus enhancing sulfate adsorption when elevated concentration of Ca 2+ is found in bulk solution. Such bridging ability of Ca 2+ implies that higher concertation of Ca 2+ (and/or Mg 2+ ) combined with higher concentration of SO4 2-(as was observed by Al Mahrouqi et al. [START_REF] Al Mahrouqi | Zeta potential of artificial and natural calcite in aqueous solution[END_REF]) should result in overall higher adsorption of sulfate and more negative zeta potential. Keeping in mind that dissolved SO4 2-could only originate from undetected minerals such as anhydrite (CaSO4) or epsomite (MgSO4•7H2O), a higher concentration of sulfate ions would always be accompanied by a higher equilibrium concentration of Ca 2+ and/or Mg 2+ as it was indeed reported by Al-Mahrouqi et al. [START_REF] Al Mahrouqi | Zeta potential of artificial and natural calcite in aqueous solution[END_REF]. Therefore, we assumed that the reactivity of calcite to SO4 2-was considerably higher for rock samples that leached more sulfate and the constraint on LogK max for R10 was removed, which resulted in a substantially better match to the experimental data (Figure 3b) and the corresponding value of the optimized equilibrium constant, LogK opt , of 11.9 (Table 2). We expect that including Ca 2+ -SO4 2-ion bridging in the surface complexation reactions would lower LogK opt . Such inclusion implemented in SCM, would require currently unavailable experimental data to validate the model and therefore, data acquisition and the corresponding SCM adjustments will be conducted in a follow-up study. After increasing the LogK opt for R10 to the value of 11.9, the match to the experimental results significantly improved as shown in Figure 3b for ionic strength ≤ 2 M. However, the simulated zeta potential at ionic strength greater than 2 M remained significantly more positive compared with the observed values, and for the Ketton sample it became positive as opposed to the negative zeta potential obtained from the experiments. Therefore, additional modification to the model was made.

The second step optimization was initiated by fixing the optimal values of all equilibrium constants obtained from the first step optimization. Then, using the same objective function as in the 1 st optimization step (Eq. 3), the 2 nd step optimization was implemented for each rock sample and each concentration of NaCl solution using the following variables:



The concentration of Ca 2+ and SO4 2-was allowed to vary within the reported experimental uncertainty.



A variable capacitance of the Stern layer was assumed and allowed to vary between 0.2 F•m - 2 and 1.4 F•m -2 (see more detailed discussion in section 2.2.5).

The resulting from the 2 nd optimization step pH and concentration of SO4 2-and Mg 2+ remained the same compared with the results of the 1 st optimization step. However, to obtain a better match to the experimentally measured zeta potential required a non-monotonic change of the optimized Stern layer capacitance (Figure 5a) and increased concentration of Ca 2+ (see the yellowred shift for Estaillades and greenblue shift for Ketton in Figure 5b). The computed zeta potential was found to be in a good agreement with the experimental data after the variable capacitance was implemented (see the yellowred shift for Estaillades and greenblue shift for Ketton in Figure 6). , [START_REF] Al Mahrouqi | Zeta potential of artificial and natural calcite in aqueous solution[END_REF]. The yellowred and greenblue shift represent manual adjustment of pCa to obtain a better match to the computed zeta potential in Fig. 6.

Figure 6. Calculated zeta potential after the first and the second step optimization. Black symbols correspond to the experimental data of Al Mahrouqi et al. [START_REF] Al Mahrouqi | Zeta potential of artificial and natural calcite in aqueous solution[END_REF]. The yellowred (Estaillades) and greenblue (Ketton) shifts demonstrate a significant improvement of the match after applying variable Stern layer capacitance.

Variable Stern layer capacitance

As mentioned in the previous section, in order to obtain a good match with the experimentally measured zeta potential it was required to allow the Stern layer capacitance (CSL) to vary and to decrease globally. Following the definition of CSL (Eq. 2), its optimized concentration dependence (Figure 5a) can be explained by either: i) constant distance between the 0-plane and b-plane, x, and decreasing relative permittivity of the Stern layer, r, with increasing NaCl concentration; ii) constant relative permittivity and increasing with salinity x; iii) increasing with salinity x and decreasing r; or iv) decreasing x and r with increasing salinity so that the rate of decrease of r exceeds that of x.

Previously published studies have suggested that hydration diameter of cations decreases considerably with increasing salinity, hence x might become smaller as smaller in diameter ions may have a closer distance of approach to the mineral surface (see for example, results obtained from an analytical model for Mg 2+ in Afanas'ev and Ustinov [START_REF] Afanas'ev | Hydration numbers and the state of water in hydration spheres of magnesium chloride and magnesium sulfate solutions[END_REF] and for Na + in Afanasiev et al. [START_REF] Afanasiev | State of Hydration Shells of Sodium Chloride in Aqueous Solutions in a Wide Concentration Range at 273[END_REF] using the same approach, both verified by various experimental data) thus restricting the possible explanation for the decreasing CSL to option iv). This assumption is also supported by Brown et al. [START_REF] Brown | Effect of Electrolyte Concentration on the Stern Layer Thickness at a Charged Interface[END_REF] who suggested a substantial compression of the Stern layer with increasing electrolyte salinity, which was interpreted from experimental observations including EPM. In contrast, some molecular dynamics studies (e.g., [START_REF] Bourg | Molecular dynamics simulations of the electrical double layer on smectite surfaces contacting concentrated mixed electrolyte (NaCl-CaCl2) solutions[END_REF]) suggest that there is no change in x with increasing salt concentration due to the presence of the hydrolysis layer. At the same time, an electric field within the Stern layer that is exerted on polar water molecules results in its polarization and may lead to a substantial decrease in r with increasing salinity (e.g., [START_REF] Yeh | Dielectric constant of water at high electric fields: Molecular dynamics study[END_REF]), consistent with options i), iii) and iv). In the light of the above arguments, we assumed option i) in our model, implying constant x and decreasing with salinity r. Note, that we tested hypothesis iv) and found that CSL is significantly less sensitive to the variation of x (even assuming the maximum possible range between fully hydrated Na + radius of 4.5

Å and crystallographic Na + radius of 0.5 Å at 5 M [START_REF] Kielland | Individual Activity Coefficients of Ions in Aqueous Solutions[END_REF]) compared with allowed variation of r (between the value that corresponds to diluted electrolyte, ~80 at a temperature of 20°C and that of structured water in the Stern layer, ~6 [START_REF] Heberling | Structure and reactivity of the calcite-water interface[END_REF]), but this sensitivity analysis is not presented here.

For simplicity we used a linear regression through optimized values of CSL as shown in Figure 7:

𝐶 𝑆𝐿 = -0.2232 × 𝐼𝑆 + 1.357 [𝐹 • 𝑚 -2 ] ( 4
)
Where IS is the ionic strength of the solution of interest (M). To verify the regression is physically meaningful we assumed a constant distance between 0-plane and β-plane of x = 2.3 Å [START_REF] Heberling | Reactivity of the calcite-water-interface, from molecular scale processes to geochemical engineering[END_REF] and used a published equation for the Stern layer relative permittivity [START_REF] Conway | Ionic hydration in chemistry and biophysics[END_REF][START_REF] Hiemstra | On the relationship between charge distribution, surface hydration, and the structure of the interface of metal hydroxides[END_REF]:

𝜀 𝑟 = 𝜀 𝑧 -𝜀 𝑠 1 + 𝑏 (- 𝑑𝜙 𝑥 ) 2 + 𝜀 𝑠 (5)
where 𝜀 𝑧 is the relative permittivity of the bulk electrolyte (at a given ionic strength and temperature), 𝜀 𝑠 is the relative permittivity near the mineral surface (~6 [START_REF] Hiemstra | On the relationship between charge distribution, surface hydration, and the structure of the interface of metal hydroxides[END_REF]), 𝑑𝜙 is the electrical potential difference (V) between 0-plane with β-plane calculated by our model, and b = 1.2 × 10 -17 m•V -1 is a constant [START_REF] Hiemstra | On the relationship between charge distribution, surface hydration, and the structure of the interface of metal hydroxides[END_REF]. The values of CSL calculated using r computed from Eq. 5, and x = 2.3 Å are plotted in Figure 7 and are in good agreement with the optimized values thus confirming the validity of our approach. Note, that 𝑑𝜙 for Eq.5 were computed by PHREEQC using the linear regression of 𝐶 𝑆𝐿 (blue line in Figure 7) and provided an excellent match to the zeta potentials, so that non-linear regressions for 𝐶 𝑆𝐿 were not tried as this was not the focus of this study.

Calcite surface complexation model implementation

To use our model for predicting the zeta potential of calcite in contact with aqueous solutions requires identification of input parameters and model options that accurately replicated the reported experimental conditions. Therefore, the following steps should be taken to obtain an accurate prediction of the zeta potential:

 Define the saturation index and the amount (typically, 10 mol to 20 mol to allow sufficient amount of calcite to equilibrate with the solution of interest) in 'EQUILIBRIUM_PHASES' section (see Appendix A). A non-zero amount should be chosen to replicate chemical equilibration between calcite and water of the experiment of interest, while the saturation index should be adjusted so that the computed Ca 2+ concentration matches the experimental value.



Define input concentrations of all ions in the modelled solution consistent with the reported values. Note, that the input concentration of Ca 2+ should be kept at zero and dissolution/precipitation of calcite switched on, when simulating the zeta potential experiments under equilibrium conditions, in which case content of Ca 2+ in the solution is computed by PHREEQC and cross-compared with the experimental data. Otherwise, the measured concentration of Ca 2+ should be used while switching off the calcite dissolution option (see examples in the Appendix A).



Use the optimized equilibrium constants for R1-R11 (Table 2).



For a given ionic strength, calculate the Stern layer capacitance using the linear regression equation (Eq. 4) and allow PHREEQC to compute the zeta potential. To validate the result cross-compare the simulated pH and concentrations of all ions against the measured values.

Model validation, results and discussion

Our model described in previous sections was tested against the experimentally measured zeta potential in intact carbonate core samples saturated with aqueous solutions reported by Al Mahrouqi et al. [START_REF] Al Mahrouqi | Zeta potential of artificial and natural calcite in aqueous solution[END_REF], Jackson et al. [START_REF] Jackson | Zeta potential in oil-water-carbonate systems and its impact on oil recovery during controlled salinity water-flooding[END_REF] and Li et al. [START_REF] Li | Influence of surface conductivity on the apparent zeta potential of calcite[END_REF]. The developed model simulated the zeta potential from SPM conducted on three different intact carbonate samples saturated with equilibrated NaCl solutions [START_REF] Al Mahrouqi | Zeta potential of artificial and natural calcite in aqueous solution[END_REF], one intact carbonate sample saturated with equilibrated artificial solutions [START_REF] Jackson | Zeta potential in oil-water-carbonate systems and its impact on oil recovery during controlled salinity water-flooding[END_REF], and crushed Iceland spar sample saturated with unequilibrated NaCl solutions [START_REF] Li | Influence of surface conductivity on the apparent zeta potential of calcite[END_REF].

3.1.

Estaillades and Ketton samples saturated with equilibrated NaCl solution [START_REF] Al Mahrouqi | Zeta potential of artificial and natural calcite in aqueous solution[END_REF] We began modelling the zeta potential with defining the input equilibrium concentration of Na + , Cl -, Mg 2+ , SO4 2-as reported in the paper. The initial input concentration of Ca 2+ was set to zero M, partial CO2 pressure is set to 10 -3.44 atm (consistent with atmospheric CO2), and dissolution of calcite is switched on (zero default value of the saturation index and 20 mol amount) to replicate the equilibrium experimental conditions. Applying the specific capacitance values calculated using Eq. 4 for each value of tested ionic strength our model computed the zeta potential and the resulting equilibrium concentration of Ca 2+ , and the results are shown in Figures 8a and8b (note, that simulated pH and concentration of all ions except Ca 2+ were found to be identical to the experimental values used as input parameters). It can be seen from Figures 8a and8b that simulated results for the Ketton sample underestimate the concentration of Ca 2+ and consequently overestimate the negative zeta potential. To be consistent with the experimental data, the calcite dissolution rate for Ketton was enhanced by setting the saturation index to 0.8, so that the modelled equilibrium concentration of Ca 2+ became comparable with the measured one, which resulted in an excellent correlation between the modelled and measured zeta potential (Figures 8c and8d).

Portland sample saturated with equilibrated NaCl solution [12]

Our model was then applied to simulate the zeta potential obtained from Portland limestone. To evaluate the equilibrium concentration of Ca 2+ dissolved from the Portland sample we implemented a similar procedure by comparing the computed and measured equilibrium concentration of Ca 2+ , which was reported to be slightly lower than that for the Ketton sample (compare diamonds and squares in Figure 10a of Al Mahrouqi et al. [START_REF] Al Mahrouqi | Zeta potential of artificial and natural calcite in aqueous solution[END_REF]). A good agreement between the computed and measured Ca 2+ concentration corresponded to the saturation index of 0.6. However, a significantly higher, compared with either Estaillades or Ketton, concentration of SO4 2-was reported for equilibrated solutions in contact with Portland (compare diamonds with squares in Figure 10b of Al Mahrouqi et al. [START_REF] Al Mahrouqi | Zeta potential of artificial and natural calcite in aqueous solution[END_REF]), while the reported concentration of Mg 2+ in all rock samples was the same (pMg = 4.19). Since SO4 2-in the equilibrated solution could only originate from dissolution of sulfate containing minerals (anhydrite or epsomite), the equilibrium concentration of SO4 2-would not be expected to exceed that of the total concertation of Ca 2+ and Mg 2+ , which was not the case reported for Portland sample. We could not explain this discrepancy, and thus assumed that the reported concentration of Ca 2+ was inaccurate and contained an experimental uncertainty, which resulted in underestimated concentration of the cation. Therefore, we artificially added Ca 2+ to our simulated aqueous solution, so that the concentration of the added cation was equal to the difference in sulfate concentration across the salinity range:

∆[𝑆𝑂 4 2-] = [𝑆𝑂 4 2-] 𝑃𝑜𝑟𝑡𝑙𝑎𝑛𝑑 -[𝑆𝑂 4
2-] 𝐾𝑒𝑡𝑡𝑜𝑛 . Note that the difference in sulfate concentration depends on the salinity as shown in Figure 9c, so that a single value of the saturation index could not be used thus justifying our approach of adding Ca 2+ to the solution artificially. Keeping the saturation index equal to 0.6, the total computed concentration of Ca 2+ increased by the amount equal to ∆[𝑆𝑂 4

2-] (as denoted by the yellowgreen shift in Figure 9b), and we successfully modelled the zeta potential as presented in Figure 9a (the yellowgreen shift in Figure 9a shows the effect on the computed zeta potential made by adding Ca 2+ to the solution).

3.3.

Estaillades sample saturated with artificial natural solutions [START_REF] Jackson | Zeta potential in oil-water-carbonate systems and its impact on oil recovery during controlled salinity water-flooding[END_REF] Jackson et al. [START_REF] Jackson | Zeta potential in oil-water-carbonate systems and its impact on oil recovery during controlled salinity water-flooding[END_REF] measured zeta potential on Estaillades rock sample saturated with three different equilibrated solutions: low salinity 20dSW (20 times diluted seawater), SW (seawater) representing moderate ionic strength of 0.749 M, and 3.537 M high salinity FMB (formation brine), with detailed composition provided in Table 1 in Jackson et al. [START_REF] Jackson | Zeta potential in oil-water-carbonate systems and its impact on oil recovery during controlled salinity water-flooding[END_REF]. The purpose of these experiments was to simulate the conventional and/or inverted low salinity waterflooding used to maximize oil recovery.

To simulate the results, we modelled equilibrium conditions between calcite and the aqueous solutions by using zero initial concentration of Ca 2+ and the same calcite dissolution rate as in Section 3.1. for Estaillades (default zero value of saturation index). The results of modelled zeta potential for all three solutions are presented in Figure 10. Our model could successfully reproduce the measured zeta potential (Figure 10a). However, the modelled equilibrium pH was significantly higher with 20dSW and significantly lower with FMB experimental data (Figure 10b) despite the fact that modelled concentrations of all ionic species were found to be identical to the measured ones. We hypothesize that such discrepancy could result from the experimental protocol reported by Jackson et al. [START_REF] Jackson | Zeta potential in oil-water-carbonate systems and its impact on oil recovery during controlled salinity water-flooding[END_REF] who reported the initial (at time of preparation) pH and ionic concentration of solutions but not the final values established after equilibration with the rock sample. It is worth mentioning, that equilibrium pH computed by our model in all tested aqueous solutions and rock samples matched the experimentally measured values apart from the results presented in Figure 10, where the simulaed equiibrium pH was higher than experimental for 20dSW and SW and lower than experimental for FMB. Inclusion of protonation of > CO 3 -0.5 surface sites without changing the optimized equilibrium constants resulted in increased simulated equilibrium pH compared with the experimental values. Thus, the assumed exclusion of protonation of the > CO 3 -0.5 surface sites was confirmed for all tested parameters except for FMB, for which such reaction might need to be included, and this will tested in a follow-up study.

3.4.

Iceland spar saturated with NaCl at non-equilibrium conditions [START_REF] Li | Influence of surface conductivity on the apparent zeta potential of calcite[END_REF] Finally, our model was tested against the data reported by Li et al. [START_REF] Li | Influence of surface conductivity on the apparent zeta potential of calcite[END_REF] on crushed Iceland spar saturated with 10 -3 M, 10 -2 M, and 5×10 -2 M NaCl not equilibrated with the mineral. Description of the experimental procedure reported by Li et al. [START_REF] Li | Influence of surface conductivity on the apparent zeta potential of calcite[END_REF] did not provide detailed information on the duration of the streaming potential measurements. However, considering the expected high permeability of the crashed samples (comparable to permeability of sandpacks reported by Vinogradov et al. [START_REF] Vinogradov | Zeta potential in sandpacks: Effect of temperature, electrolyte pH, ionic strength and divalent cations, Colloids Surfaces A Physicochem[END_REF]) and low salinity solutions used by Li et al. [START_REF] Li | Influence of surface conductivity on the apparent zeta potential of calcite[END_REF], it was assumed that the reported streaming potential measurements on crushed Iceland spar did not last hundreds of hours required for complete equilibration between the mineral and the tested solutions. Therefore, we assumed that only partial equilibration was reached during the experiments and for that reason we modelled the data reported by Li et al. [START_REF] Li | Influence of surface conductivity on the apparent zeta potential of calcite[END_REF] as obtained at non-equilibrium conditions so that any equilibration of calcite was disabled in the model (zero amount was set, Section 2.3 and Appendix A). However, due to partial equilibration expected in the experiments, we artificially added some non-zero initial concentration of Ca 2+ to the modelled solution. The concentration of dissolved Ca 2+ in 0.05 M NaCl was adjusted to 10 -4.2 M (expressed as pCa = 4.2 in Figure 11b) at pH 9, consistent with observations reported by Alroudhan et al. [START_REF] Alroudhan | Zeta potential of intact natural limestone: Impact of potential-determining ions Ca, Mg and SO4, Colloids Surfaces A Physicochem[END_REF] . Firstly, we assumed that the amount of dissolved (and therefore artificially added to the model) Ca 2+ should depend on pH (lower pH would lead to higher concentration of dissolved Ca 2+ , consistent with higher dissolution rate reported by Anabaraonye et al. [START_REF] Anabaraonye | Brine chemistry effects in calcite dissolution kinetics at reservoir conditions[END_REF] . Furthermore, the pH dependence of pCa in 0.05 M NaCl experiments was extended throughout the entire range of pH (5-11) using a linear slope of ∆𝑝𝐶𝑎 ∆𝑝𝐻 ⁄ = 0.1, consistent with calcite dissolution rate reported by Chou et al. [START_REF] Chou | Comparative study of the kinetics and mechanisms of dissolution of carbonate minerals[END_REF] .

Secondly, we assumed that dissolved Ca 2+ concentration should also depend on NaCl concentration (higher NaCl concentration would lead to higher concentration of dissolved Ca 2+ as reported in Al Mahrouqi et al. [START_REF] Al Mahrouqi | Zeta potential of artificial and natural calcite in aqueous solution[END_REF]). Therefore, for NaCl solutions of 10 -2 M and 10 -3 M the end-points of pCa were moved up the vertical scale by 0.4 units relative to one another to reflect on the salinity dependence of calcite dissolution. The value of 0.4 was approximated from the reported calcite dissolution rate as a function of NaCl concentration [START_REF] Ruiz-Agudo | An atomic force microscopy study of calcite dissolution in saline solutions: The role of magnesium ions[END_REF] and assuming duration of the steaming potential measurement experiments conducted by Li et al. [START_REF] Li | Influence of surface conductivity on the apparent zeta potential of calcite[END_REF] to be approximately 10 hours. Ultimately, ascribing the end-point values for Ca 2+ concentration for each tested NaCl solution on the pH scale (Figure 11b) and applying linear variation of pCa between these end-points, our model successfully predicted the measured zeta potential (Figure 11a). Moreover, our model predicts that Ca 2+ must be added to the modelled solution in order to replicate the experimentally measured zeta potential. Note, that C(4) of the amount equal to the added Ca 2+ must also be added to the solution in addition to the initial concentration of C(4) that reflects the equilibrium content of dissolved CO2. Despite the fact that the computed from calcite dissolution or artificially added Ca 2+ and C(4) to the bulk electrolyte was not reported by Li et al. [START_REF] Li | Influence of surface conductivity on the apparent zeta potential of calcite[END_REF], the adsorption reactions for these ions were included in their SCM. Therefore, we are confident that our approach is consistent with Li et al. [START_REF] Li | Influence of surface conductivity on the apparent zeta potential of calcite[END_REF]. Our results also explain how the experimental data from Li et al. [START_REF] Li | Influence of surface conductivity on the apparent zeta potential of calcite[END_REF] was successfully matched by means of non-zero concentration of Ca 2+ and C(4) in the solution, both of which should increase with decreasing pH. If such pH dependence of Ca 2+ and C(4) concentration was not applied to the model of Li et al. [START_REF] Li | Influence of surface conductivity on the apparent zeta potential of calcite[END_REF], the simulated zeta potential would become more positive with increasing pH, which is an inverted trend relative to that reported in the paper.

Model capabilities, limitations and implication to carbonate-water subsurface settings

The developed robust SCM uses a set of justified and well-defined calcite surface sites and surface reactions' equilibrium constants combined with physically meaningful salinity dependence of the Stern layer capacitance. The model is capable to accurately predict zeta potential of a variety of carbonate-aqueous solution systems as long as calcite is the dominating mineral comprising natural rocks and the aqueous solution is dominated by NaCl. Moreover, an additional requirement for predictive capability of our model relates to the naturally occurring concentration of SO4 2-in the solution, which has to reflect dissolution of sulfate containing minerals such as anhydrite or epsomite over geological deposition timescale. This requirement implies that the equilibrium concentration of SO4 2-should not exceed the combined concentration of Ca 2+ and Mg 2+ .

The set of equilibrium constants proposed in this paper for ambient conditions (25 o C and 1 atm) has not been tested to model experimental results obtained with aqueous solutions with artificially added SO4 2-beyond natural occurrence, so that some amendments could be required to model such experiments. Moreover, we have not considered elevated temperature experimental conditions, which would require additional adjustments to the model: i) the equilibrium constants should be adjusted for elevated temperature by using for example the Van't Hoff equation; ii) the relative permittivity of the Stern layer should be recalculated; iii) the distance between the 0-plane and the b-plane should be reduced to account for reduced hydrated diameter of cations at elevated temperature (consistent with [START_REF] Vinogradov | Zeta potential in sandpacks: Effect of temperature, electrolyte pH, ionic strength and divalent cations, Colloids Surfaces A Physicochem[END_REF][START_REF] Kuyucak | Temperature dependence of conductivity in electrolyte solutions and ionic channels of biological membranes[END_REF]); iv) the salinity dependence of the Stern layer capacitance should be amended taking into account ii) and iii) so that a new temperature specific regression is used.

Capabilities of our model have a broad range of applications since the approach reported here has demonstrated to provide an accurate prediction of the zeta potential of natural carbonate rock in contact with aqueous solutions, subject to the above conditions. Our model works best for equilibrium experimental conditions as long as concentration of all major ionic species is known. In this sense, our model can significantly improve our understanding of streaming potential measurements and associated flows in shallow aquifers where temperature and pressure are low and therefore, laboratory measurements of pH and concentration of all ions of the solution at aquifer conditions are straight forward. Laboratory experiments should, in this case, assure establishment of full equilibrium between rock and the solution of interest and fluid samples should be analyzed for molar concentration of all ions to be used as input parameters of the model. The modelled zeta potential can then be used to predict flow patterns in critical zones (e.g., [START_REF] Jougnot | Modeling Streaming Potential in Porous and Fractured Media, Description and Benefits of the Effective Excess Charge Density Approach[END_REF]), permeability heterogeneities (e.g., [START_REF] Jackson | Spontaneous Potentials in Hydrocarbon Reservoirs During Waterflooding: Application to Water-Front Monitoring[END_REF]) or even serve as an early warning of saline intrusion (e.g., [START_REF] Graham | Self-Potential as a Predictor of Seawater Intrusion in Coastal Groundwater Boreholes[END_REF][START_REF] Macallister | Remote Detection of Saline Intrusion in a Coastal Aquifer Using Borehole Measurements of Self-Potential[END_REF]).

Our model is useful in situations where the exact concentration of all ions is unknown. For example, the model can accurately assess the initial estimate of CO2 geological storage efficiency. The model would require input from the experimental data, given that laboratory experiments are conducted to establish full chemical equilibrium between carbonate rock and a CO2-saturated aqueous solution at the target formation conditions of temperature and pressure. Such experiment should report at least one measured concentration (of constituent ions such as Ca 2+ , Mg 2+ , SO4 2-or protons, pH) to validate the simulated by our model concentrations (similar to the experiments reported by Li et al. [START_REF] Li | The pH of CO2-saturated aqueous NaCl and NaHCO3 solutions at temperatures between 308 K and 373 K at pressures up to 15 MPa[END_REF]). The modelled zeta potential then can be used to interpret the wetting state using Derjaguin-Landau-Verwey-Overbeek (DLVO) theory (e.g., [START_REF] Tokunaga | DLVO-Based Estimates of Adsorbed Water Film Thicknesses in Geologic CO2 Reservoirs[END_REF]) and the resulting residual trapping of the gas.

Conclusions

We report a robust surface complexation model of calcite-water interfaces. The model was developed adopting a two-step optimization, in which the equilibrium constants of surface reactions were initially optimized to match two experimental datasets, while the ionic strength dependence of the Stern layer capacitance was obtained in the second optimization step. The model was applied to successfully predict multiple experimental datasets with an excellent quality of match, and the modelling results demonstrate that:

 The unique set of the optimized equilibrium constants can be used universally to simulate the calcite-water zeta potentials obtained from the streaming potential measurements for all tested systems and conditions.

 The Stern layer capacitance should decrease with increasing salinity to replicate high salinity zeta potentials; the range of varying capacitances was found to be consistent with analytically predicted values [START_REF] Hiemstra | On the relationship between charge distribution, surface hydration, and the structure of the interface of metal hydroxides[END_REF] and with values used in previously published numerical studies [START_REF] Yeh | Dielectric constant of water at high electric fields: Molecular dynamics study[END_REF].

 Using experimental composition of investigated solutions as the only input, our model accurately predicts zeta potentials of all tested systems and conditions including natural limestones equilibrated with simple salt and complex solutions of ionic strengths between 0.05 M and 5 M, and crushed calcite not equilibrated with NaCl solutions of ionic strength between 0.001 M and 0.01 M.

 To simulate zeta potentials of natural carbonate-water systems at equilibrium conditions, the developed model requires knowledge of concentration of SO4 2-and Mg 2+ that leach from dolomite and anhydrite inclusions, while the equilibrium concentration of Ca 2+ is produced by the model.

 To simulate zeta potentials of calcite in contact with water under non-equilibrium conditions requires knowledge of either measured in real-time Ca 2+ concentration or duration of the simulated experiment.

 Our model is fully predictive given the required input parameters (e.g., ion concentration) are provided. However, additional surface reactions and/or model adjustments might be required to simulate the zeta potentials for carbonate rock in contact with aqueous solutions in which SO4 2- concentration exceeds that of Ca 2+ and Mg 2+ combined, this will be investigated in a follow-up study.

Future work will also aim at acquiring additional experimental data obtained with varying concentration of C(4) and SO4 2-ions under equilibrium conditions to inform the surface complexation model, which will be updated and/or modified to include the hypothesis of ion bridging outside OHP. ----------------------------------------------------------------------------------------------------------------------------------

EQUILIBRIUM_PHASES

Figure 1 .

 1 Figure 1. Electrical potential distribution at the calcite-water interface: a) -BSM (no distinction between IHP and OHP, the slip plane may or may not coincide with OHP); b) -TLM (IHP and OHP are considered separately with different surface complexation reactions taking place at each plane, the slip plane coincides with OHP); c) -QLM (IHP and OHP are considered separately with different surface complexation reactions taking place

Figure 3 .

 3 Figure 3. Calculated zeta potential after the first step optimization a) with the constraint -2.10 ≤ log 𝐾 R10 ≤ 3.30, b) without the constraint -2.10 ≤ log 𝐾 R10 . Black symbols denote the experimental results from AlMahrouqi et al. [12]. The calculated concentration of Ca 2+ and SO4 2-(note, that input pH and concentration of Mg 2+ were constant at 8.3 and pMg = 4.19, respectively, as reported by Al-Mahrouqi et al. [12], pMg = -logCMg with C denoting the ion concentration in M) for both rock samples is shown in Figure 4. The computed solution pH and concentration of all other ions (those that are not shown in Figure 4) remained unchanged and equal to the input values regardless of LogK opt for R10, hence the corresponding plots for these ions are not included in the figure.

Figure 4 .

 4 Figure 4. Computed pCa (a) and pSO4 (b) from the first optimization step. Black symbols represent the experimental data [12]. The results in color correspond to the modelled zeta potential that appears in Figure 3b.

Figure 5 .

 5 Figure 5. Optimized Stern layer capacitance (a) and pCa as a function of NaCl concentration. Black symbols in b) correspond to the experiments of Al Mahrouqi et al., [12]. The yellowred and greenblue shift represent

Figure 7 .

 7 Figure 7. Optimized Stern layer capacitance compared with values obtained using Eq. 5 and approximated by the linear regression fitted to the optimized values in blue (Eq. 4; quality of match to the optimized values of capacitance in blue is R 2 = 0.91) as a function of ionic strength.

Figure 8 .

 8 Figure 8. Experimentally measured (black symbols) and simulated zeta potential (a, c) and pCa (b, d) at equilibrium conditions for Estaillades and Ketton samples saturated with NaCl solutions. The modelled results were obtained assuming identical saturation index for both rock samples toward NaCl solutions (a, b), and higher saturation index for Ketton (0.80) sample relative to the default value (0) of Estaillades (c, d).

Figure 9 .

 9 Figure 9. Measured and modelled zeta potential (a), pCa (b) and SO4 2-(c) of three natural rock samples. The amount of artificially added Ca 2+ to the model of Portland equals the reported difference in concentration of SO4 2-as shown in (c) and explained in Section 3.2. Yellow diamonds represent the modelled zeta potential (a) and pCa (b) assuming saturation index of 0.6 and zero artificially added Ca 2+ . The results of the model with artifically added Ca 2+ are denoted by green diamonds.

Figure 10 .

 10 Figure 10. Simulated zeta potential (a) and pH (b) of three different brine compositions for Estaillades rock sample plotted as a function of ionic strength. The experimental data denoted by black symbols are extracted from Jackson et al. [49].

Figure 11 .

 11 Figure 11. Modelled zeta potential (a) and pCa (b) of crashed Iceland spar samples saturated with three different NaCl solutions under non-equilibrium conditions. The modelled zeta potential is plotted in comparison with the experimental data [11] denoted by black symbols. The slopes and intercepts of the pCa as a function of pH are discussed and defined in Section 3.4.

#example 3 -

 3 simulation of no established equilibrium between calcite and water (Iceland Spar), default saturation index non-zero input concentration of Ca 2+ . 0000631 #equivalent to pCa=4.2 interpreted from Figure 2a (Alroudhan, et al., 2016) assuming 10 hours of partial equilibration during the experiment, which results in pH=9. SAVE SOLUTION 3 -

Table 1 .

 1 Petrophysical and mineralogical properties of three different carbonate rock samples. Ketton and Estaillades were used in the experiments of by Al Mahrouqi et al.[START_REF] Al Mahrouqi | Zeta potential of artificial and natural calcite in aqueous solution[END_REF], Portland was used by Alroudhan et al.

	/rock	Ketton	Estaillades	Portland
	Porosity (%)	23.0 ± 0.5	28.0 ± 0.5	20.0 ± 0.5
	Permeability (Darcy) 1.4 ± 0.4	0.13 ± 0.2	0.005 ± 0.001
	Formation factor (F) 13.87 ± 0.5	12.92 ± 0.5	22.04 ± 0.5
		97 a -calcite	97 a (95 b ) -calcite	96.6 a -calcite
	Composition (%)	3 a magnesium	3 a (4 b ) -magnesium 3.4 a -quartz
			(1 b ) -anhydrite

[START_REF] Alroudhan | Zeta potential of intact natural limestone: Impact of potential-determining ions Ca, Mg and SO4, Colloids Surfaces A Physicochem[END_REF]

.

Table 2 .

 2 Adsorption reactions and their equilibrium constants. Equilibrium constant for R1 was not optimised.

	No	Reactions	LogK min LogK max	LogK opt
		0-plane reactions
	R1	> CaOH 2 +0.5 ⟺> CaOH -0.5 + H +	-0.5	-0.5	-0.5
		β-plane reactions
	R2	> CaOH -0.5 + Na + ⟺> CaOH -0.5 -Na +	0.09 a	0.56 b	0.56
	R3	> CaOH 2 +0.5 + Cl -⟺> CaOH 2 +0.5 -Cl -	-2.10 c	0.45 d	-0.64
	R4	> CaOH -0.5 + Ca 2+ ⟺> CaOH -0.5 -Ca 2+	1.68 b	3.40 e	3.40
	R5	> CaOH 2 +0.5 + HCO 3 -⟺> CaOH 2 +0.5 -HCO 3 -	0.04 b	12.50 f	10.65
	R6	> CaOH 2 +0.5 + CO 3 2-⟺> CaOH 2 +0.5 -CO 3 2-	-7.07 b	6.00 g	-4.59
	R7	> CO 3 -0.5 + Ca 2+ ⟺> CO 3 -0.5 -Ca 2+	1.68 b	3.40 h	3.40
	R8	> CO 3 -0.5 + Na + ⟺> CO 3 -0.5 -Na +	0.56 b	3.40 i	0.56
	R9	> CaOH -0.5 + Mg 2+ ⟺> CaOH -0.5 -Mg 2+	1.66 a	3.40 j	2.81
	R10	> CaOH 2 +0.5 + SO 4 2-⟺> CaOH 2 +0.5 -SO 4 2-	-2.10 g	3.30 k	3.30  11.9
	R11	> CO 3 -0.5 + Mg 2+ ⟺> CO 3 -0.5 -Mg 2+	1.68 l	3.40 l	2.81

Acknowledgments

Miftah Hidayat was supported by the Aberdeen-Curtin PhD studentship. David Vega-Maza is funded by the Spanish Ministry of Science, Innovation and Universities ("Beatriz Galindo Senior" fellowship BEAGAL18/00259). The research work of Philippe Leroy was conducted within the framework of the O-ZNS project which is part of PIVOTS project (financial support provided by the Région Centre-Val de Loire and the French Ministry of Higher Education and Research).

Appendix A -PHREEQC code examples

In examples below, text that follows the # symbol is not part of the code, but comments aimed at clarifying the meaning of program keywords and values EQUILIBRIUM_PHASES #basic description Calcite 0 20 # 0 is the default saturation index; 20 is the available amount of calcite in mol. The saturation index can be either negative or positive to hinder or enhance calcite dissolution, respectively. To prevent any dissolution or precipitation of calcite during the equilibration, the amount should be set to zero. CO2(g) -3.44# defines partial CO2 pressure at equilibrium with atmospheric air equal to 10 -3.44 atm ---------------------------------------------------------------------------------------------------------------------------------- --------------------------------------------------------------------------------------------------------------------------------- ----------------------------------------------------------------------------------------------------------------------------------