A multi-species reactive transport model based on ion-solid phase interaction for saturated cement-based materials

This study established a reactive transport model (RTM) for saturated cement-based materials based on multi-ionic transport through an uncharged, ‘free water ’ and charged, ‘the electrical double layer ’ coupled with a chemical equilibrium model for interaction between pore solution and solid phase. A set of modified Poisson-Nernst-Planck equations, including the effect of temperature, pore solution property, and pore structure changes, have been introduced for multi-ionic transport. A comprehensive method has been illustrated for determining the diffusion coefficient at infinite dilution of the aqueous species present in the pore solution in a cementitious material. Numerical predictions were compared to experimental data to demonstrate the applicability of the RTM. For mortar specimens exposed to NaCl solution and seawater, numerical predictions agree well with experimentally determined Portlandite and carbonate in the phase assemblage, elemental distribution such as chloride, sodium, potassium, magnesium, and sulfur adjusting the initial tortuosity factor.


Introduction
Reinforced concrete materials are widely used for civil infrastructure such as buildings, bridges, tunnels, roads, offshore structures, etc.However, in their service life, reinforced concrete structures located in a harsh environment are severely attacked by the ingress of substances from the service environment [1][2][3][4].Furthermore, research and experience over the past 20 years show that planning, material, and construction costs are often dwarfed by the costs associated with material deterioration, i.e., costs of inspecting, repairing, and maintaining civil infrastructure over its useful life [5].Therefore, numerical methods capable of predicting the changes in concrete materials over the service life help select the cement composition and mix design, optimize the cost for infrastructure, and predict the structure's lifetime.
Several numerical models were developed to predict deterioration, especially for cement-based materials, based on this perspective over the last two decades.Initially, deterioration models were developed only considering multi-ionic transport models [6][7][8][9].Later, multi-ionic transport models were coupled with chemical equilibrium [10][11][12][13].Currently, RTMs have been developed, including a multi-ionic mass transport model coupled with a chemical equilibrium including surface complexation [14,15].Here, the interactions between pore solution and solid phases and the interactions between the C-S-H phase and the surrounding ions of the pore solution were included.In the current RTM, the Poisson-Nernst-Planck (PNP) equation derived using the hybrid mixture theory [16,17] is commonly used for the multi-ionic transport model in porous media.PHREEQC [18][19][20] is adopted to perform chemical equilibrium between pore solution and solid phases, including surface complexation.However, the modelling approaches are restricted by simplifying assumptions and a large number of input parameters concerning the processes associated with the deterioration of cementbased materials.Predictive capabilities are also lacking due to the absence of essential relations between chemical reactions, pore structure changes, and ion mass transport.Without considering these interactions, the predictive capabilities of such modelling approaches will always be limited.
The electrical double layer (EDL) will form on the pore surface of cementitious materials due to interactions between the surface site of the C-S-H phase and the surrounding ions of the pore solution.As a result of the EDL formation, the diffusion of anions and cations will take a different pathway in the pore solution.Most of the anions diffuse through free water space, while the diffusion of cations occurs through the EDL [15,21,22].In the current RTMs [10,[12][13][14], the ion mass transport is performed only through free water space, whereas the ions in the EDL are considered stagnated ions, which means that part of cations is neglected in the ion mass transport calculation.In addition, the ion mass transport in porous materials strongly depends on temperature [23], composition and ionic strength of pore solution [24,25], and pore structure [26,27].Pore solution composition, phase assemblage, and pore structure will continuously vary in a cement-based material due to exposure to long-term environmental actions.The effect of these changes is not taken into account by the traditional PNP equation.Moreover, most approaches accounting for surface complexation assume a fixed specific surface area of the C-S-H phase in the cementitious materials over time [13,14].Due to the reaction between pore solution and solid phases, the formation of new phases, as well as the precipitation and dissolution of solid phases, occur on the pore surface.As a result, the specific surface area of the C-S-H phase on the pore surface continuously changes in cementitious materials.Therefore, the ion absorption capacity of the pore surface by the silanol sites (≡SiOH) varies over time.
In order to accurately predict the durability of cement-based materials exposed to various environmental conditions, RTMs should be developed with the following considerations: i) the flux of ions in both free water and the EDL must be considered for the actual multi-ionic transport calculation through a charged surface, ii) the effect of temperature, pore solution property, and pore structure changes must be accounted for in ion transport through a porous medium, iii) experimental results strongly indicate that C-S-H absorbs part of the ionic species in pore solution on the pore surface [28][29][30][31][32][33][34][35][36].Thus, the exchange reaction between an ion in the pore solution and silanol sites (≡SiOH) needs to be included in the ion absorption on the pore surface in the surface complexation model, and v) for computation of the variation of the specific surface area of the C-S-H phase and the effect of pore structure changes on ion mass transport, a proper micro-pore structure model must be adapted as part of RTMs.

Outline of reactive transport model (RTM)
The present study outlines the framework for a RTM, which can predict the spatial and temporal variation of pore solution composition and solid phases for various cement compositions and different exposure conditions.A schematic description of the phenomena and processes included in the RTM is shown in Fig. 1.The available space in pore distribution for ion mass transport and the schematic illustration of ion mass transport due to interaction between the ions and solid phase is described in Fig. 1a and Fig. 1b, respectively.The underlying framework of the RTM is thereby dealt with within two different time domains, i.e., the early-stage development of cementitious materials during the curing period and the deterioration processes of cementitious materials due to exposure to long-term environmental actions.The modelling part of early-stage material development and material deterioration is illustrated through regions A and B, respectively, in Fig. 1c.The RTM is established for saturated conditions within the present study, including multiple constitutive models, i.e., cement hydration, micropore structure, multi-component ionic transport models through free water and the EDL, and chemical equilibrium models.
During early-stage material formation, the incorporated hydration model determines the degree of cement hydration from input data of the mix design, composition of cement, and curing conditions.Subsequently, results from the hydration model are used to calculate the amount of pure phases (EQ) and solid solution phases (SS), the composition of pore solution, amount of absorbed ions on pore surface, and water content.During the deterioration process, the hydration model determines the degree of hydration for the remaining unhydrated cement with the exposure condition for each time step and the new  equilibrium state of solid phases, the concentration of ions in free water and the EDL, amount of absorbed ions on the pore surface, and water content is computed applying the equilibrium phase theory.The geochemical model PhreeqcRM determines the chemical equilibrium, accounting for pore solution, pure and solid solution phases, and surface complexation in this study.For the computation of the material's microstructure, the total porosity of interlayer, gel, and capillary pores and the distribution of gel and capillary pores are determined from the phase assemblage provided by the geochemical computation.The effective diffusion coefficient of ions is subsequently calculated based on computed pore structure and pore solution property, whereas threshold porosity is also determined with updated pore structure.Finally, multicomponent ion transport through free water and the EDL is established with updated boundary conditions.
In the present study, the model proposed by Parrot and Killoh [37], a set of equations to describe the hydration rate of individual major clinker at each time step, is used for the hydration model.The reliability of the Parrot and Killoh approach for different types of water-cement ratio, blended cement, relative humidity, and temperature is described in [38][39][40][41].The pore structure model proposed by Maekawa et al. [42] is adapted to compute the microstructure in this study.The adaption provides a more reliable modification, adapted to various cement and water-cement ratios considering the dissolution and precipitation of solid phases through chemical equilibrium calculations [42].The details of the ion mass transport model, chemical equilibrium model, and numerical studies of the developed RTM for saturated cement-based material are illustrated in the following sections.

General outline
Under the assumption of saturated conditions, all pores of cementitious materials are assumed to be filled with the pore solution during transport calculation.Thus, in the present study, a moisture transport model is not adapted, and the effect of moisture transport on the flux of ions is not included in the ion mass transport model.An electrical double layer (EDL) will form in cementitious materials due to the interaction between ions in the pore solution and silanol sites (≡SiOH) of the C-S-H phase in the pore surface.The water in the EDL (charged water) and water in the middle of the pore (uncharged water) are defined as EDL water and free water, respectively (see Fig. 1b).The ion diffusion and pathway of ions along an EDL water space are different from free water space.Therefore, the multi-ionic transport model will be developed using a one-dimensional FEM method in free water and the EDL separately.In addition, the effect of temperature, pore solution properties, and pore structure changes are implemented in the multi-ionic transport model.Moreover, it is considered that pores below a certain threshold radius are not accessible for multi-ionic transport (even if the pores are saturated) due to the formation of an EDL and high electrical potentials (see Fig. 1a).In the present study, the proposed model based on threshold radius (r th ) by Yuya Takahashi et al. [43] is used to calculate threshold porosity (φ th ).

Balance equations
This section will discuss the balance equations for the ion constituents in free water and the EDL separately.The mass balance equation for the j th ion constituents in the free water of pore solution is, where φ gc is the sum of gel and capillary porosity, φ th is the threshold porosity, ε f l is the saturation level of free water in the pore volume, ρ f lj is the mass density of the j th ion constituent in the free water, which is related to the mass and volume as ρ f lj = M f lj /V f , where M f lj and V f are the mass of j th ion constituent in the free water and volume of free water, respectively, and v f lj is the absolute velocity of the j th ion constituent in the free water.The property q f lj denotes the mass gain of the j th constituents in free water from all other constituents present in the phase due to chemical reactions.The mass balance equation for the j th ion constituent in the EDL of the pore solution can be written as, where A s is the surface area of the pore occupied by the pore solution, t d is the thickness of the EDL, A th is the surface area of the threshold porosity, ρ d lj is the mass density of the j th ion constituent in the free water, which is related to the mass and volume as ρ d lj = M d lj /V d , where M d lj and V d are the mass of j th ion constituent in the EDL and volume of the EDL, respectively, and v d lj is the absolute velocity of the j th ion constituent in the EDL.The property q d lj denotes the mass gain of the j th ion constituent in the EDL from all other constituents present due to chemical reactions.

Constitutive equations
This section discusses the constitutive equations leading to the governing equations for ion mass transport in free water and the EDL based on the approach presented in [15,16].The cement-based material is considered as a non-deformable rigid porous material with constant porosity within the transport calculations.Moreover, for fully saturated conditions, all pores in the material are considered filled with liquid water.Thus, the following relations are valid in this case: where V s and V p is the volume fraction of solids and total porosity in the material, respectively, φ lr is interlayer porosity, and ε d l is the saturation level of water in the EDL.The following constitutive equation for the diffusion flows v lj, l as derived in Bennethum and Cusham [17], is subsequently adapted for multi-ionic transport.r lj ρ lj v lj,l = − ρ lj ∇μ lj + ρ lj g lj − ρ lj z j E T (3) where r lj is the resistivity tensor for the j th ion constituents, g lj is the gravity for the j th ion constituents, μ lj is the chemical potentials for the j th ion constituents, and E T is the total electrical field.In the constitutive equation for the diffusion flow of ions, the charge character of all ionic constituents of the solution is allowed to affect the diffusion of the considered j th ion constituent, which is included in the last term on the right-hand side.Moreover, the ionic diffusion flow depends on the constitution of the chemical potentials of the diffusing species in the solution.In the present approach, the chemical potential, which was introduced without the electrical field term in Bennethum and Cusham [17], is used, and the activity coefficient is adapted for the influence of the composition of the pore solution as follows: where μ 0 lj is the reference chemical potential for the j th ion constituent in the liquid phase, γ lj is the activity coefficient for the j th ion constituent in the liquid phase, and R m lj is the gas constant for the j th ion constituent in the liquid phase.The constitutive equation for the diffusive flow, i.e., Eq.
(3), is subsequently modified with the introduced chemical potential in Eq. ( 4), where the total electric field intensity is described as the gradient of the total electric intensity potential Φ T as E T = ∇ Φ T .Since the effect of gravity, g lj on the diffusive flow of ions is negligible, the term is omitted in Eq. (5).

Governing equations
The governing equations for the multi-ionic transport in free water S. Sharmilan et al. and the EDL are derived using the introduced balance equation and appropriate constitutive equations.In fully saturated conditions, it is assumed that the moisture in the porous media is in a static condition, and thus, the velocity of the liquid phase can be neglected.Hence, the diffusion velocity of the j th ionic species v lj, l is considered equal to the absolute velocity of the j th ionic species v lj in the porous material.Moreover, it is assumed that the following relations are valid for ion flow in free water and the EDL, respectively; is the diffusion velocity of the j th ionic species in free water and v d lj, l is the diffusion velocity of the j th ionic species in the EDL.To transform the mass flow equation into moles, the mole mass term, M lj of the ionic species is introduced.The mass concentration ρ lj and the mole concentrations C lj are related through ρ lj = M lj C lj and the gas constant, R, which Here, v j is the valence of the j th ionic species and the charge z j is related to the valance by the Faraday constant F C and the mole mass through v j F C = z j M lj .For multi-ionic transport in free water, Eqs. ( 1) and ( 5) are combined via a velocity condition and represented in mole quantities, assuming no chemical reactions during ionic transport occur, i.e., q f lj = 0, where C f lj is the mole concentrations of ions in free water.Similarly, Eqs.
(2) and ( 5) are combined based on the velocity condition and formulated in mole quantities for multi-ionic transport in the EDL, assuming no chemical reactions during ionic transport, i.e., q d lj = 0, where C d lj is the mole concentrations of ions in the EDL.D lj and A lj are the diffusion constant and ionic mobility, respectively, and are identified in the multi-ionic transport equation for free water and EDL as The first term on the right-hand side of Eqs. ( 6) and ( 7) describes the self-diffusion of ions.Here, the term ∂ ln γ j /∂ ln C lj represents the effect of the mixture composition of all ions on the diffusion of individual ions.
The second term accounts for the effect of the electrical potential on the ion diffusion in the mixture.The electric intensity potential in free water (Φ f T ) and the EDL (Φ d T ) are determined using Gauss's law separately.The obtained electric potential equation for free water is, While the electric potential equation for the EDL can be written as follows, where ζ 0 and ζ r are the relative dielectricity coefficient and the dielectricity coefficient of vacuum, respectively.The finite element formation for the set of governing equations of the multi-ionic transport is described in Appendix A.

Diffusion coefficient of ions
The diffusion coefficient of ions in a porous material depends on pore structure [26,27] and pore solution properties such as the composition of the solution, ionic strength [24,25], and temperature [23].In the present approach, the effective diffusion coefficient is used after applying several corrections to the diffusion coefficient of an ion at infinite dilution.The diffusion coefficient at infinite dilution of individual ions at a reference temperature (T 0 = 298.15K) is adjusted by the temperature, composition, and ionic strength of the solution and tortuosity factor to determine the effective diffusion coefficient.The tortuosity factor depends on the pore structure of the material.The effective diffusivity D i, T eff of the i th ion for the multi-component ion transport model is given as where D i, T 0 ∞ is the diffusion coefficient at infinite dilution of the i th ion at a reference temperature (T 0 ), f τ is the tortuosity factor, and f Tem and f I s are correction factors for temperature and ionic strength, respectively.
The effect of the solution composition (1 + ∂ ln γ i /∂ ln C li ) on the diffusion coefficient of ions is implicitly accounted for in the ion mass transport equation.The diffusion coefficient at infinite dilution is obtained from the limiting conductivities of ions (λ i, T ∞ ) by the following Eq. [44], Due to the presence of a charge-coupling mechanism between the ions, the diffusion coefficient at infinite dilution of many of the aqueous species present in the pore solution of cementitious materials is not determined experimentally.Various models for determining the diffusion coefficient at infinite dilution have been proposed over the past several decades.However, none of the models has led to a comprehensive set of equations suitable for estimating diffusion coefficient at infinite dilution for cementitious ion species [45][46][47].Therefore, in this study, the values of unknown D i, T0 ∞ of ions are estimated empirically as a function of both the size and charge property of the ion.According to Nernst, Stokes, and Einstein [44,46], the ionic stokes radius can be expressed in terms of the limiting conductivities where N A is the Avogadro's number and η w the viscosity of water.From Eqs. ( 12) and ( 13), the ionic diffusion coefficient can be derived as inversely proportional to the formal charge and stokes radius of ions, where: To determine an empirical function for the diffusion coefficient at infinite dilution of ions present in cement-based materials, experimentally reported diffusion coefficient at infinite dilution at a reference temperature of 298.15 K, D i, T=298 ∞ , and the van der Waals radius of ions (calculated using the Calculator Plugins module of the chemical structure visualization tool MarvinSketch, Version 21.10, ChemAxon (htt ps://chemaxon.com/), is used.It was found that the experimentally measured value of D i, T=298 ∞ is strongly correlated to the approximate van der Waals radius and charge of the ion, as shown in Fig. 2.
From Fig. 2, the diffusion coefficient at infinite dilution, D i, T=298 ∞ can be expressed as a function of 1/(|z i | r st i ) for neutral, cations, and anions as The estimated diffusion coefficient of typical ions in the pore solution of a cementitious material at infinite dilution and for a temperature of 298.15K is presented in Table B.1.The molar conductivities of ions change with the ionic strength of the pore solution at a constant temperature.To account for the influence of the ionic strength of solution on the diffusion coefficient at infinite dilution, the model proposed by Appelo [48] is used where: where a 1i and a 2i are coefficients, A the Debye-Hückel parameter (0.51 (mol/dm 3 ) − 0.5 at 25 • C), κ the Debye length, and I s the ionic strength.The limiting ionic conductivities and, equivalently, the diffusion coefficient at infinite dilution are dependent on temperature.For a wide range of temperatures, the temperature dependence of limiting ionic conductivities is reported by Smolyakov's equation [23] ) exp where η T 0 is the viscosity of pure water (Pa⋅s) at T in Kelvin and d i a coefficient.The tortuosity is calculated based on changes of the pore structure due to dissolution and/or precipitation of solid phases determined by chemical equilibrium as follows [13], where f τ is the tortuosity factor at time step t, f τ, 0 the initial tortuosity factor, c a shape factor, and P t a penalty factor, which is described as the ratio of the total porosity at time step t, φ t , and initial porosity, φ 0 .

Chemical equilibrium
In the chemical modelling, the mass exchange terms q f lj and q d lj from the governing Eqs. ( 1) and ( 2) for chemical interactions between ions and solid phases are solved by the external geochemical PHREEQC code [18,19].The PHREEQC module is commonly used as a general-purpose geochemical reaction model to simulate interactions between water, minerals, gases, ion exchangers, surface complexation, and solid solutions [10,12,14,15].However, in reactive transport modelling, transferring cell solution data between multi-ionic transport calculations and chemical equilibrium calculations is tedious and inefficient for running millions of calculations.Therefore, for the compatibility between the transport module and the chemical module with sufficient fast parallel processing in all time steps, the PHREEQC_RM module with an operatorsplitting approach for the reactive transport model is used [49].In both modules (i.e., PHREEQC and PHREEQC_RM) the mass action law is applied for aqueous species, pure phase, and solid solution reaction to obtain chemical equilibrium between aqueous species and solid phases.
The stoichiometric coefficient and equilibrium constant for the reaction of aqueous species, phases, and solid solution phases are defined through a thermodynamic database.

Surface complexation model
It is widely accepted that the pore surface in a cement-based material binds a part of ions present in the pore solution physically see, e.g., [28][29][30][31][32][33][34][35][36].The physical binding of ions affects the pH value of the pore solution by absorbing alkali ions (Na + and K + ) as well as the mass transportation of Cl − and various other ions.Surface complexation allows for accounting for two mechanisms to describe the physical binding of ions.Absorption of ions to a pore surface by an exchange reaction between ionic species in pore solution and the silanol sites (≡SiOH) of the C-S-H, and non-adsorbed ions captured in the EDL present at the solid-solution interface.In the present modelling approach, the total moles of the surface site are calculated by assigning a value to the surface site density, the specific surface area, and the C-S-H product in the phase assemblage.In the pore surface, a charge imbalance will be developed due to the reaction of ion species in the pore solution and the silanol sites (≡SiOH) of the C-S-H, which results in the formation of two distinct layers in the pore solution, i.e., free water and the EDL (see Fig. 1b).In the present study, the charge density on the surface and surface potential are related using the Dzombak and Morel approach [50].The Donnan approach is used to determine the diffuse-layer composition by balancing the surface charge accumulation and computed by the Boltzmann equation.

Numerical examples
The numerical study compares experimental conditions and results from [51] to model predictions to test and validate the developed reactive transport model.The experimental study compares, among others, Cl profile, total elemental contents, i.e., sodium (Na), potassium (K), magnesium (Mg), and sulfur (S), as well as solid phases (Portlandite and carbonate phase) over time for samples exposed to various solutions.A brief overview of the investigated materials and measurements is provided in the following, while more detailed information on sample preparation, exposure, and measurement campaign can be found in [51].

Materials
Mortars samples were prepared in 125 ml cylindrical sealed HDPE bottles using ordinary Portland cement (OPC) type CEM I 42.5 R according to EN197-1 with 6 % silica fume (SF).The PC mortar was proportioned with a water-to-binder mass ratio of 0.40 and a sand-tobinder mass ratio of 2.5:1.The mix design details of the mortar are given in Table 1, while the chemical composition of the OPC, SF, and sand is given in Table 2.

Exposure conditions
Two different exposure solutions were used in the experimental study [51], i.e., NaCl solution and seawater.The seawater sample was collected from the Trondheim fjord, Norway, whereas the NaCl solution was prepared to have the same chloride concentration as the seawater.The measured element composition of the exposure solution is given in Table 3.After 21, 90, and 180 days of exposure, samples were investigated to determine the elemental mass content of chloride (Cl), potassium (K), sodium (Na), magnesium (Mg), and sulfur (S) relative to the

Input parameters and databases for numerical studies
Experimental samples were assumed to be fully saturated during the numerical studies, and the temperature was kept constant throughout the exposure period.Input for the hydration model was determined based on the hydration of each clinker phase composition in the binder materials.The content of main clinkers phases alite, belite, aluminates, and ferrites was determined by XRD-Rietveld [51].The distribution of alkalis content between sulfates and oxides in OPC was calculated based on the measured data of total and soluble alkalis presented in Table 2.The content of calcite is calculated by CO 2 content in OPC.The sulfate content as a substituent in the major phases, especially alite and belite in OPC, is determined using the proposed method of Taylor [52].It is assumed that the remaining sulfates other than alkali-sulfates and sulfates in the major clinkers are accumulated as gypsum phases.The resulting clinker composition of the OPC cement is shown in Table 4.
The ionic composition of the boundary solution (i.e., either NaCl or seawater) was determined using PHREEQC based on the measured elemental composition presented in Table 3, which is presented in Table B.1.In total, 74 ions were included in the multi-ionic mass transport model during the numerical investigations.The CEMDATA18 database [53] (http://www.empa.ch/cemdata) is used to describe reactions and thermodynamic properties for various substances and phases found in cementitious systems.In a CaO − systems, C-S-H gel is formed as the major hydrate described by the solid solution CSHQss model proposed by Kulik [54] in the present study, whereas two end members for KSiOH and NaSiOH are provisionally added to improve predictions of pH and composition of the OPC pore water [53].Moreover, Portlandite (CH), Ettringite (AFtss), monocarbonate (CO 3 -afm), hydrogarnet (HGss), hydrotalcite (Ht), and calcite (CC) are included in the thermodynamic database for numerical investigations.In addition, Kuzel's salt and Friedel's salt phases, primarily for chemical binding of Cl − , magnesium silicate hydrate (MSHss), ferrihydrite-mc (FeOOH), brucite, gibbsite (Gb), gypsum (Gp), straetlingite, and natrolite phases are added in the database to simulate exposure to NaCl and seawater solution.An overview of the dissolution reactions and corresponding equilibrium constants used in the numerical study are shown in Table C.1.
The selected ideal solid solution, CSHQss, contains a mix of CSHQ-TobH, CSHQ-TobD, CSHQ-JenH, CSHQ-JenD, NaSiOH, and KSiOH as endmembers, which have been found as stable phases and in good agreement with measured and calculated alkali concentrations in chemical equilibrium calculations for OPC cement [53].The surface site of the crystal structure of the C-S-H phase has not been entirely determined.Various models can be found in the literature, including Viallis-Terrisse et al.'s model [31], considering one kind of surface site called silanol (≡SiOH).Pointeau et al. [29] mentioned that C-S-H consists of two different sites, i.e., silanol (≡SiOH) and silandiol (=Si(OH) 2 ) sites, while Heath et al. [55] considered that silicon (≡SiOH) and calcium (− CaOH) sites dominate the surface of C-S-H.Based on the availability of surface site data, in the present study, the proposed surface site model by Viallis-Terrisse et al. [31] is used in the surface complexation model.The density of silanol sites of the C-S-H phase is assumed to be equal for all CSHQss endmembers, with a value of 4.878 sites/nm 2 .
For absorption reactions between ions in the pore solution and surface sites, Ca 2+ , (SO 4 ) 2− , (OH) − , Cl − , and alkali ions Na + and K + are included in the surface complexation model.Reaction equations and thermodynamic parameters, log(K), for ion adsorption on the C-S-H surface are shown in Table C.2 [21].The equilibrium constants for ion reactions onto silanol sites of CSHQss endmembers are assumed to be equal.In addition, the specific surface area per mole of C-S-H, A ss , and the thickness of the EDL, t d , are required for the surface complexation model.The specific surface area, A ss , is calculated by the following

equation,
n CSH (19) where A int is the total surface area of the interlayer, dA gl r the surface area of a gel pore for a pore radius, dA cp r the surface area of a capillary pore for a pore radius r, f CSH Ain and f CSH Agl_cp the fraction of the C-S-H surface in the interlayer, gel, and capillary pore surface, respectively, and n CSH the total moles of the C-S-H product.The pore structure model determines the surface area of the interlayer, gel, and capillary porosity.It is thereby assumed that the C-S-H phase occupies 65 % of the gel and capillary pore surface, and C-S-H forms the interlayer.Therefore, the values for f CSH Ain and f CSH Agl_cp are set to 1.0 and 0.65, respectively.The thickness of the EDL, t s , depends on the ionic strength of the pore solution and is set to be proportional to the Debye length, κ, (t d = χ κ, where χ is the proportionality coefficient) in the present study.The Debye length (for electrolyte solution) is determined as follows where, k B is the Boltzmann constant and e the elementary charge.In this study, the proportionality coefficient, χ, is set to 1.55 after [14].Details of spatial discretization and time step used in the numerical study, together with assumed threshold pore radius, initial tortuosity, shape factor, and other constant values, are shown in Table 5.In general, using a relatively high number of elements to discretize the domain in numerical studies will increase the accuracy; However, computational costs will increase.Therefore, to balance accuracy and computational cost in the present investigations, the 30 mm 1D domain is discretized using a total 50 of elements with a growth factor to create a finer mesh near the surface of the exposed boundary.The time step is decided through sensitivity studies for ion mass transport based on comparing total exposure time, computational costs, developed charged imbalance in the pore solution, and the truncated error in the boundary node [56].Chemical equilibrium is determined every second time step, which is sufficient to reach a minimum level of accumulated ion concentration in the pore solution.

Phase assemblage of reference mortar samples
The numerical results for the hydration of major clinker phases, the average level of hydration of cement, and the phase assemblage for the reference mortar sample after three years of sealed curing at 20 • C are shown in Fig. 3. Results of the hydration model indicate that >90 % of alite and aluminate clinker are hydrated after three years of sealed curing and that the dissolution rates for belite and ferrite are significantly slower in the reference mortar specimen.The predicted phase assemblage contains C-S-H, Portlandite, Ettringite, monocarbonate, hydrogarnet, hydrotalcite, calcium carbonate, and unhydrated cement.The simulated phase assemblage for OPC cement is thereby in good agreement with reported results in, e.g.[53,57].Nearly 31 % of the total mortar volume consists thereby of both hydrated and unhydrated cement, whereas the numerical results indicate that approximately 12 % of the total volume is pore space after three years of sealed curing.Due to the equilibrium between pore solution and hydrated phase of a standard Portland cement, pore solution is formed with high alkalinity, which maintains the pH value of the pore solution generally between 12.5 and 14 [58].The simulated results show that the pH value of the pore solution of the reference mortar is approximately 13.22.
Dissolved ions in the pore solution will be transported among others by electrochemical potential forces, and these transported ions will subsequently affect the dissolution and precipitation of phases.Thus, validating the pore solution composition against experiment data is crucial for RMT modelling.More than 70 ions are considered in the present study, for which experimental determination is cumbersome.Hence, the elemental concentration of selected ions in the pore solution is compared with the available experimental results for the reference mortar sample in Fig. 4. Presented results indicate that numerical results are in good agreement with experimental results obtained through ICP-MS [51], apart from Fe and Mg, for which the model predicts lower concentrations in the pore solution.The discrepancy may be due to the influence of other factors such as the concentration of other elements, ionic strength, etc., on the dissolution rate of Fe and Mg containing solid phases, i.e., hydrogarnet and hydrotalcite.Therefore, the choice of phases and the dissolution rate of phases within the thermodynamic modelling framework strongly affect and control the composition and concentration of the pore solution.In the present study, incorporating a hydration model and the selected phase assemblage (see Fig. 3) allows capturing the trend for most of the elements in the pore solution.

Mortar samples exposed to NaCl and seawater solution 4.2.1. Solid phase composition
The presented modelling framework was used to determine the effect of NaCl and seawater on the solid phases of exposed mortar samples.Results of the numerical simulations, i.e., phase assemblage for OPC mortar specimens immersed in NaCl and seawater solution for 180 days, are shown in Fig. 5.The right-hand side of the figures displays the paste composition in the unaffected core, while the left-hand side illustrates the effect of exposure on the mortar samples.For the non-exposed core, the predicted phase composition includes CSH, Portlandite, Ettringite, monocarbonate, hydrogarnet, hydrotalcite, calcium carbonate, and unhydrated cement.
Upon exposure of the mortar sample to NaCl solution, Friedel's salt and ettringite (Aftss) are formed in the vicinity of the exposed surface, whereas Portlandite and monocarbonate decompose with a continuous increase in calcite content (see Fig. 5a).The observed changes in the phase assemblage result from the presence of Cl − ions from the boundary solution, which start to substitute carbonate ions (CO 3 ) 2− in the monocarbonate to form Friedel's salt, see, e.g., Balonis et al. [59].The released carbonate ions subsequently react with Ca 2+ ion in the pore solution to form calcite.The Ca 2+ deficit in the pore solution causes an increase in the dissolution rate of Portlandite and monocarbonate until both phases are decomposed near the exposed surface.Moreover, leaching of selected ions is observed from the exposed mortar specimen as the concentration of most ions (considered in the RTM model) is close to zero in the boundary solution.In addition, a minute decrease in the volume of solid phases is observed near the exposed surface due to the leaching of ions.
In the case of seawater exposure, the formation of Friedel's salt as observed for exposure to NaCl was not predicted (see Fig. 5b).On the other hand, the formation of additional phases such as ettringite, ferrihydrite-mc, brucite, and gypsum is predicted near the surface of the exposed mortar specimen due to the presence of sulfate, carbonate, and magnesium in the exposure solution.A similar trend was noticed for the mortar specimen exposed to NaCl solution in connection to monocarbonate and Portlandite.Both phases continuously decompose with exposure time until fully dissolved near the exposed surface.The decomposition of monocarbonate and Portlandite is thereby accompanied by an increase in calcite content as observed for the mortar specimen exposed to NaCl solution.In addition to the decomposition of monocarbonate and Portlandite, the decalcification of C-S-H and the dissolution hydrogarnet are observed near the exposed surface.The porosity near the exposed surface decreases mainly with the ettringite formation because the molar volume of ettringite is comparably higher.For both exposure condition and except the region near the exposed surface (0-2 mm section), the gradual decreases in pH value was observed together with the leaching of some of ions and decomposition of phases such as monocarbonate and Portlandite, whereas the decalcification of higher Ca/Si ratio C-S-H phase (CSHQ-JenD) due to seawater exposure causes a sudden drop in pH value at the vicinity of the exposure surface (approximately 11.83 at the exposed surface).Fig. 6 illustrates a comparison between predicted and experimentally determined (through TGA) Portlandite and carbonate profiles for specimens exposed to NaCl and seawater, respectively, after various times of exposure.The Portlandite and carbonate contents are presented in total phase content as the mass percentage of the dried mortar at 900 • C. Results illustrate an increase in the carbonate content in the outermost section (0-5 mm) for both exposure solutions over time (see Fig. 6b).The formation of carbonate phases may result from either carbonation during sample preparation or transport of CO 2 dissolved in the boundary solution inside the mortar sample during exposure.In the present study, the dissolution of CO 2 in the exposure solution was considered under atmospheric conditions for which a good agreement between experimental results and numerical predictions is observed.Moreover, a reduction in Portlandite content near the exposed surface is observed for both exposure conditions over time.As previously discussed, the dissolution of Portlandite is increased by calcite formation.A gradual decrease in Portlandite is also found in the deeper, uncarbonated sections, i.e., both experimentally and numerically (see Fig. 6a).The dissolution of Portlandite in the uncarbonated area results thereby from the leaching of (OH) − ions, which is accompanied by the transport of Ca 2+ ions inside the mortar.

Pore structure
Fig. 7a and Fig. 7b show the simulated cumulative pore volume and pore size distribution of capillary, gel, and total pores of OPC mortar at a depth of 0.1 mm from the exposed surface before and after exposure to NaCl and seawater solution for 180 days.Numerical results for the cumulative pore volume indicate seawater exposure appear to fill pores with newly formed phases, whereas NaCl exposure shows leaching behaviour.For both exposures scenarios, a slight decrease in the total pore volume of gel pores is found because decomposition of the C-S-H phase occurs due to exposure (see Fig. 5).The effect is more pronounced for seawater exposure, for which nearly half of the capillary pore space is filled by newly formed phases, while an increase in total porosity of capillary pores is observed due to NaCl exposure.The peak radius has a considerable effect on the ionic transport in the RTM model, influencing ion flux through a change in threshold porosity.The pore size distribution shows that the peak lies between 10 nm and 1 μm for all non-exposed and exposed mortar specimens.A more significant reduction in peak radius is observed for seawater exposure than NaCl exposure.For seawater exposure, the formation of additional phases affects the total porosity distribution and reduces the peak radius ten times smaller than non-exposed OPC mortar, whereas a considerable change is not found for NaCl exposure.
Numerical predictions of the tortuosity factor for the OPC mortar specimens after exposure to NaCl and seawater solution for 180 days are given in Fig. 8. Results of the numerical simulations illustrate changes in the tortuosity factor near the exposed surface for both NaCl and seawater exposure.For seawater exposure, the tortuosity factor is significantly more reduced than NaCl exposure, indicating an increase in the length of the ion pathway.Consequently, the ion flux through the pore space near the mortar's surface is reduced, forming new solid phases, although the boundary solution provides sufficient ion supply.A peak reduction in the tortuosity factor was observed at a certain depth from the exposed surface, marking the transition between decrease and increase in the total volume of solid phases.

Chloride profiles
The predicted and experimentally measured distribution of chlorides in mortar exposed to NaCl and seawater solution after 21, 90, and 180 days are shown in Fig. 9.The chloride profiles are thereby presented as total chloride contents per mass percentage of the dried mortar (at 105 • C).Within the presented modelling framework, chlorides may be stored in free water and the EDL of the pore solution, physically bound on the pore surface, and chemically bound in, e.g., Kuzel's salt and Friedel's salt.The sum of all chlorides, i.e., free and physically and chemically bound, is then provided as the total chlorides in the predicted result.Generally, the presented results indicate an excellent agreement between simulation and the experimental results for most parts of the spatial and time domains for all exposed mortar specimens.Deviations between experimental and numerical results are observed for the surface region (0-2 mm section).These deviations may be caused by troublesome sample preparation from the mortar surface for chloride content measurements as well as the casting of mortar, including variations in the mix proportion at the mortar surface due to preparation and bleeding.For all exposure times, the predicted and measured total Cl results indicate that seawater exposure seems to lead to a higher chloride content than exposure to NaCl solution at the same depth from the exposed surface (see Fig. 9b).
On the other hand, the results show the peak behaviour in chloride profile near the exposed surface (0-3 mm section) of the mortar due to the NaCl exposure because the Cl from NaCl boundary solution was immediately stored as Friedel's salt near the exposed surface (see Fig. 9a).In addition, over time, the predicted and measured chloride [NaCl] [Pore structure] [Seawater] profiles for both exposure conditions show a reduction in the chloride content at the surface (0-1 mm section).The decrease in chloride content near the exposed surface is caused by the release of initially physically absorbed Cl − ions from the surface site of C-S-H due to the decalcification of higher Ca/Si ratio C-S-H phase (CSHQ-JenD), whereas the mass of dried mortar at 105 • C increased in the surface section.Therefore, the normalized total chloride contents with the mass of the dried mortar at 105 • C reduced after some exposure period.

Elemental profiles
Results of the numerical simulations comprising profiles of total elements, the composition of free water, composition of the EDL, absorbed element on the surface sites, and composition of surface sites for the mortar specimens exposed to NaCl and seawater solution for 180 days are given in Fig. 10.The quantity for all compositions is expressed as moles of an element per m 3 volume of mortar.As illustrated in Fig. 10a, changes in the elemental composition are observed mainly near the exposed surface, while the ingress of Cl and leaching of K are observed up to a depth of approximately 15 mm.The ingress of Mg and S is identified only for seawater exposure as the exposure solution contains certain amounts of Mg and S, whereas a unique shape in Ca and Na profiles were found from the simulation results.
The free water and EDL compositions are expressed as the content of an element instead of the ionic concentration in Fig. 10b and Fig. 10c, respectively.For example, Ca concentration in solution indicates the sum of the concentration of Ca 2+ in different species such as Ca 2+ , Ca (OH) + , etc.For NaCl exposure, the ingress of Na and Cl ions through free water is observed; simultaneously, Na ions in the EDL are moving outwards.Moreover, leaching is observed for K and S, while the concentration of Ca and Mg increases near the exposed surface, despite a lack of supply from the exposure solution.Upon mortar sample exposure to seawater solution, Na and Cl profile tendency is the same as for NaCl exposure.The leaching of K is lower compared to samples exposed to NaCl due to the presence of K in the seawater solution.S, Ca, and Mg profiles indicate ingress inside the mortar due to higher concentrations in the exposure solution than the pore solution.
Fig. 10d and Fig. 10e show the amount of absorbed elements onto surface sites and the composition of surface sites in the C-S-H phase after 180 days of exposure to NaCl and seawater solution.The right-hand side of the figures displays unaffected surface sites, while the left-hand side of the figures shows the effect of transported ions on the silanol sites.Initially absorbed alkali ions, i.e., for both NaCl and seawater exposure, Na and K, onto the surface sites are released to the pore solution due to leaching of ions and decalcification of C-S-H phase near exposure surface (see Fig. 5).On the other hand, the absorption reaction of Ca and Cl increased on the surface sites compared to the initial stage of absorption.A peak in the absorbed element profile of Ca and Cl is observed at a certain depth, caused by a release of absorbed Ca and Cl due to the decalcification of the C-S-H phase near the exposure surface.In the case of NaCl exposure, some of the absorbed S is released and leached through the exposure surface.In contrast, a certain amount of transported S from the boundary solution is absorbed onto the surface sites for seawater exposure.
A comparison between numerical and experimental results for sodium (Na), potassium (K), magnesium (Mg), and sulfur (S) after 21, 90, and 180 days of exposure to NaCl and seawater is given in Fig. 11.Total elemental distributions are presented as mass percentage of the dried mortar (at 105 • C).In general, a good agreement between numerically predicted and experimentally determined elemental distribution is found for most parts of the spatial as well as time domain for all exposed mortar specimens.
Both exposure solutions, i.e., seawater and NaCl solution, contain >3-times the Na concentration than the pore solution of the mortars (see Table 3 and Fig. 4); therefore, Na is expected to ingress into the samples due to electrochemical and potential gradients.However, except for the outermost section (0-1 mm) of the mortar specimen exposed to NaCl, the measured and predicted sodium profiles indicate enrichment in sodium only in deeper sections from the exposed surface (3-15 mm).At the same time, the results illustrate a reduction in sodium content in the section located 1-3 mm from the exposed surface.The presented modelling framework provides an explanation for the observed (experimentally and numerically) sodium distribution.Initially, transported sodium ions are absorbed on the pore surface through the reaction between ions and surface sites of the C-S-H.Later, the absorbed sodium ions are released from the pore surface due to the lowering pH value, and the released sodium ions are subsequently transported inside the mortar due to the concentration gradient between the boundary and pore solution.Thus, sodium enrichment is observed in the deeper section.
The exposure solutions contain only limited potassium concentrations (in the case of seawater and none for NaCl) compared to the pore solution (see Table 3 and Fig. 4), leading to continuous leaching of potassium ions from the exposed surface.In addition, initially physically absorbed potassium ions from the pore surface are released to the pore solution (see Fig. 10a).The modelling framework captures the effect of leaching very well on the profile of potassium over time for both exposures.
The presence of magnesium in the seawater solution leads to the ingress of magnesium in the mortar samples.However, the magnesium content increased only in the outer sections (0-2 mm) of samples  exposed to seawater.Numerical results predict the formation of hydrotalcite and brucite in the outer sections of the exposed samples (see Fig. 5) as the transported magnesium ions react with other ions in the pore solution and precipitate as solid phases.For samples exposed to NaCl solution, leaching of magnesium is expected due to the lack of Mg in the exposure solution.However, the leaching is restricted due to a slight concentration gradient between boundary and pore solution as well as the dissolution of the initially formed hydrotalcite (see Fig. 5 and Fig. 10a), resulting in a temporal increase of magnesium in the pore solution.Thus, the magnesium concentration in samples exposed to NaCl remains unchanged, as indicated by experimental results and numerical predictions.
Sulfur is not present in the NaCl solution, whereas the seawater solution contains slightly more sulfur (see Table 3 and Fig. 4) than the pore solution of the unexposed mortar.Due to the presence of sulfur in the seawater solution, experimental and numerical results show an increased sulfur content in the outer section (0-4 mm) of the exposed sample.As a result, additional ettringite (solid solution) is formed, as indicated in Fig. 5.For samples exposed to NaCl, sulfur leaching occurred due to concentration gradients in the pore solution.However, only a minor reduction in the sulfur content is observed in the outer section (0-4 mm) of the exposed mortar.

Conclusions
A reactive transport model (RTM) was developed for saturated cement-based materials based on multi-ionic transport through an uncharged, 'free water' and charged, 'EDL' coupled with a chemical equilibrium model for interaction between the pore solution and solid phase.A set of modified Poisson-Nernst-Planck (PNP) equations, including the effect of temperature, pore solution property, and pore structure changes, were introduced for multi-ionic transport through free water and the EDL.The diffusion coefficient at infinite dilution of many of the aqueous species present in the pore solution in cementitious material was expressed as a function of van der Waals radius and charge of the ion.
Numerical predictions were compared to experimental data to demonstrate the potential of the RTM.Comparisons include the variation of total element content, solid-phase composition, free water and EDL composition, absorbed elements onto surface sites, and composition of surface sites of C-S-H phase for mortar specimens exposed to NaCl solution and seawater from Trondheim fjord in Norway.Excellent agreements (in both spatial and time domains) between experimentally determined and numerical predictions were found comparing Portlandite and carbonate in the phase assemblage, total chloride (Cl), elemental distribution of sodium (Na), potassium (K), magnesium (Mg), and sulfur (S) for all exposed specimens adjusting only the initial tortuosity factor.Furthermore, the developed RTM provides an explanation for the observed elemental distribution, phase assemblage, and pore structure in all exposed mortar specimens.The results indicate the importance of considering multi-ionic transport through free water and the EDL for cement-based material due to exposure to the environment.Moreover, considering the effect of temperature, pore solution properties, and pore structure changes on multi-ionic transport in porous media significantly improve model prediction.

CRediT authorship contribution statement
All persons who meet authorship criteria are listed as authors, and all authors certify that they have participated sufficiently in the work to take public responsibility for the content, including participation in the concept, design, analysis, writing, or revision of the manuscript.Furthermore, each author certifies that this manuscript or similar content in manuscript has not been and will not be submitted to or published in any other publication before its appearance in the Cement and Concrete Research.
The specific contributions made by each author are listed below.

Appendix A. Finite element formulation
The present approach solves the governing equations through the finite element method (FEM) for the multi-ionic transport.For the finite element formulation, the governing Eqs. ( 6), ( 7), (9), and ( 10) are transformed into weak forms dividing the system of equations in a spatial domain and a time domain using a weight function w for the spatial domain and a weight function W for the time domain.Galerkin's approximation is used for the discretization of the one-dimensional spatial scheme using linear spatial elements [60,61].The Green-Gauss theorem is used in the ion mass transport terms to separate the boundary flow conditions.An implicit single-step approach is used to solve the problem in the time domain.The complete coupled set of equations is solved simultaneously to avoid staggering schemes between separate solutions of state variables.A modified Newton-Raphson iteration scheme is used to improve the FEM solution of the non-linear governing equations for the employed single time-stepping scheme.The modified iteration scheme reduces computational time compared to the complete scheme while minimizing the residual of the problem [9,62].The weak forms of the local differential equations, i.e., Eq. ( 6), for the concentration of the ionic species in the free water solution of the porous material, are where a time increment is used in the integration, starting with t1 and ending with t2, J f lj the mass density flow of ions at the boundary surface dS of free water, and n the outward drawn normal to this surface.The weak forms of the local differential equations, i.e., Eq. ( 7) for the concentration of the ionic species in the EDL of the porous material, are where C represents the total damping of the complete coupled set of equations and K and f, are the total stiffness of the problem and the load vector, respectively.Δt is time increment, and Θ is the time integration parameter.The time integration parameter Θ is restricted to 0 ≤ Θ ≤ 1, in which Θ = 0 represents a truly explicit scheme, and Θ = 1 is a truly implicit scheme.The total damping matrix of the total coupled problem (A.7) can be written as  In order to shorten the notation, the abbreviation â = a n + Θ(a n+1 − a n ) is introduced.The stiff-ness part of the total coupled problem described in Eq. (A.7) can be formed as S. Sharmilan et al.
The submatrices of K in Eq. (A.10) are identified as (using the weak Eqs.(A.1)-(A.4)) Ion Transport through Free water and EDL water Update Boundary condition -Boundary solution Con for BC_Free water -EDL Con for BC_EDL -Temperature of boundary

Fig. 1 .
Fig. 1.Schematic diagram of the proposed multi-species reactive transport model.

Fig. 3 .Fig. 4 .Fig. 5 .
Fig. 3. Numerical results for reference OPC mortar sample.(a) Hydration of major clinkers phases and average degree of hydration over three years, (b) phase assemblage for mortar specimens after three years of sealed curing at 20 • C.

Fig. 6 .
Fig.6.Comparison of experimentally determined[51] and numerical predictions of Portlandite and carbonate content after 21, 90, and 180 days of exposure to NaCl and seawater solution, Phase contents were measured through TGA on profile ground OPC mortar and expressed as wt% of mortar dried at 900 • C (left: NaCl, right: seawater).

Fig. 7 .
Fig. 7. Cumulative pore volume and pore size distribution for OPC mortar specimens at a depth of 0.1 mm from the exposed surface after 180 days (left: NaCl, right: seawater).

Fig. 8 .
Fig. 8. Tortuosity factor of OPC mortar specimens exposed to NaCl and seawater solution for various exposure times (left: NaCl, right: seawater).

S
.Sharmilan et al.

(A. 6 )
T A lj v j (A s − A th )t d lj is the mass density flow of ions at the boundary surface dS of the EDL.The weak form of the Poisson type of Eq. (9) determining the electrical potential in the free water is given by ∫ Φ is the electrical displacement at the boundary surface of the free water space.The weak form of the Poisson type of Eq. (10) determining the electrical potential in the EDL is given by ∫T (A s − A th )t d ζ 0 ζ r ∇Φ T d dV dt = − Φ isthe electrical displacement at the boundary surface of the EDL.The state variables in the weak formulations are approximated by the general expansion N • a, where N is the global shape function, and a contains the nodal state variables.In the one dimensional case, the local linear shape function is given as [1 − x/l e , x/l e ], where l e is the element length.The arbitrary spatial weight function w is approximated with the same general expansion following Galerkin's method.The state variables are approximated by the shape functions N, i.e.w = c T N T ; c T is an arbitrary matrix of the spatial weight function, and yields c T = c.The weak formulations are formed using the gradient of the shape functions N, the approximate gradient of the state variables is denoted by B. The gradient of the weight function and the state variables may be written as follows: ∇w = c T B T ; ∇C lj f = Ba j f ; ∇Φ T f = Ba Φ f ; ∇C lj d = Ba j d ; ∇Φ T d = Ba Φ d Applying assumptions (A.5) and (A.6), the weak formulations (see Eqs. (A.1)-(A.4))can be brought to the form C(a n+1 − a n ) Δt + K(a n + Θ(a n+1 − a n ) ) − (f n + Θ(f n+1 − f n ) ) = 0 (A.7)

8 )
According to the weak form Eqs. (A.1)-(A.4),sub-matrices of C in Eq. (A.8) are identified, asC i f = ∫ V N T φ gc ε l f NdV; C j d = ∫ V N T A s t d N dV (A.9)

Table 5
Model parameters for the numerical investigations.