Model for alkali-silica reaction expansions in concrete using zero-thickness chemo-mechanical interface elements

Alkali-Silica Reaction (ASR) is a particular type of chemical reaction in concrete, which produces cracking and overall expansion of the affected structural element due to the formation of expansive reaction products within the cracks. This paper develops the formulation of a coupled Chemo-Mechanical (C-M) Finite Element (FE) model for simulating ASR expansions in soda-lime glass concrete at the meso-scale. The model considers several C-M coupling mechanisms, including a reaction-expansion mechanism qualita- tively proposed by the authors elsewhere on the basis of experimental results, which is introduced in order to reproduce the effect of compressive stresses on the development of ASR expansions. The model has the characteristic ingredient of using zero-thickness interface elements for modelling the C-M mechanisms leading to the propagation of cracks due to the formation of ASR products within them. This fact has required the development of: (i) a new FE formulation for diffusion-reaction processes occurring within discontinuities represented by interface elements, and (ii) a new mechanical constitutive law for interface elements, which is able to reproduce the propagation of a crack induced by the development of an internal pressure exerted by solid reaction products formed within it. In addition, the numerical implementation of the diffusion-reaction formulation has been advantageously performed with clear separation between the boundary-value or ‘structural’ governing equations (i.e. continuity and concen- tration gradient relations), and the ‘constitutive’ (i.e. chemical) equations. The model is illustrated with some application examples. (cid:1) 2020 The Author(s). Published by Elsevier

Model for alkali-silica reaction expansions in concrete using zero-thickness chemo-mechanical interface elements

Introduction
Alkali-Silica Reaction (ASR) is a particular type of chemical reaction in concrete involving the alkaline pore solution of the Hydrated Cement Paste (HCP) and various metastable forms of silica present in many natural and synthetic aggregate particles used in concrete (Poole, 1992;Hobbs, 1988). The silica structure is dissolved by the attack of hydroxyl anions, passing to the pore solution in aqueous forms which later recombine with calcium and alkali cations (always present in the concrete pore solution) to form silicate hydrates of variable stoichiometric composition, the so-called 'ASR gel'. This reaction product is usually hygroscopic and swells as it absorbs water. The degree of hygroscopy and the associated volume change of the ASR gel depends mainly on its chemical composition and the relative humidity in the concrete pores. If the space available at the reaction site, i.e. pores, cracks, etc., is not enough to allocate the ASR gel, the 'swelling pressures produced by the gel induce the formation of microcracks close to the reaction sites, and these propagate and coalesce to produce cracking within the fabric of the concrete and overall expansion of the structural element affected' (Poole, 1992).
The ASR expansions in concrete are very much affected by external factors, such as temperature, moisture and external mechanical loads. Of these three factors, the latter is by far the one which has received the least attention, in spite of being of direct interest for practical engineers dealing with the prognosis of structures affected by ASR. There is field and experimental evidence indicating that the concrete stress state affects the magnitude, distribution and kinetics of the ASR expansions, e.g. (Charlwood et al., 1992;Herrador et al., 2008;Multon and Toutlemonde, 2006;Gautam et al., 2017). However, this evidence seems sometimes contradictory since the Chemo-Mechanical (C-M) mechanisms underlying these effects are not yet completely understood.
Along this line, the authors have undertaken a research project aiming at improving the understanding of the mechanisms by which the stress state affects the development of concrete expansions and cracking due to ASR. In the framework of this research project, the authors and co-workers have reported experimental results obtained with concrete specimens made with soda-lime (SL) glass as reactive aggregate Liaudat et al., 2019). Based on these results and on results of other researchers found in the literature, a qualitative meso-scale reaction-expansion mechanism has been proposed which may explain the effects of stress on ASR expansions in concrete (Liaudat et al., 2019). The validation of this mechanism requires to prove its ability to reproduce the observed experimental behaviour. For this purpose, in the present paper, this mechanism is theoretical and numerically formulated, obtaining a coupled C-M Finite Element (FE) model for ASR expansions in SL glass concrete.
In the last three decades, a wide variety of numerical models have been proposed for simulating ASR expansions in concrete, e.g. (Leger et al., 1996;Bažant and Steffens, 2000;Ulm et al., 2000;Suwito et al., 2002;Li and Coussy, 2002;Steffens et al., 2003;Capra and Sellier, 2003;Farage et al., 2004;Poyet et al., 2006a;Poyet et al., 2006b;Fairbairn et al., 2006;Herrador et al., 2009;Multon et al., 2009;Comi et al., 2009;Sellier et al., 2009;Comi and Perego, 2011;Comi et al., 2012;Pignatelli et al., 2013;Puatatsananon and Saouma, 2013;Charpin and Ehrlacher, 2014;Saouma et al., 2015;Morenon et al., 2017;Esposito and Hendriks, 2017;Huang et al., 2018;Iskhakov et al., 2019). For extensive literature reviews on ASR models, the reader may refer to Pan et al. (2012), Esposito and Hendriks (2017) and Saouma (2014). Different scales of analysis have been proposed in those models: micro-, meso-and macro-scale. From our point of view, the study of the C-M mechanisms leading to ASR expansions requires a meso-scale analysis, i.e. to explicitly consider the reactive aggregates embedded in the cement paste. However, to the best of the authors' knowledge, only three models in the literature have been proposed at this scale: the model by Comby-Peyrot et al. (2009), the model by Dunant et al. (Dunant and Scrivener, 2010;Dunant and Scrivener, 2012a;Dunant and Scrivener, 2012b;Giorla et al., 2015), and the model by Alnaggar et al. (2013), Alnaggar et al. (2015) and Alnaggar et al. (2017). An important conceptual limitation shared by these three models is that the proposed C-M coupling is unidirectional. That is, the formation of ASR products provokes internal swelling of the aggregates (Comby-Peyrot et al., 2009;Alnaggar et al., 2013;Alnaggar et al., 2015;Alnaggar et al., 2017) or of 'gel pockets' (Dunant and Scrivener, 2010;Dunant and Scrivener, 2012a;Dunant and Scrivener, 2012b;Giorla et al., 2015) within the aggregates, but the material stress state and/or cracking have no effect on the rate of formation of the ASR products. Among other contributions, the model proposed in this paper overcomes this limitation by explicitly considering the effect of cracking on the diffusion-reaction processes leading to the formation of ASR products, as well as the effect of compression stresses at the reaction sites on the chemical kinetics.
In this model, the diffusion-reaction problem and the mechanical problem are treated separately in two partial FE models linked together by means of a number of coupling terms. These models are implemented in two different codes, which are called sequentially in an iterative procedure. Both codes are used to analyse the same FE mesh, in which the reactive aggregates are explicitly represented and they are embedded in a cementitious matrix phase representing non-reactive mortar or HCP. This mesh includes zero-thickness interface elements which are inserted a priori along all the aggregate-matrix contacts and also along some predefined matrix-matrix and aggregate-aggregate inter-elements boundaries in order to represent the main potential crack paths. In the case of the aggregate-matrix contacts, the interface elements also make it possible to represent the specific properties of the Interfacial Transition Zones (ITZ).
The diffusion-reaction formulation consists of three primary diffusion/reaction fields for the concentrations of aqueous silicate, calcium and alkalis in the pore solution, complemented by a number of chemical kinetics and chemical equilibrium equations. The amounts of solid chemical species (reactive silica, portlandite, and ASR products) evolve with the progress of the reactions, eventually leading to the development of internal pressures. In analogy to traditional formulations of the mechanical problem, the resulting diffusion-reaction formulation is structured in two parts. On one hand, the general balance, compatibility, and discretization equations, developed for a generic a diffusion-reaction problem, and, on the other, the particular 'constitutive' equations developed for the ASR model proposed in this paper. In this way, different constitutive (chemical and diffusivity) laws can be implemented as a library and therefore the same code may be used for analysing different modelling assumptions or different diffusion-reaction problems.
With the purpose of reproducing the C-M mechanisms leading to the propagation of cracks due to the formation of ASR products within them, a new FE formulation for diffusion-reaction processes occurring within discontinuities represented by interface elements is developed in this paper. In addition, an existing mechanical constitutive law for interface elements based on the fracture mechanics theory is modified in order to model the interactions between a crack and the solid reaction products formed within it.
The content of this the paper is organized as follows. After this introduction, in Section 2 the main modelling assumptions are established, regarding concrete micro-structure, reactivetransport processes, and chemo-mechanical couplings. In Section 3, these assumptions are used for developing the formulation and numerical implementation of the FE model for the diffusionreaction problem, with the distinctive feature of using zerothickness interface elements with double nodes representing diffusion-reaction processes in discontinuities. The mechanical problem is treated in the following Section 4. Since its formulation is based on an existing and well documented meso-scale model, the discussion in this section is focused on the development of a mechanical constitutive law for the propagation of cracks due to the formation of ASR products within discontinuities. Once the diffusion-reaction model and the mechanical model have been presented, the staggered procedure used for establishing the coupling between them is discussed in the following Section 5. In Sections 7 and 8, the capabilities of the model are illustrated with some application examples. Finally, the paper concludes in Section 9 with a series of remarks on the most distinctive aspects of the model developed.

General
The model is conceived to be used to simulate the experimental tests such as (Liaudat, 2018;Liaudat et al., 2019). In those tests, two stages have been distinguished: the curing stage, at room temperature for a period of time greater than 28 days, and the exposure stage at 60°C. Since the extent of the ASR developed during the curing stage is much smaller than the one developed during the exposure stage, the simulations will only consider the ASR occurring during the exposure stage. It is also assumed, for the sake of simplicity, that before starting the exposure stage, the Portland cement in the specimens had reached its maximum possible hydration degree, and, therefore, no changes in the HCP microstructure due to clinker hydration during the exposure stage are considered. Moreover, the specimen pores and cracks are assumed to be fully saturated with water at all times during the tests.

Concrete micro-and meso-structure
SL glass concrete is considered as formed by two homogeneous continuum components. On one hand, the SL glass aggregates, and, on the other, the 'cementitious matrix' (HCP or mortar). In the analysis, each of these components is represented by means of continuum finite elements, and in between them zero-thickness interface elements are inserted to represent the particular characteristics of the cement/glass contacts. Interface elements are introduced also in between continuum elements of the same component in order to simulate the propagation of cracks.
The SL glass is considered to be non-porous, impervious, and constituted by a unique volume fraction. The cementitious matrix, in contrast, is considered to have three different volume fractions: capillary pores, portlandite, and Inert Cement Paste (ICP), the latter corresponding to cement hydration products which are considered not to have any chemical interaction with ASR, e.g. C-S-H, AFt, or AFm. The volume of ICP includes both the volume of solid hydration products and the volume of gel water absorbed by them. The volume of gel water is considered to be part of the volume of cement pore solution, or, what is the same, given that the material is assumed to be fully saturated, the volume of gel pores is considered to be part of the total volume of pores. The remaining volume occupied by pore solution is referred as 'capillary pores' even if it may include other voids such as macro-pores or cracks. The diffusion-reaction processes are allowed to occur in the gel pores, but the reaction products are not allowed to precipitate within gel pores. All these volume fractions, as well as the concentration of aqueous species in pore solution, are assumed to be known at the beginning of the exposure stage.

Transport processes
The only transport process considered is the Fickian diffusion of aqueous species due to concentration gradients. These diffusion processes are allowed through the continuum FEs representing the cementitious matrix and along the interface elements representing glass/cement interfaces and open cracks.
Bulk glass is assumed to be non-porous and impervious. Only if the glass is cracked, the pore solution can penetrate the glass particle boundaries in the form of localized fluxes along the open cracks. Otherwise, only the external surface of the glass particle, in contact with the cement paste, is exposed to the cement pore solution. The effective diffusivity and porosity of the cementitious matrix are assumed to remain constant during the simulations.
The dissolution/formation of solid species will alter the local balance of pore water by two (related) mechanisms. On one hand, the dissolution/formation of solid species will release/consume water to/from the pore space (see Section 2.4). On the other, the dissolution of solids will increase the capillary pore fraction, while the formation of solids will reduce it. It is assumed that the water flows associated to these mechanisms occur very rapidly in such a way that the fully saturated condition in the capillary pores can be considered to be restored instantaneously, as well as the condition of negligible water pressure which is basic assumption in this study.

Chemical reactions
It is assumed that ASR in SL glass concrete can be suitably represented by considering the following chemical reactions: iii) Dissolution/formation of high-calcium ASR product (also denoted as RPA) iv) Dissolution/formation of low-calcium ASR product (also denoted as RPB) Note that, in the above equations, sodium (Na) and potassium (K) are treated indifferently as 'alkalis' (R). The effect of the remaining chemical species usually present in concrete on the development of ASR, is neglected.
The SL glass is considered to be composed only by SiO 2 , CaO, and R 2 O. Other minor components are neglected. For this study, a typical SL glass composition of 76.7 wt% SiO 2 , 14.1 wt% CaO, and 9.2 wt% R 2 O is adopted.
It is assumed that, in contact with alkaline solutions, the SiO 2 contained in SL glass dissolves according to the process described by Rajabipour et al. (2015), and that Ca 2+ and R + are congruently released. For the sake of simplicity, the whole dissolution process is reduced to a single reaction with a fixed stoichiometry. Moreover, only one aqueous silica species is considered, a simplification that gives good results for solutions with pH 11 to 12. Out of this range, other aqueous silica species (SiO 4 H 2 2À above, and SiO 4 H 4 0 below this range) should be considered for an accurate estimate of silica solubility. The ASR product is assumed to be a simple mixture of two components of constant composition, as suggested by Helmuth and Stark (1992). These two components, referred as Reaction product A (RPA) and Reaction Product B (RPB), are assumed to be the result of Reactions (3) and (4) (toward the left). The stoichiometric coefficients were obtained by fitting experimental measurements of Si, Ca and R content of the ASR products, following a method similar to the one proposed by Helmuth et al. (1993). Note that the RPA has higher content of calcium and lower content of alkalis than the RPB. See Liaudat (2018) for further details. Coefficients x A and x B , which can be negative, make it possible to adjust the content of nonevaporable water of the RP. Since in the current version of the model the water balance is not explicitly considered, the stoichiometric coefficients x A and x B do not need to be specified.

Reactive transport processes at cracks and interfaces
In the SEM images of SL glass concrete affected by ASR (performed in the framework of the research project referred in the introduction Liaudat, 2018;Liaudat et al., 2019), ASR products could only be found filling cracks, both inside glass particles and in the cementitious matrix, and in between the cementitious matrix and glass aggregate particles. This, of course, does not imply that the precipitation of ASR products could not have occurred dispersed in the porous cementitious matrix, since, in this case, it would be probably impossible to detect in the SEM images. In any case, it seems clear that the precipitation reactions occur mainly at these locations. Moreover, the glass zones which appear to be altered by dissolution processes are reduced to a few microns inside of the original glass surface. For a detail discussion on these issues refer to Liaudat et al. (2019).
With this experimental evidence in mind, and with the aim of reducing the computational cost, in the proposed model Reactions (1)-(4) are assumed to occur only at the glass-cement interfaces and within opened cracks, i.e. dissolution/formation reactions are allowed within interface elements, but not in the continuum FEs. This assumption is also supported by simulations of ASR at a SL glass-cement interface performed with previous versions of the model (Liaudat et al., 2013). In those simulations, although the dissolution/formation reactions were allowed to occur anywhere in the continuum, they only took place along narrow bands of a few tens of microns on both sides of material interfaces. Note that modelling the processes occurring in this narrow band using continuum elements requires very fine meshes (element size % 1 lm) close to the glass-cement contacts (Liaudat et al., 2013) (and probably also close to potential glass-glass and cement-cement crack paths). This mesh refinement is substantially relaxed by the adopted assumption of restricting dissolution/formation reactions to occur only at interface elements, as it will be appreciated in following Sections 7 and 8.
The proposed reaction mechanism is schematically summarized in Fig. 1, where the interfacial zone between HCP and a SL glass particle, as well as a crack within that particle, are represented. As already mentioned, material pores and cracks are assumed to be fully saturated with water at all times during the reaction. Moreover, the pore water of the HCP has a certain content of aqueous calcium and alkalis, which determine the highly alkaline character of the pore solution, and which eventually diffuse into the pre-existing glass cracks. The SL glass in contact with this highly alkaline solution, both at the glass-cement contact and within the pre-existing glass crack, dissolves into silicate, calcium, and alkali ions (Reaction (2)). Wherever (within cracks or in between glass and cement) the three reactants (aqueous calcium, alkali and silicate ions) are available in sufficient concentration, Reactions (3) and (4) may occur forming RPA and RPB. Simultaneously, at the glass-cement interfacial zone (and eventually at both sides of HCP cracks), portlandite is dissolved (Reaction (1)) in order to restore the equilibrium calcium concentration of the pore solution.
Both, glass and portlandite dissolutions are assumed to occur as dissolution fronts moving away from the discontinuity (crack or interface) mid-plane and into the glass and the HCP, respectively. At the glass-cement interface and at HCP cracks, where calcium ions resulting from portlandite dissolution are more readily available, the formation of RPA is favoured in detriment of the formation of RPB. Consequently, the resulting ASR product has higher calcium content than the ASR products formed within the glass crack.
At a given location, if the volume of precipitated reaction products is greater than the available space within the crack/interface (plus the additional space liberated by silica and portlandite dissolution), a localized internal pressure is induced, which eventually may lead to the propagation of cracks. Additionally, the volume balance of the solid constituents plus the crack aperture, determines the transport properties of the crack/ interface.
It shall be noted that, as most studies in the literature, the mechanism described is mainly based on experimental expansion tests conducted at high temperature (60 and 80°C). Since temperature does not influence in the same manner the different diffusion-reaction processes involved, ASR expansions occurring in SL glass concrete at lower temperature might not exhibit exactly the same characteristics as the expansions occurring at high temperature. For a more extensive discussion on this topic, the reader is referred to Lothenbach and Leemann (2020). The model, as it is formulated in the following sections, would be able to predict these differentiated effects of temperature by merely introducing the temperature dependence of each model parameter (chemical equilibrium constants, diffusion coefficients, etc.). However, a deeper study of other qualitative effects of temperature is out of the scope of the present paper.

Mechanical properties of the ASR products
The ASR products are assumed to precipitate at the reaction site (where Reactions (3) and (4) occur) and remain there. That is, no flow of the ASR products due to applied stress is considered.
The volume occupied by a mole of solid ASR product 'a' (where 'a' may refer to either RPA or RPB) is referred as 'specific molar vol- Fig. 1. Schematic representation of ASR mechanism in SL glass concrete.
Joaquín Liaudat, I. Carol and C.M. López International Journal of Solids and Structures 207 (2020) 145-177 ume' (g a ). The 'apparent molar volume' (x a ) is defined as the sum of the volume of the solid ASR product plus the volume occupied by the absorbed water (gel water). The ratio between the volume of gel water and the apparent volume of ASR products is referred as 'gel porosity' (/ a ). In unstressed condition, the gel porosity and, consequently, the apparent molar volume have their maximum values / a o and x a o , respectively.
If a progressively increasing compressive stress r g (negative for compression) is applied to the ASR product, the gel porosity is assumed to monotonically decrease from / a o until zero for stresses beyond a certain threshold value r a th (Fig. 2, left). As the result, the apparent molar volume is monotonically reduced from x a o until reaching (if the solid deformations are neglected) the solid volume g a (Fig. 2, right). Note that parameter r a th determines the Maximum Swelling Pressure that can be exerted by the ASR product 'a' (Liaudat et al., 2019. In comparison with the volume changes induced by the variation of gel water content, the volume change of the solid particles due to applied stresses is regarded as negligible. The values of g a ; / a o , and r a th are assumed to depend on the composition of the reaction products, in such a way that lowcalcium products (RPB) have much higher values of / a o and r a th than high-calcium products (RPA). Finally, the properties of the mixture of RPA and RPB will result of the weighted average of the properties of each component.
The mechanical properties of the ASR products have been adopted on the basis of the reaction-expansion mechanism proposed in Liaudat et al., 2019, where a detailed discussion of the supporting experimental evidence can be found.

Diffusion-reaction problem
The formulation and numerical implementation of the diffusion-reaction problem may be structured in two parts, in analogy to traditional formulations of the mechanical problem. The first part deals with the general mass balance (continuity) equations in their discrete forms (Secs 3.1-3.4), while the second part deals with the particular 'constitutive' equations of the problem in hand (Sections 3.5 and 3.6). The general equations (in their generic version) would be valid not only for the ASR model proposed in this article but also for other diffusion-reaction problems in concrete (or another fractured porous medium). In contrast, the 'constitutive' equations, are developed for the particular case of the reaction-expansion mechanism for ASR in SL glass concrete proposed in Liaudat et al., 2019.

Solute mass balance
Solute mass balance equations are stated for the primary aqueous species (SiO 4 H 3 À , Ca 2+ , and R + ) in both the continuum porous medium and the discontinuities. The concentrations of the secondary aqueous species (OH À and H + ) are obtained as functions of the concentration of the primary species via chemical speciation calculations (Section 3.5.4). For the continuum porous medium, the approach proposed by Bear and Cheng for contaminant transport in groundwater (Bear and Cheng, 2010, Ch. 7) has been followed, whereas for the discontinuities, the formulation of Segura and Carol (2004) has been adapted making it consistent with the approach of Bear and Cheng. Second order effects in diffusion and chemical kinetics due to geometry changes induced by the mechanical loads are neglected for the porous continuum medium, but explicitly considered for discontinuities in terms of the normal aperture of the zerothickness interface elements.
The same formulation is used for simulating both types of discontinuities, cracks and ITZs. The only difference lies in the fact that cracks only appear when the mechanical strength of the material is overcome, while ITZ are always present in the HCP-aggregate interfaces. As developed below, from the formulation point of view that means that the discontinuity width (w) will be always w > 0 in the case of ITZ, whereas in the case of a crack it will remain w ¼ 0 while the mechanical strength is not overcome.
Finally, sink/source terms are introduced in the balance equations of interface elements to take into account the increase of the amount of aqueous species due to the dissolution of solids, and, vice versa, the reduction of the quantity of aqueous species due to the precipitation of solids.

Continuum porous medium
Let us consider a 2D domain X of continuum water saturated porous medium. The standard mass balance equation of the bspecies over this domain is obtained under the assumptions stated in Section 2 (Fick's diffusion, constant porosity, no chemical reactions), resulting in: where / [m 3 ps m À3 ] is the total porosity of the porous medium, c b [mol m À3 ps ] is the concentration of the b-species in the pore solution, expressed as quantity of solute per unit volume of pore solution, D b [m 2 s À1 ] is the isotropic effective diffusivity of b-species, and r Á ð Þ ¼ @ Á ð Þ=@x @ Á ð Þ=@y ½ T . Two possible boundary conditions are considered: i) Prescribed concentartion condition (Dirichlet condition): ii) Prescribed flux condition (Neumann condition): where c b is the prescribed concentration at the domain boundary @X c b and j b b is the prescribed normal flow at the domain boundary @X j b , with outer normal direction n. The last term in Eq. (8) expresses the solute diffusive flow through the border due to the concentration difference between the boundary c b ð Þ and the external media c b a À Á ; h [m s À1 ] is the 'convection' coefficient.
If now the three primary aqueous species of the problem are considered together and the dot notation ( _ X ¼ @X=@t) is used, the resulting system of equations can be written in matrix notation as follows: where The superscripts s; c, and r stand for the primary aqueous species SiO 4 H 3-, Ca 2+ , and R + , respectively. I 3 is the 3 Â 3 identity matrix.

Discontinuities
Let us consider a 2D discontinuity of width w [m] surrounded by continuum porous medium, and a local orthogonal coordinate system (l; n) defined on the plane tangent to the discontinuity surface (Fig. 3). The continuum and the discontinuity are assumed to be fully saturated with water. Diffusion of a generic aqueous b-species may occur both in the longitudinal and in the transversal directions to the discontinuity. Additionally, the bspecies may be consumed or produced locally due to chemical reactions.
For the sake of simplicity, it is assumed that the longitudinal diffusion along the discontinuity, as well as the chemical reaction rates, depend on the concentration of the b-species in pore solution at the centre (mid-plane) of the discontinuity width w and that is equal to: where c b top and c b bot [mol m À3 ps ] are the concentrations at both sides (top and bottom) of the discontinuity. Consistently, the remaining parameters/variables of the discontinuity are also computed at the mid-plane, but assumed to be valid across the discontinuity width w.
Then, considering a REV along the discontinuity with volume U ¼ wdl, and enforcing conservation of the solute mass we obtain: where w [m] is the discontinuity width, J b l [mol m À1 s À1 ] is the longitudinal local flow of b-species, j b bot and j b top [mol m À2 s À1 ] are leakoff terms incoming to the discontinuity from the surrounding continuum medium, / mp [m 3 ps m À3 ] is the total porosity of the discontinuity, i.e. the ratio between the void volume (saturated with pore solution) and the total volume of the discontinuity, and q b mp [mol m À2 s À1 ] is the chemical sink/source term of the b-species, which generally is a function of the concentration of the b-species itself and also of the concentrations of all the other aqueous species present in the pore solution. The effect of the time rates of w and / mp are neglected because it is assumed that the pore solution incoming into the discontinuity in order to keep it water saturated has the same concentration of b-species as the pore solution at the midplane.
The longitudinal local flow in Eq. (11) is obtained from the following generalized Fick's law: where D b l [ m 3 ps m À3 m 2 s À1 À Á ] is the effective longitudinal diffusivity of the discontinuity. Besides the longitudinal transmissivity, the existence of a discontinuity may also represent an obstacle or resistance to solute flow in the transversal direction, for instance due to the transition from a pore system into an open channel and back into a pore system or due to the existence of infill material. This resistance may reduce the diffusive flow in the transversal direction and results in a localized concentration drop across the discontinuity. In this formulation, a transversal concentration drop c b mp [mol m À3 ps ] between the discontinuity boundaries is considered, which is related to a transversal flux j where the subscripts bot and top stand for bottom and top discontinuity boundaries. If now Eqs. (11)-(13) are considered for the three primary aqueous species of the problem, the resulting system of equations can be written in matrix notation as follows:

Spatial discretization (FE formulation)
In this section, the above obtained governing equations are discretized in space by means of the Finite Element Method (FEM). For the sake of brevity, only the main steps of the derivation are described here, and a detailed description can be found in Liaudat (2018).

Continuum porous medium
The 2D domain X is discretized using finite elements with n nodes. The real distribution of concentrations c within this element is approximated by the following interpolation of the nodal values: where x and y are the global spatial coordinates, t stands for elapsed time, c i is the concentration vector at node i and N e i ¼ N i x; y ð ÞI 3 is the corresponding shape function matrix.
After applying the standard Galerkin's Weighted Residual method, the following matrix equation is obtained: where E e [m 3 s À1 ] is the diffusivity matrix, S e [m 3 ] is the capacity matrix, F e [mol s À1 ] is the nodal flux vector corresponding to prescribed flux over boundary @X j e , and f e [mol s À1 ] is the vector of nodal fluxes coming from neighbouring elements via shared nodes.

Discontinuities
For the FE discretization of discontinuities, zero-thickness interface elements with 'double nodes' proposed by Segura and Carol (2004) for diffusion problems are considered. These are interface elements with 4 or 6 nodes, depending on its linear or quadratic character respectively (see for instance the linear element sketched in Fig. 4).
In order to obtain the FE formulation for these elements, some auxiliary points are defined on the interface mid-plane and 'in between' nodes facing each other (Fig. 4), which in practice means the same coordinates for both if (as usual) the two nodes in the pair have same location. Concentrations at these mid-plane points are considered only as auxiliary variables, which will be finally eliminated in terms of the 'external' node variables, and therefore will not appear in the final element formulation. In the following, the formulation will be developed for linear elements although an analogous procedure would be valid for quadratic order elements.
The 'mid-plane points' define a 1D auxiliary element which is used to discretize the mid-plane of the discontinuity as sketched in Fig. 5a. The real distribution of concentrations along the discon-tinuity mid-plane (c mp ) is approximated through interpolation of the concentration values at the mid-plane points: where c mp i is the concentration vector at the mid-plane point i and N mp i ¼ N mpi l ð ÞI 3 is the corresponding shape function matrix. For a linear element, i.e. n mp ¼ 2; N mp ¼ N mp 1 N mp 2 Â Ã , and In order to express c mp in terms of the main variables at the actual nodes 1-4 ( Fig. 4), it is assumed that the concentration field along the mid-plane of the discontinuity is the average of the concentration field on the upper and lower sides of the interface, i.e.
which is a generalization of Eq. (10) for the three primary aqueous species. This assumption is considered to be valid all along the mid-plane including the mid-plane points and may be expressed as: where c j ¼ c 1 c 2 c 3 c 4 ½ T is the nodal concentration vector, and I 6 is the 6 Â 6 identity matrix. Introduction of Eq. (22) in Eq. (21) leads to the expression of concentration at any point of the midplane in terms of the nodal concentration values: Similarly, the transversal concentration jump along the discontinuity mid-plane ( c mp ) is approximated through interpolation of the concentration jumps at the mid-plane points ( c mp i ): where the shape function matrix (N mp ) is the same as in Eq. (21), and c mp is the vector of concentration jumps at mid-plane points.
In order to relate the nodal diffusion flows to the nodal concentration values in terms of the real nodes, the concentration jump at the mid-plane points are defined as the difference of concentration between the two surrounding real nodes: Then, introducing Eq. (25) in Eq. (24) one obtains: Now, following similar steps to those proposed by Segura and Carol (2004) for diffusion only problems on the basis of the analogy between the mechanical and the diffusion problems, and making use of the Principle of Virtual Work (Segura and Carol, 2004) the following FE equation for the diffusion-reaction problem is obtained: where E j [m 3 s À1 ] is the diffusivity matrix, S j [m 3 ] is the capacity matrix, Q j [mol s À1 ] is nodal flux vector corresponding to chemical source per unit area of discontinuity, and f j [mol s À1 ] is the vector of nodal fluxes coming from neighbouring elements via shared nodes.

Discontinuous porous medium
After assembly of the continuum and zero-thickness interface elements, the nodal flows f e and f j cancel out and a single system of equations in terms of the concentration of the primary aqueous species in pore solution at the nodes is reached for the whole problem domain. Assembly of Eq. (17) for the continuum elements and Eq. (27) for the interface elements gives the following system of equations for the boundary value problem under analysis after spatial discretization:

Time discretization
The generalized trapezoidal rule is considered in order to discretize Eq. (37) in time. If the problem is to be solved between an initial time t 0 and a final time t f , a partition of the time interval t 0 ; t f Â Ã is considered, which results in a series of time increments Dt nþ1 ¼ t nþ1 À t n . During a given time step Dt nþ1 , the variation of the nodal concentration vector c is linearly approximated by means of the following expressions: where h is a weighting factor that may vary between 0 and 1 depending on the time integration scheme considered. Analogous expressions are used for approximating E; S, and F in Dt nþ1 . The chemical source vector, in contrast, is obtained at time t nþh as a function of the linearly interpolated concentration vector c nþh , i.e.

Numerical solution
3.4.1. Residual vector and iteration strategy By considering function c t ð Þ to be an approximation to the solution of Eq. (37) in the time interval Dt nþ1 , then the corresponding residual vector is defined as: Note that these integrals represent the net amount of solute incoming/outgoing to the FE nodes in the time interval Dt nþ1 , the first integral corresponding to 'internal solute fluxes' (diffusion and chemical reactions), and the second integral corresponding to 'external solute fluxes' (imposed boundary conditions).
Then, using the time discretization proposed above with h ¼ 1 (Backward Euler scheme), and assuming that the values of E n ; S n ; F n ; F nþ1 , and c n are known, the residual vector i I for given trial values i E nþ1 ; i S nþ1 , and i c nþ1 may be approximated as follows: In order to improve the initial approximation to the solution, the Newton-Raphson method is applied, leading to the following system of linear equations: is the tangent 'stiffness' (jacobian) matrix. Solving this linear system of equations we obtain an improvement iþ1 dc nþ1 of the initial approximation such as iþ1 c nþ1 ¼ i c nþ1 þ iþ1 dc nþ1 . The calculation of @I Q =@c is discussed in the following section.

Vector I Q
In a general case, the global vector I Q will include the contributions of the local vectors I Q e (continuum elements) and I Q j (interface elements). However, the simplifying assumption of considering that no chemical reactions occur in the continuum medium (Section 2.5) implies that I Q e ¼ 0 for all the continuum elements. Therefore, only the expressions for I Q j are developed in the following. The analogous expressions for I Q e can be found in Liaudat (2018).
The most general expression of vector I Q j is: Introducing Eqs. (32) and (36), assuming linear interface elements, and doing some algebra, one obtains: where vector i q j is the time integral of the chemical source vector (q mp ) when a linear variation of the concentration vector c is assumed. For the three primary species of the model, the components of vector i q j are: The derivative of I Q j with respect to c j , needed for computing K (Section 3.4.1), is given by: By applying Leibniz's rule and the chain rule, and introducing Eq. (23), Eq. (47) can be rewritten as follows: Depending on the adopted chemical law, the components of vector i q j and of the 3 Â 3 Jacobian matrix @i q j =@c mp may be computed via closed-form expressions or may require to use numerical methods. For the ASR model proposed in this a paper, the components of i q j are calculated with a mixed analytical and numerical procedure, which is discussed in Section 3.5.5. Then, the components of @i q j =@c mp are obtained numerically via a forward difference implementation, which requires to compute the vector i q j three additional times per each time increment.

Implementation
The procedure implemented in the code DRACFLOW to solve the diffusion-reaction problem is summarized in Algorithm 1. The general scheme is similar to traditional implementations for mechanical problems. It consists of three main nested loops. The outermost loop (steps) allows the model to consider changing boundary conditions, as well as the update of the coupling variables (the mechanical nodal displacements u in the ASR model at hand). The second loop goes over the time increments assigned to each loading step. Within each time increment, first a trial value of the concentration vector is tested by computing the residuals vector (I). Then, if the euclidean norm of the residuals vector exceeds a certain tolerance value, an iterative Newton-Raphson procedure is activated in order to reduce the residuals, configuring the innermost loop. Each computation of the residual vector and of the tangent 'stiffness' matrix (K) requires the previous computation of intermediate matrices (E; S, and @I Q =@c) and vector I Q , which is achieved by means of two additional nested loops over the mesh elements and Gauss points. For the sake of brevity, the procedures in these loops have been only described in Algorithm 1 for the Algorithm 1. Pseudo-code for the diffusion-reaction problem 1: loop loading steps 2: Set boundary conditions 3: Compute F for the time increments in the loading step 4: if coupled analysis then 5: Get u at the beginning and at the end of the step (from the mechanical analysis) 6: Get u for each time increment by linear interpolation 7: else 8: u ¼ 0 9: end if 10: loop time increments 11: Get last converged c n 12: Set i ¼ 0 13: i c nþ1 ¼ c n

14:
Compute vector i I Q and matrices i E nþ1 , i S nþ1 , and i @I Q =@c À Á nþ1 (*) see below 15:

20:
Compute vector iþ1 I Q and matrices iþ1 E nþ1 , iþ1 S nþ1 , and iþ1 @I Q =@c À Á nþ1 (*) see below 21: end while 25: end loop time increments 26: end loop loading steps (*) Compute vector and matrices 27: loop elements 28: Extract c j j n from c n , and c j j nþ1 from i c nþ1 29: if coupled analysis; then Extract u j j n from u n , and u j j nþ1 from i u nþ1 ; end 30: loop Gauss points 31: Get last converged M mp j n 32: Get c mp j n and c mp j nþ1 interpolating nodal values from c j j n and c j j nþ1 33: if coupled analysis; then Get a cr j n from u j j n , and a cr j nþ1 from i u j j nþ1 ; end 34: call CHEMICAL LAW (IN: M mp j n , c mp j n , a cr j n , c mp j nþ1 , a cr j nþ1 ; OUT: M mp j nþ1 , i q j , @i q j =@c mp , U mp j nþ1 , U cp mp j nþ1 , / mp j nþ1 ) 35: call DIFFUSIVITY LAW (IN: U mp j nþ1 , U cp mp j nþ1 , a cr j nþ1 ; OUT: T l j nþ1 , D n j nþ1 ) 36: Compute contributions of the Gauss point to E j , S j , I Qj , @I Qj =@c j 37: end loop Gauss points 38: Assemble E j in i E nþ1 , i S j in i S nþ1 , I Qj in i I Q , and @I Qj =@c j in i @I Q =@c À Á nþ1 39: end loop elements interface elements in the mesh. For continuum elements, the procedures are analogous to those followed for the interface elements. However, for the particular case of the ASR model proposed in this paper, no chemical reactions are allowed in the continuum elements and, consequently, the chemical source terms are null and the diffusion coefficients remain constant.
Note that within the loop over the Gauss points the computation of the chemical source terms and of the related evolution of the amounts of solid species (saved in the local vector M) have been encapsulated in a subroutine called 'chemical law'. Similarly, the computation of the evolution of the diffusion coefficients with the solid volume fractions has been encapsulated in a subroutine called 'Diffusivity law'. In this way, different chemical and diffusivity laws can be implemented as a library and therefore the same code may be used for analysing different modelling assumptions or different diffusion-reaction problems.
The 'chemical' and 'diffusivity' laws for the ASR model proposed in this paper are developed in Sections 3.5 and 3.6, respectively, according to the modelling assumptions stated in Section 2.

Source terms
The net production rates at of primary aqueous species in the discontinuities are obtained by considering the stoichiometry and the reaction rates of Reactions (1)-(4) as follows: q r mp ¼ À 0:24 where the superscripts S; C; A and B indicate SL glass, portlandite, RPA and RPB, respectively, and C a mp [mol m À2 mp s À1 ] is the net rate of formation of a generic solid a-species per unit area of discontinuity.

Mass balance of solid species
The mass balance equation for solid a-species in a discontinuity is given by the expression: where M a mp [mol m À2 mp ] is the amount of solid a-species per unit area of discontinuity. The initial content of a-species in the interface elements is always null, i.e. M a mp t ¼ 0 ð Þ¼0. However, depending on the type of discontinuity being represented, dissolution of SL glass and/or portlandite is allowed leading to negative contents. A negative M S mp or M C mp indicates the amount of SL glass/portlandite dissolved at one or both sides of the interface. Therefore, M S mp can only be negative if at least one side of the interface element is shared with a continuum element representing SL glass. Similarly, M C mp can only be negative if at least one side of the interface element is shared with a continuum element representing cementitious matrix (mortar or HCP).

Volumetric balance
The initial volumetric fractions of solid species in the continuum medium will remain constant, since precipitation/dissolution reactions are not allowed to occur there. In contrast, the precipitation/dissolution of solid species is allowed at discontinuities (glasscement contacts or cracks) and consequently the volumetric fractions assigned to them may change.
The volumes of solid species considered for discontinuities are obtained as functions of their molar concentration as follows: where g a [m 3 mol À1 ] and / a [m 3 ps m À3 ] are the 'specific' molar volume and the gel porosity of the solid a-species, respectively.
The gel porosity is generically assumed to depend on the applied stress r g , according to: where / a o is the gel porosity value for the unstressed state. f a / r g À Á is a continuous, non-negative function defined for r g 2 À1; 0 ð as follows: where p > 1 and k-0 are dimensionless shape parameters, and r th [Pa] is the threshold stress below which the gel porosity becomes null (see Fig. 2). Note that f a / r g À Á is an interpolation function ranging from 1, for r g ¼ 0, to 0, for r g 6 r th . Since SL glass and portlandite have no gel porosity, function f a / is only applicable for Reaction Products A and B. The same threshold stress is used in functions f A / and f B / , which is obtained as follows: where r A th and r B th are the threshold stresses of pure Reaction Products A and B, respectively.
The 'specific' molar volume g a is treated differently in the diffusion-reaction formulation than in the mechanical one. In the first case, it is considered to remain constant, being unaffected by r g , i.e. SL glass, portlandite and the reaction product particles are treated as incompressible solids. This is justified because the volume change of these solids due to the applied stress is expected to be much smaller than the overall volume changes induced by the reduction of the gel water content. Therefore, neglecting the change of specific molar volume is not expected to significantly modify the chemo-transport results.
In the case of the mechanical formulation, in contrast, considering the reaction products filling a discontinuity as incompressible solids is a very restrictive condition which would hinder the numerical solution of the problem. This problem is avoided by assigning, in the mechanical formulation, the following stress dependence to the specific molar volume of Reaction Products A and B: where g a o is the specific molar volume in unstressed state and f a g r g À Á is a continuous, non-negative function defined for r g 2 À1; 0 ð as follows: where k g [Pa À1 ] is a positive (very small) compliance constant. Since r g 0, the smaller k g , the higher the stiffness of the reaction product particles. Note that f a g ranges from 1, for r g ¼ 0, to 0, for r g ! À1. k g plays the role of a 'penalty coefficient' that needs to be sufficiently small to induce negligible volume changes from the chemo-transport point of view, but sufficiently high to facilitate the converge to the numerical solution of the mechanical problem. This issue is further discussed in Section 4.2.2.
As explained above (Section 3.5.2), M S mp and M C mp (and consequently U S mp and U C mp ) can be negative, indicating the amount of SL glass/portlandite dissolved at one or both sides of the discontinuity. Portlandite is finely dispersed in the HCP and, consequently, its dissolution leaves behind a solid skeleton of inert cement paste which continues delimiting the discontinuity. In contrast, SL glass dissolution, being a congruent dissolution process, moves forward the border of the discontinuity. In a simplified approach, these observations lead us to consider the volume of dissolved SL glass (but not the volume of dissolved portlandite) as part of the total volume of discontinuity. Then, the total volume and the volume of 'capillary pores' (volume filled with free water) assigned to a discontinuity are, respectively: where X h i À and X h i þ are Macaulay brackets standing for 'the negative part of X' and 'the positive part of X', respectively, U o mp [m 3 m À2 mp ] is an initial volume of capillary pores, and a cr n [m] is the cracking normal aperture. The latter is a coupling variable obtained from the mechanical analysis, as it is described below (Section 4), which is updated at each iteration of the staggered procedure.
Finally, under the assumption of water saturation, the total volume of pore solution (U ps mp ) and the total porosity in the discontinuity (/ mp ) are given by: In order to facilitate the interpretation of previous equations, Fig. 6 depicts a schematic representation of a SL glass-HCP interface under different conditions. The first case (Fig. 6a) represents a sound (not cracked) interface before any dissolution/formation reaction occurs. In the proposed modelling framework, reactions in discontinuities can occur only if they contain certain amount of pore solution (/ mp > 0). Therefore, in order to simulate the ASR occurring at a SL glass-HCP interface, an initial volume of pore solution needs to be considered in the corresponding interface elements. This is done by assigning an initial volume of 'capillary pores' (volume filled with free water) U o mp . This initial volume is only considered in the diffusion-reaction problem, i.e. for the mechanical problem the interface is considered closed (a cr n ¼ 0). Note that U o mp will be null for interface elements inside bulk glass or cementitious matrix.
The second case (Fig. 6b) represents the same interface with its faces separated a distance a cr n after a fracturing process. The total volume of the interface in this condition is U mp ¼ a cr n þ U o mp , which coincides with the volume of capillary pores (U cp mp ) since the solid fractions are null.
The third case (Fig. 6c) represents a sound (not cracked) interface where dissolution of glass and portlandite, as well of precipitation of reaction products, have occurred. Dissolution fronts of glass and portlandite are indicated with dotted lines in the figure.
Since glass is assumed to dissolve congruently, the advancement of the dissolution creates a volume of 'capillary pores' which is incorporated to the interface volume. In contrast, the dissolution of portlandite in the adjacent HCP increases the porosity but leaves a residual solid skeleton which continues delimiting the discontinuity. Therefore, the volume of the dissolved portlandite is not incorporated to the interface volume. Simultaneously to the dissolution of glass and portlandite, reaction products have been formed filling in part the 'capillary pores'.
As the reaction progresses, eventually the volume of 'capillary pores' is exhausted and the precipitation of any additional amounts of reaction products will exert an internal pressure at both sides of the interface. The mechanical reaction to this pressure (r g ) will reduce the gel porosity and, consequently, the apparent volume of the reaction products. If the internal pressure is high enough (Fig. 6d), the tensile strength of the interface may be overcome, and the interface is cracked developing a normal aperture a cr n . The volume created by the crack aperture is incorporated to the total interface volume as capillary pore volume. The crack aperture decompresses the reaction products in such a way that the newly created capillary pore space is immediately filled by them. On the other hand, if the interface cannot crack either because of high strengths or because of normal compressive stresses, the gel porosity is progressively reduced and eventually the reaction progress may be stopped altogether.

Reaction rates
The calculation of the sink/source terms q b mp according with Eqs. (49)-(51) requires establishing the kinetic laws of Reactions (1)-(4) in order to obtain the corresponding reaction rates where C a mp [mol m À3 ps s À1 ] is the net formation rate of a-species per unit volume of pore solution in the discontinuity.
It is assumed that the driving force of the dissolution/formation reactions of reactive a-species is the distance from equilibrium expressed as w a À 1 ð Þ , where w a is the dimensionless saturation index of the pore solution with respect to the solid a-species. Dissolution and precipitation reactions are assumed not to occur simultaneously but alternately depending on w a . If w a > 1, the solution is over-saturated with respect to the a-species and, consequently, the reaction progresses in the precipitation direction. If w a < 1, the solution is under-saturated and solid dissolution occurs. If w a ¼ 1, the solid and the solution are in thermodynamic equilibrium. The saturation indexes adopted for Reactions (1)- (4) are: where superscripts h and oh stand for H + and OH À , respectively, a b [mol m ps

À3
] is the thermodynamic activity of the aqueous b-species in pore solution, and K a sp is the saturation product constant of aspecies dissolution. Note that the saturation index of SL glass has been assumed to depend only on the activity of silicate ions and the pH of the solution, as it has been proposed by Maraghechi et al. (2016). The saturation product constants of portlandite (K C sp ) and of SL glass (K S sp ) as functions of temperature may be obtained from the literature as it is detailed in Liaudat (2018). The saturation product constants of RPA and RPB (K A sp and K B sp ), in contrast, need to be calibrated by fitting ASR experimental results.
Then, the above saturation indexes are used for defining the following kinetic laws: where k a f and k a d are kinetic constants for the a-species formation and dissolution, respectively. The kinetic laws are usually consid-ered to be proportional to the distance w a À 1 ð Þ , e.g. (Icenhower and Dove, 2000;Galí et al., 2001;Maraghechi et al., 2016). However, from the numerical point of view, this has the disadvantage that the maximum precipitation rate is not bounded. It has been found that in some cases this may lead to divergence of the Newton-Raphson iterative procedure. In order to avoid this problem, the reaction rates which have been adopted are proportional to the following exponential function of the saturation index w a : This function varies monotonically from f w ¼ À1, for w a ¼ 0, to f w ¼ 0, for w a ¼ 1, and asymptotically f w ! 1 for w a ! þ1. In all cases, the formation reactions are considered to occur uniformly in the pore solution within the interface FE and, therefore, the kinetics of the a-species formation has been made proportional to U ps mp . Having no experimental data on the specific surface of ASR products and for the sake of simplicity, the dissolution reactions of RPA and RPB are also considered to be proportional to U ps mp . In contrast, as mentioned above (Section 2.5), the dissolution reactions of SL glass and portlandite are considered to occur in dissolution fronts, and, therefore, to have dissolution rates proportional to the surface area of the corresponding dissolution front. This surface area changes depending on the location of the interface element (at the cementitious matrix, inside a glass particle, or at a glass-cement contact) and is not constant but evolves with the cracking of the interface and the precipitation of solids within it. These effects have been introduced in the kinetic laws of SL glass and portlandite dissolution by means of the dimensionless factors is the specific fracture energy in mode I, and W cr [J m À2 mp ] is a history variable of the mechanical constitutive law of the interface which saves the amount of work spent on fracture processes during the formation of the crack (see following Section 4.2.1).
Finally, the dissolution rate of SL glass is affected by two additional factors. The first one is the penalization factor f S 2 which asymptotically tends to 0 for growing amounts of dissolved glass where A S [m 2 mol À1 ] is a fitting parameter. This factor accounts for the increasing distance of the glass dissolution front to the midplane of the discontinuity. The second factor, the activity of the hydroxyl ions with an exponent p > 0 (Eq. (68)), intensifies the effect of the pH of the pore solution on the dissolution rate and it is inspired in Maraghechi et al. (2016). Thermodynamic activity In order to account for deviations from ideal behaviour of the aqueous b-species in pore solution due to the interactions with other aqueous species, activity a b is considered in chemical equilibrium equations instead of concentration c b . Activity and concentration are related by means of a dimensionless factor c b (activity coefficient) as follows The activity coefficients are calculated using Eq. (76), where z b is the number of electric charges of the b-species considered, I [mol/L] is the ionic strength of the solution, and A c is the Debye-Hückel parameter (at 60°C, A c ¼ 0: 546 Pitzer, 1991). This equation is the well-known Davies Equation with the modification of the second term proposed by Samson and Lemaire (1999), in order to extend its applicability to solutions at ionic strengths up to 1.20 mol/L.
Chemical speciation The concentration of secondary species in pore solution is determined locally as a function of the concentration of primary species. This is achieved by simultaneously enforcing the following equations: The first equation enforces the electrical charge neutrality of the solution. The second one is the equilibrium equation of water selfdissociation (Reaction (5)), where K w eq is the water self-ionization constant.
Considering previous Eqs. (77) and (78) plus Eqs. (75) and (76), it is possible to determine the equilibrium concentrations of both secondary species as a function of the concentrations of the primary ones by solving the resulting nonlinear system of equations. This is achieved by a mix of analytical and numerical procedures. First, a unique polynomial equation of second order in terms of c oh has been obtained by assuming constant values for activity coefficients, i.e. ignoring Eq. (76), and operating with the rest of the equations. Then, trial values of the activity coefficients are adopted and the polynomial equation is solved analytically. Once the value c oh is known, the c h is readily obtained from Eq. (78), and the corresponding, new, activity coefficients are calculated with Eq. (76). If the new activity coefficients differ less than a pre-established tolerance from the initially adopted ones, the calculation concludes. Otherwise, the procedure is repeated until convergence is reached.

Time integration of the source terms
In order to obtain the components of the vector i q j the time integrals between t n and t nþ1 of the reaction rates of the solid species need to be calculated (Eqs. (46)-(49), (47)-(51)). Since according to Eq. (52), and since each reaction rate is a function of the amounts of the four solid species (M S ; M C ; M A ; M B ) through the total porosity / (see Sections 3.5.3 and 3.5.4), these time integrals cannot be calculated separately.
To solve the resulting system of equations, a mixed numericalanalytical procedure is proposed, which is based on assuming that (i) the reaction rates during the considered time lapse remain constant and equal to their respective initial values (explicit forward Euler approximation), and (ii) the crack normal aperture (a cr n ) and the concentration of the primary species (c) vary linearly within the sub-increment. With these assumptions, it is possible to analytically calculate the local solid concentrations (M a ) at the end of the considered time lapse. The error introduced by the second assumption, can be reduced as much as required by subdividing the time increment D nþ1 into smaller sub-increments, as it is described below. Further details of this procedure can be found in Liaudat (2018).

Implementation
The procedure to compute the chemical constitutive law is summarized in Algorithm 2. It consists of two nested loops. The inner loop calculates de evolution of the solid species during the time increment Dt nþ1 . To do so, Dt nþ1 is subdivided into N equal sub-increments. The reaction rates are evaluated at the beginning of the sub-increment according to the expressions in Section 3.5.4, and the mixed numerical-analytical procedure described in Section 3.5.5 is used for the time integration. Once the vector of local i T has been obtained, the vector i q j and the volumetric fractions are calculated.
The outer loop is implemented for computing the Jacobian matrix @i q j =@c mp with a forward difference procedure. In the first run of the loop, i q j is computed using vector c as it is introduced in the subroutine ( c nþ1 ¼ c nþ1 ). In the following three runs, a slight concentration increment k is selectively applied to each component (species) of vector c nþ1 obtaining in this way three additional Once out of the loop, @i q j =@c mp is approximated with a forward difference scheme, by using the four vectors i fd ð Þ q j as indicated in the algorithm.
Algorithm 2. Pseudo-code for the chemical constitutive law.
Subindexes mp have been removed for clarity.
procedure CHEMICAL LAW IN: M n , c n , a cr n j n , c nþ1 , a cr n j nþ1 , Dt nþ1 ; OUT: 3.6. Diffusivity constitutive laws 3.6.1. Continuum porous medium As stated above (Section 2.3), the effective diffusivity of primary aqueous species in the continuum medium (D b ), as well as its porosity /, are assumed to remain constant. Since the direct introduction of these conditions as D b ¼ 0 and / ¼ 0 would cause indetermination of the global diffusivity and capacity matrices, an internal algorithm was implemented in the diffusion-reaction code (DRACFLOW) that, besides assigning D b ¼ 0 and / ¼ 0 to the SL glass elements, identifies the problematic nodes and imposes to them null concentrations (c ¼ 0). In that way, the DOFs corresponding to these nodes are 'excluded' from the solution of the system of equations corresponding to the diffusion-reaction problem.

Discontinuities
The longitudinal (T b l ) and transversal (D b n ) diffusivities of discontinuities will vary as ASR progresses, due to dissolution/formation of the solid species assigned to them, and also due to the possible changes in the mechanical aperture due to fracturing processes. These effects are introduced by means of the following expressions: are the diffusion coefficients of the b-species in the filled sector of the discontinuity, in bulk water, and in the Interfacial Transition Zone (ITZ), respectively; w crit % 1 Â 10 À4 m is a critical aperture above which the discontinuity walls have no effect on the diffusivity of the b-species (Idiart et al., 2011); and w ITZ [m] is the width of the ITZ, which is greater than zero for interface elements located in between SL glass and the cementitious matrix and equal to zero otherwise. The derivation of Eqs. (80) and (81) can be found in Liaudat (2018).

General description
The mechanical problem is treated with an existing model which uses the same FE mesh as the previously described diffusion-reaction model. In this mechanical model, tractionseparation constitutive laws based on principles of non-linear fracture mechanics are used for the interface elements in order to represent the main potential crack paths, while the bulk aggregate and the bulk cementitious matrix are considered to be linear elastic. The model formulation is based on the criterion that tension is positive and compression negative, and on the assumptions of smallstrain theory, isothermal equilibrium and negligible inertial forces. With this conceptually simple model, the mechanical behaviour of concrete has been successfully reproduced under a variety of uniaxial, biaxial and triaxial conditions (Carol et al., 2001;Caballero et al., 2006;Caballero et al., 2007;López et al., 2008a;López et al., 2008b). This model is based on a well established mechanical formulation for zero-thickness interface elements with double nodes within the discrete crack approach (Carol et al., 1985;Gens et al., 1990;Carol et al., 1997;Carol et al., 2001).
In comparison with previous implementations, the only substantial change is the modification of the constitutive law for inter-face elements in order to consider the effect of the ASR products formed within discontinuities, as it is discussed below.

Constitutive law for interface elements
Following the concepts proposed in Ulm et al. (2000), the mechanism of ASR expansion of a discontinuity is schematically represented in Fig. 7a. The precipitation of reaction products within the discontinuity after the volume of 'capillary pores' (U cp mp ) has been filled up, induces a local pressure r À g that tends to separate the discontinuity sides. As a consequence, the solid bridges that keep both sides of discontinuity together are tensioned developing local stress r À s . The average of these local pressure/stress in the discontinuity surface area are referred to as r g and r s , respectively.
The conjugate variable for both stresses is the relative displacement of both sides of the discontinuity a mp . That means that the reaction products and the discontinuity interact in a parallel scheme and, therefore, the constitutive stress vector of the discontinuity filled with reaction products is given by: By differentiating previous equation, and taking into account that r g depends on the aperture as well as on the amount of reaction products accumulated within the discontinuity, one obtains: where concentrations. The mechanical constitutive model considered for a discontinuity filled with reaction products is schematically represented by the parallel device shown in Fig. 7b. The left branch represents the constitutive law of the discontinuity, while the right branch represents the constitutive law of the reaction products within it. The constitutive law considered for the discontinuity, briefly described in the Section 4.2.1, is the one proposed in Carol et al. (2001) and López et al. (2008a) without introducing any modification. The constitutive law of the reaction products within the discontinuity (Section 4.2.2), in contrast, has been newly derived from the balance of volumetric fractions proposed in Section 3.5.3.

Constitutive model for cracks
The constitutive model for discontinuities is formulated in terms of the normal and tangential stress components in the mid-plane of the element r s ¼ r n s l ½ T and the corresponding relative displacements a mp ¼ a n a l ½ T . The formulation borrows the mathematical structure from the theory of elasto-plasticity, and introduces nonlinear Fracture Mechanics concepts in order to define the softening behaviour due to the work dissipated in the fracture process (denoted as W cr ). The elastic stiffness matrix is diagonal and the elastic constants K n and K l can be regarded simply as penalty coefficients. This means that their values should be as high as possible in such a way that interface elements in the elastic regime, before cracking, do not add significant elastic compliance (Carol et al., 1997). The upper bound of K n and K l is set as a compromise between the added compliance and the numerical illconditioning found for very high values.
The fracture surface F ¼ 0 is defined by a hyperbola of three parameters in the stress-traction space given by: where v [Pa] is the vertex of the hyperbola representing the tensile strength, c [Pa] is the cohesion, and u is the asymptotic friction angle (see curve '0' in Fig. 8).
Cracking starts as soon as the fracture surface is reached. Once the fracture process has started, the failure surface contracts due to the decrease of the main parameters according to some evolution laws based on the work dissipated in the fracture process (W cr ). In order to control the fracture surface evolution, the model introduces two additional parameters: the fracture energy in mode I, G f i (pure tension), and a second fracture energy in mode IIa, denoted as G IIa f , defined by a shear state under high compression, so that dilatancy is prevented (Carol et al., 1997). Total exhaustion of tensile strength (v ¼ 0) is reached for W cr ¼ G i f (curve '1' in Fig. 8), and the residual friction state (c ¼ 0 and u ¼ u r ) is attained when Fig. 8), state for which the hyperbola has degenerated into a pair of straight lines representing pure residual friction. A non-associated plastic potential Q ¼ constant is adopted in order to progressively eliminate dilatancy for high levels of compression (dilatancy vanishes progressively for r n ! r dil ). The dilatancy angle also decreases during the fracture process, vanishing for W cr ¼ G IIa f . The numerical implementation is based on an implicit backward Euler integration scheme with a consistent tangent matrix operator, as described in Caballero et al. (2008).
The constitutive behaviour of the model may be directly compared with experiments to determine some of the parameters (v o ; c o ; u o ; u r ; G I f ) as it was done for the cases of uniaxial tensile tests and mixed mode shear/compression tests in concrete in Carol et al. (2001), Carol et al. (1997), Caballero et al. (2008) and López et al. (2008b). Regarding the rest of parameters (G IIa f ; r dil ), they can be estimated from inverse analysis when fitting the overall response of concrete specimens subject to compression and tension (especially including the softening branch) with a meso-scale model including zero-thickness interface elements (Carol et al., 2001;Caballero et al., 2006;Caballero et al., 2007;López et al., 2008a;López et al., 2008b).

Constitutive law of the reaction products
A mechanical constitutive law for the reaction products filling a discontinuity is derived in this section from the balance of volumetric fractions proposed in Section 3.5.3. The precipitated reaction products are assumed not to exhibit any Poisson's effect and to have null tangential stiffness. Out of the two components of the general strain/stress vector in the interface local axes n; l ð Þ, only the component in the normal direction is assumed to play a role. In this way, the stress state of the reaction products can be suitably presented in the same stress space used for the formulation of zero-thickness interface elements with the vector r g ¼ r g 0 ½ T .
Moreover, the precipitation of the reaction products is assumed to occur homogeneously in the normal direction of the discontinuity, and, therefore, the reaction products stress r g and the conjugated reaction products strain e g are assumed to be uniform in this direction. These assumptions have two main consequences: (i) the strain in the normal direction e g coincides with the volumetric strain of the reaction product; (ii) for a given content of reaction products, e g can be univocally determined by the relative normal displacement of the interface a cr n . To establish this last relation, it is convenient to define the auxiliary variable U max [m 3 m À2 mp ] as follows: Obtained from Eq. (53), U max is the volume occupied by the reaction products filling the discontinuity in unstressed state, and, therefore, it is the reference volume from which the strain of the reaction products has to be computed, i.e.: The Macaulay Brackets are used here for considering the case in which the volume U of the discontinuity is greater than U max . Since the reaction products are assumed not capable of developing tensile strain/stress, U > U max implies that U cp > 0 and, therefore, that the reaction products are unstrained. Otherwise, if U 6 U max , then U cp ¼ 0 and the reaction products have been compressed (e g < 0 and r g < 0) until the volume of reaction products and the volume of discontinuity were equal. In this last condition and based on Eqs. (59) and (53), the volume U can be expressed either as Note that in Eq. (87) it has been assumed that the elastic relative displacements are negligible and, consequently, a cr n % a n h i þ . Replacing these latter equations in Eq. (86) and for a given set of values M a mp , the strain of the reaction products can be obtained indistinctly either as a function of a n , or as function of r g , The first equation establishes a 'compatibility condition' by relating the reaction products strain to the relative normal displacement of the interface element and, indirectly, to its nodal displacements. The second equation, in contrast, establishes the mechanical 'constitutive relation' between strain and stress of the reaction products.
In order to combine the mechanical response of the discontinuity and the reaction products in the above-proposed parallel scheme (Fig. 7), it is convenient to obtain an expression to calculate r g as a function of a n . This is achieved by combining Eqs. (89) and (90) as follows: From this implicit equation, depending on the adopted functions f a / and f a g , the calculation of r g for a given a n may be done analytically or may require to use numerical methods.
The tangential stiffness matrix of the constitutive equation of the reaction products within the discontinuity is defined as follows: Applying the chain rule, D g can be rewritten as: The derivatives on the right-hand side of this equation are obtained by differentiating Eqs. (89) and (90) as follows: Finally, introducing Eqs. (94) and (95) in Eq. (93), one obtains the final expression of the tangential stiffness: In order to illustrate the resulting constitutive behaviour, Fig. 9 depicts the curves r g -e g and r g -a cr n , corresponding to a crack without initial volume (U o mp ¼ 0) and with a certain content of Reaction Product B (RPB) (M B mp > 0). At point r, crack aperture is greater than U max , and, consequently, the RPB is unstrained/unstressed. At point s, crack aperture is smaller than U max but greater than the 'specific' volume of RPB. Since the strain of solid clusters is negligible (k g % 0), this volume reduction is due to the reduction of gel water content, reflected in the change of gel porosity from / B o to a smaller value / B . Finally, at point t, stress is higher than the maximum swelling pressure of RPB (r g > r B th ) and, consequently, the gel water has been completely driven out (/ B ¼ 0). The tangential stiffness at this point is very high, since any further reduction of the RPB volume will occur at the expense of deforming the solid clusters. An extended discussion on the effect of the different parameters of the constitutive model can be found in Liaudat (2018).

Chemo-mechanical coupling
As mentioned before, the mechanical formulation and the diffusion-reaction formulation are implemented in two different codes, called DRAC and DRACFLOW, respectively. Both codes use the same FE mesh and the same integration points. Moreover, both codes have the same time discretization scheme consisting in time 'steps' composed each of them by a number of time 'increments'. In order to establish the coupling staggered scheme, the time steps in DRAC and DRACFLOW have to match each other (Fig. 10). However, the number of increments within the time steps may be different for each code, depending on the requirements of the particular problem to be modelled.
The staggered procedure is administrated by a third code named STAG, schematically illustrated in Fig. 11. For the calculation of each time step, DRAC and DRACFLOW are called alternatively Fig. 9. Mechanical constitutive curves of Reaction Product B filling a crack. Joaquín Liaudat, I. Carol and C.M. López International Journal of Solids and Structures 207 (2020) 145-177 in order to independently solve the mechanical problem and the diffusion-reaction problem, respectively. After each call, the mechanic code DRAC writes a file with the nodal displacements obtained at the end of the time step. This file is read by DRACFLOW and the displacement values are used to calculate the crack normal aperture a cr n of the interface elements at each integration point at the end of the current step. The value of a cr n at each increment within the time step is linearly interpolated between the values at the beginning and at the end of the step. This value is then used in the chemical constitutive law (Section 3.5) for calculating the evolution of the amount M a of each solid a-species, as well as of the associated volumetric fractions U mp ; U cp mp , and U ps mp . These volumetric fractions, in turn, determine the transport properties of the interface elements (/ mp ; T l , and D n ), via Eq. (62) and the diffusivity law (Section 3.6). On the other hand, after each call, The nodal concentration files and the nodal displacement files are used by STAG to assess the convergence of the staggered procedure by comparing their values obtained at two successive iterations of the staggered procedure (Fig. 11).
6. Summary of model parameters i) Parameters that can be obtained from chemistry manuals or geochemistry databases as functions of temperature: K w eq ; K S sp ; K C sp ; g S ; g C ; A c ; D s w ; D c w ; D r w . ii) Parameters that can be directly measured in standard laboratory tests or estimated using expressions from the literature as functions of the chemical and mix composition of the cementitious iii) Parameters that can be determined by purely mechanical inverse analysis or direct comparison with standard compression  (68)). v) Penalty coefficients that should be as high as possible without causing numerical ill-conditioning: K n ; K l ; k À1 g .

The experiments
The Interfacial Expansion Tests (Liaudat et al., 2019) aim at facilitating the study of ASR-induced expansions by reducing the complexity of concrete meso-structure. For this purpose, instead of studying ASR in regular concrete specimens, a new type of specimen with a single, disc-shaped reactive aggregate is used. In these specimens, the reactive aggregate disc is placed as a sandwich in between two layers of cementitious matrix. In this way, two plane, symmetrical and geometrically defined interfaces between the reactive aggregate and the cementitious matrix are obtained. See Fig. 12a for a schematic representation of the reactive specimens. In this configuration, the ASR products mainly precipitate along these interfaces, tending to separate the aggregate disc from the cementitious matrix. In order to accelerate the ASR expansion   rates, the specimens are exposed to an alkaline solution at 60°C after a curing period of 28 days at room temperature. Besides the specimens with reactive aggregate disc, similar specimens without reactive aggregate are simultaneously tested under the same exposure conditions, in order to assess the deformation occurring in the cementitious matrix. With these tests, it is expected that an estimate may be obtained of the expansion developed at a single glass-cement interface, by measuring the length change of the reactive specimen and subtracting the length change of the cementitious matrix tested simultaneously. Furthermore, thanks to the simplicity of the geometry it is hoped that the measured expansion may be related to the thickness and composition (measured by means of SEM and EDS analysis) of the layer of ASR products. See Liaudat et al. (2019) for further methodological details.
Preliminary results obtained with this methodology suggested that the mass exchange between the alkaline bath and the specimens may play a significant role on the development of the interfacial expansions (Liaudat, 2018). In order to further assess this possible relation, a series of tests was performed with two different sealing schemes, one with the lateral surface of the specimens completely sealed (Fig. 12b), and the other only covering the interfacial zone (Fig. 12c). These tests were performed as part of the bachelor's thesis of Martínez (2014) under the supervision of the authors. An essential underlying assumption in the approach described is that the glass-matrix interfaces constitute preferential paths for water and/or ion fluxes due to the higher permeability of the ITZ (Scrivener and Nemati, 1996). Consequently, one would expect that the specimens with the interfacial zone sealed expand similarly to those with complete lateral sealing.
For these tests, specimens made with SL glass as reactive aggregate and plain cement paste as cementitious matrix were used (Fig. 12a). The Portland cement was CEM I 42,5 N-SR5 (AENOR, 2011) and the water/cement (w/c) weight ratio was 0.47. Sodium hydroxide was added to the mixing water in order to rise the alkali content in terms of equivalent sodium oxide (Na 2 O e ) up to 1.1 wt% of the cement weight. Once the curing period was completed, the specimens were sealed laterally according to the scheme shown in Fig. 12, by using a self-curing rubber tape. The alkaline solution for the curing and exposure baths consisted of 40 g NaOH and 2 g of Ca(OH) 2 per kg of water. Further details of the preparation of the specimens and curing and exposure conditions can be found in Liaudat (2018).
The interfacial expansions developed during the exposure stage by the specimens of the three sets are given in Fig. 13, in terms of the absolute expansion of a single glass-cement interface d I (Liaudat et al., 2019). The tests of Set #1 ended when one of the glass-cement interfaces got detached due to the internal pressure exerted by the ASR products. The tests of Sets #2 and #3 ended after two months of exposition, except for specimen F17 which was prematurely retired at day 21 because its sealing got loose. In the plots of Fig. 13, the effect of lateral sealing on the development of the interfacial expansion can be appreciated. As the degree of lateral sealing increases, the expansion rate decreases. The mechanisms behind this effect are assessed by numerically simulating each of the three sets of Interfacial Expansion tests as it is described below.

Model description
The model geometry is schematically represented in Fig. 14a. The geometry considered takes advantage of the specimen symmetry: a horizontal symmetry plane normal to the axis of the specimen, and a vertical symmetry axis coincident with the specimen axis. The latter assumption allows us to reduce the 3D problem to a 2D axisymmetric problem. (The FE formulation of the model for axisymmetric problems is not explained in this paper, but it can be readily obtained by analogy with axisymmetric formulations from the literature, e.g. Zienkiewicz et al., 2013). This geometry was discretized with 378 quadrilateral quadratic continuum elements and 14 quadratic interface elements (Fig. 14b), giving a total of 1246 nodes. The interface elements were placed in between the continuum elements representing the HCP and those representing SL glass.
The examples in this section (as well as in the following Section 8) have been analysed with different FE meshes. The mesh finally presented here is the one leading to the best balance between improving the numerical convergence, the stability of the results (for the diffusion-reaction problem), and the spatial resolution of the field results, which usually require finer meshes, and the limitation of the computational cost which grows rapidly with the number of nodes. In the adoption of the FE mesh, one should be aware that, in general, the optimum mesh for the mechanical will be different from the optimum mesh for the diffusion-reaction problem. An advantage of the proposed model is that it is automatically regularized from the point view of objectivity of the mechanical (fracture) problem with mesh refinement. Intrinsic regularization is due to using a cohesive-type constitutive relation to represent fracture, which is formulated in terms of relative displacements (rather than strains). In this type of formulations, the fracture energy parameter of the interface constitutive model is directly G f , energy per unit of crack area (rather than energy per unit of 'cracked volume', g f , in continuum-type softening models).
The time discretization of the staggered scheme consisted of 896 steps of 1.5 h each. These steps were further subdivided in 12 and 100 equal increments for the diffusion-reaction analysis with DRACFLOW and for the mechanical analysis in DRAC, respectively.
The boundary conditions for the mechanical problem, intended to materialize the vertical symmetry axis and the horizontal symmetry plane, are schematically represented in Fig. 14c.
The parameters of the mechanical constitutive law of the interface elements (Section 4.2.1) are given in part in Table 1. The rest of the parameters were: G IIa f ¼ 10 Â G I f ; tan u r ¼ tan u o ; r dil ¼ 40 MPa. The continuum elements were assumed to be linear elastic, with elastic modulus E and Poisson's coefficient m given in Table 2.
The boundary conditions for the diffusion problem are schematically represented in Fig. 14d, were the thick lines indicate the zones of the surface area of the specimen exposed to the alkaline solution. At these zones, the concentration of primary aqueous species in pore solution was imposed to be the same as in the alkaline bath. These concentrations were estimated by performing preliminary speciation calculation as it is described in Section 3.5.4, resulting in c s ¼ 0; c c ¼ 0:034 mol/m 3 , and c r ¼ 1000 mol/m 3 . Depending on the type of lateral sealing (Fig. 12), the zones exposed to the alkaline bath change from one set of specimens to the other. For simulating the specimens with complete lateral sealing (Set #2), concentrations were imposed only at zone A; for specimens without any lateral sealing (Set #1), at zones A, B, and C; and for the specimens of Set #3, at zones A and B.
The concentration of aqueous silica in the pore solution of mature HCP is usually very low (< 0:1 mol/m 3 ) (Berner, 1992;Lothenbach et al., 2007), and, it is not expected to play a significant role in the development of ASR. Therefore, the initial concentration of silicate ions in the pore solution of the HCP was assumed c s ¼ 0. The initial concentration of aqueous alkalis was fixed in c r ¼ 600 mol/m 3 by trial and error in order to best fit the experimental results. The initial concentration of calcium ions was assumed to be the equilibrium calcium concentration of portlandite in a 0.6 M NaOH solution, estimated to be c c ¼ 0:08354 mol/m 3 by means of a speciation calculation (Section 3.5.4).
The effective diffusivity (D b ) and the total porosity (/) considered for the continuum elements are given in Table 2. Note that the SL glass was assumed to be non-porous and impervious (Section 2.3). The capillary and total porosity of the HCP was obtained by estimating the amount of cement hydration products using the expressions given by Brouwers Brouwers (2004, 2005, assuming a maturity factor of 1.00. With the resulting capillary porosity (/ cp ¼ 0:117), the effective diffusivities of the primary aqueous species (Table 2) were estimated using the expression proposed by Oh and Jang (2004).  The parameters adopted for the diffusivity law of the interface elements were: D b 1 =D b w ¼ 10; D b ITZ =D b w ¼ 7:64 Â 10 À3 ; w ITZ ¼ 2 Â 10 À5 m, and w crit ¼ 1 Â 10 À4 m. The diffusivities of the primary species in bulk water at 60°C are given in Table 3. For the aqueous alkali, the diffusivity of Na + was adopted. The initial volume of the interface elements (U o mp ) and the corresponding initial surface area of dissolution fronts (S o mp ) are given in Table 1. The specific molar volume (g a o ), the unstressed gel porosity (/ a o ), and threshold stress (r a th ) of the solid a-species (SL glass, portlandite, RPA, and RPB) are given in Table 4. The dimensionless shape parameters of the f a / functions (Eq. (55)) were the same for RPA and RPB: k ¼ 1 and p ¼ 2. The compliance constant was the same for the four solid species: k g ¼ 1 Â 10 À6 MPa À1 .
The saturation product constants (K a sp ), and the formation (k a f ) and dissolution (k a d ) kinetic constants of the solid species are given in Table 4. K S sp and K C sp were estimated based on data from a geochemical database (Blanc et al., 2012) as it is described in Liaudat (2018). The saturation index of SL glass was well under 1 in all the simulations, and, consequently, k S f had no effect on the results. The remaining saturation product constants and kinetic constants in Table 4 were adopted to fit the experimental results. Finally, the value taken for the exponent of the activity of the hydroxyl ions in the kinetic law of SL glass dissolution was p ¼ 3, and for the fitting parameter of the reduction factor f S 2 was A S ¼ 13 m 2 /mol.

Expansion curves
The interfacial expansion curves obtained with the model for the three different lateral sealing cases are plotted together with the experimental curves in Fig. 13. As it can be appreciated in Fig. 13, the model seems capable of reproducing the effect of the lateral sealing on the experimental expansion rates.

Diffusion of aqueous species
Figs. 15a-c show the concentration of aqueous silicate, calcium and alkalis in pore solution at day 20 for the three boundary conditions considered. In the three cases, a diffusion front of alkalis moved into the specimen from the boundaries exposed to the alkaline solution (Fig. 15c). In Set #2, at day 20, the diffusion front of alkalis was still far from the ITZ. Moreover, the alkali concentration along the ITZ did not change significantly from its initial value. This seems to indicate that the amount of alkalis consumed in the formation of reaction products was approximately the same as the amount of alkalis released in the dissolution of SL glass. In contrast, for Sets #1 and #3, the alkali concentration along the ITZ was not uniform, with higher values close to the end exposed to the alkaline solution.
The calcium concentration (Fig. 15b) shows an inverse diffusion pattern, with calcium diffusing from the HCP towards the borders exposed to the alkaline bath. It must be noted that, although the formation of reaction products consumed significant amounts of calcium, the high rate of dissolution of portlandite in the interface elements kept the calcium concentration close to the saturation concentration (w C % 1). Consequently, the concentration of calcium at the interface elements was mainly determined by the variations Fig. 15. Modelling of the Interfacial Expansion Tests. Concentration of (a) silicate, (b) calcium, and (c) alkali ions in pore solution at day 20.

Table 3
Diffusion coefficients in [10 À9 m 2 /s] for aqueous species in bulk water at 25 and 60°C. Diffusivities in bulk water at 25°C are from (Lide, 2004, p. 5-95) with the exception of the diffusivity of silicate ions which is from Yokoyama (2013). Diffusivities at 60°C have been extrapolated using Stoke-Einstein relation (Yuan-Hui and Gregory, 1974) with water viscosities given in Kestin et al. (1981). of the pH associated with the diffusion of alkalis, rather than by the diffusion of calcium itself. The concentration patterns of silicates in pore solution are given in Fig. 15a. Since the initial concentration in pore solution and the imposed concentration at the exposed borders were zero, the only source of the silicates in pore solution was the SL glass dissolution. At the beginning of the simulations, the rate of SL glass dissolution was very high and the amount of silicate ions produced was higher than the amount consumed in the formation of RPA and RPB. Consequently, a small amount of surplus silicates diffused towards the bulk HCP. Rapidly, the rate of SL glass dissolution was reduced and the silicate ions that had diffused towards the HCP, diffused back to the interface to be consumed in the production of RPA and RPB. In any case, the amount of silicate ions that diffuse out or into the interface elements was a very minor part of the silicate ions produced by the SL glass dissolution which were in their greater part immediately consumed in the formation of RPA and RPB.

Profiles of reaction products and internal pressure
The effect of the different boundary conditions can also be seen in the radial profiles of reactions products at the HCP-SL glass interface (Fig. 16a). In Set #2, which had uniform profiles of concentrations in pore solution in the radial direction (Fig. 15), the volume of reaction products (U RPA mp þ U RPB mp ) along the HCP-SL glass interface was uniform at all times. In contrast, the volume of reaction products at the HCP-SL glass interface in Sets #1 and #3 was higher close to the right face exposed to the alkaline solution. This occurred more markedly for Set #1, but it can also be appreciated for Set #3. Moreover, in Set #1 the boundary condition prescribed of c s ¼ 0 at the nodes of the interface elements located at the right boundary implied that no reaction products could be formed there, since silicate anions are needed to form RPA and RPB (Section 2.4).
The profiles of induced internal pressure are given in Fig. 16b. Note that the time (day 5) is different from the one in the previous figures (day 20). This is because at day 20 the cracking process of the interface element was already very advanced (W cr =G i f % 1) and, therefore, the internal pressure was practically zero. The profiles in Fig. 16b reflect the non-uniformity of the formation of reaction products along the interface in Sets #1 and #3. Despite the lack of uniformity of r g , the profiles of relative normal displacement (a n ) and normal stress (r n ) (not presented here) were sensibly uniform. This is explained by the fact that the continuum elements above and below the interface elements behaved practically as two rigid bodies in relation to the compressibility of the reaction products.

Chemical and mechanical evolution of the interface
Let us now consider the time evolution of chemical and mechanical variables at a single integration point within an interface element of the specimen totally sealed (Set #2). Note that all the integration points of the interface elements in Set #2 have the same evolution curves, since the concentration profiles of the primary aqueous species remained uniform in the radial direction (Fig. 15a).
The time evolution of the three primary aqueous species is shown in Fig. 17a-c. These concentrations locally determine the pH of the pore solution (Fig. 17d) and the saturation indices of the four solid species (Fig. 17g). The evolution of the amounts (moles) of dissolved (SL glass, portlandite) and formed (RPA, RPB) solid species is shown in Fig. 17h. The Ca/Si and R/Si molar ratios of the mixed reaction product (RPA and RPB) are given in Fig. 17f. The volume occupied or released by each solid species is shown as an area graph in Fig. 18a, as well as the corresponding evolution of the total porosity in Fig. 17e. In Fig. 18b, the evolution of the normal stress (r n ), the internal pressure exerted by the reaction products (r g ), and the work spent on fracture processes (W cr =G f i ) are plotted together. At the beginning of the test, the dissolution rate of SL is very high and, consequently, there is a sharp rise of the concentration of silicate ions (Fig. 17a). These sudden peak made the saturation indeces of RPA and RPB (Fig. 17g) to pass from unsaturated (initially, c s ¼ 0 and, therefore, w A ¼ w B ¼ 0) to oversaturated condition (w A > 1; w B > 1), triggering the formation of reaction products. In turn, since the formation of RPA and RPB consumes silicate ions, the initial peak of c s is rapidly reduced (Fig. 17a).
After a small time lapse, the formed RPA and RPB occupied all the available space and an internal pressure r g started to be exerted (Fig. 18b). r g and the equilibrating r n grew rapidly until reaching the initial tensile strength (v o ¼ 1:5 MPa). At this point, the fracturing process started (a cr n > 0 and W cr > 0) and the resistance opposed by the interface to the expansion of the reaction products progressively decreased (softening behaviour), being practically null at the end of the test (Fig. 18b). This agrees well with the fact that, after the experimental tests, the interfaces of specimen could be detached with a slight pull of the ends.
The porosity in the interface had a maximum value at the beginning of the test, then rapidly decreased as RPA and RPB were formed reaching a minimum corresponding with the starting of cracking process. As the interface softens, the reaction products decompressed ( r g is reduced) and the total porosity grew back tending to the unstressed gel porosity of the mixed reaction product as r g tended to 0.
In Fig. 17h, it can be appreciated that the rate of SL glass dissolution steadily decreased with time, mainly as a result of the f S 2 factor intervening in the dissolution rate equation (Section 3.5.4). The M C mp and M A mp curves were roughly proportional to the M S mp curve. The M B mp curve had initially the same tendency as the M A mp curve, but, eventually, the saturation index of RPB became lower than 1, the reaction rate became negative, and part of the previously formed RPB was dissolved. The related evolution curves of the Ca/Si and R/Si molar ratios of the mixed reaction product (Fig. 17f) show similar values to those experimentally measured by means of SEM/EDS analysis on similar HCP-SL glass specimens (Liaudat, 2018;Liaudat et al., 2019).
Note that the pH of the pore solution is about 12.6 ( Fig. 17d). With this pH, the modelling assumption that the only aqueous silica species in pore solution is H 3 SiO 4 À might not be totally accurate, and a more realistic prediction of the pH of the pore solution could be obtained by also considering the concentration of H 2 SiO 4 2À in the speciation analysis (Section 2.4). A more detailed discussion on the different chemo-mechanical interplays reflected in Figs. 17 and 18 can be found in Liaudat (2018).

On the role of water transport
The experimental interfacial expansion curves shown in Fig. 13 indicate that as the degree of lateral sealing was increased, the expansion rate decreased. This effect may be attributed to the restriction imposed to the access of alkali and/or water from the alkaline bath into the reaction site. However, since the present modelling framework does not consider the mass balance of water explicitly (Section 2.3), the second of those effects, i.e. the potential influence of (the lack of) water availability on the kinetics of the ASR expansions cannot be represented directly. This simplification does not seem to affect the ability of the model to reproduce the experimental expansion curves, as it can be appreciated in Fig. 13. Nonetheless, the fact that two of the fitted parameters had to be given values somewhat out of the expected range of values might be the consequence of that simplification. Indeed, the value of the initial concentration of alkalis (c r ¼ 600 mol/m 3 ) required to fit the experimental expansion curves turns out low for HCP with an alkali content of 1.1 wt% of Na 2 O e . According to experimental correlations in the literature (Lindgård et al., 2012;Diamond, 1989), the initial concentration of alkalis in pore solution is estimated to had been close to 1000 mol/m 3 . Clearly, further studies would be needed in order to elucidate the role played by water transport on the kinetics of the interfacial expansion curves.

Application example: single-aggregate specimen
The second application example consists of an idealized concrete specimen composed of a circular SL glass aggregate with a pre-existing crack, which is embedded in a square matrix of nonreactive mortar. For the non-reactive mortar the following parameters are assumed: water/cement weight ratio of 0.45, aggregate volume fraction of 0.43, and 1.25% of Na 2 O e per weight of cement. The cement is assumed to be CEM I 42,5 N-SR5 (AENOR, 2011), and the non-reactive sand is assumed to pass a 2.5 mm sieve and to be well graded. These parameters correspond to the mortar surrounding the SL glass particles in the SL glass concrete specimens used for the ASR Confined Expansion Tests reported by Liaudat et al. (2018). Moreover, whenever possible, the same model parameters calibrated for the above-presented simulations of Interfacial Expansion Tests have been also used in this case. The development of ASR expansions at 60°C is simulated under three different mechanical boundary conditions: (i) without applied loads (free expansion), (ii) compressive loading in the vertical direction, and (iii) compressive loading in both vertical and horizontal direction. These situations have been chosen because they are assumed to reproduce the most distinctive features of the ASR expansion mechanism proposed in Liaudat et al. (2019), and to do it with the simplest possible model geometry in order facilitate the interpretation of the results. In the three cases, plane stress condition is also assumed.

Model description
The model geometry is schematically represented in Fig. 19a. A stiff loading plate is placed on each side of the specimen in order to prescribe suitable boundary conditions for the mechanical analyses. The model was discretized by means of 1112 quadratic triangular elements and 268 zero-thickness interface elements (Fig. 19b), giving a total of 2876 nodes. The locations of the interface elements are indicated with thick lines in Fig. 19b. These interface elements are intended to simulate five different types of discontinuities: (1) cracks inside the SL glass particle (glass-glass interface), (2) cracks within the mortar matrix (mortar-mortar interface), (3) SL glass-mortar interface, (4) pre-existing cracks inside the SL glass particle, and (5) 'frictionless contacts' between the loading plates and the mortar specimen.
The interface elements representing the frictionless contacts were equipped with a linear elastic mechanical constitutive law, with high normal stiffness and very low tangential stiffness. For the rest of the interface elements, the mechanical elasto-plastic constitutive law for interface elements described in Section 4.2 was considered. Part of the parameters used for this constitutive law are given in Table 5 for each type of discontinuity. The rest  of the parameters were the same for all types: were assumed to be linear elastic, with elastic modulus E and Poisson's coefficient m given in Table 6.
The value of G i f used for the glass-glass interface elements was based on values reported in Linger and Holloway (1968) and Wiederhorn (1969) ranging between 2.0 and 3.8 J/m 2 . In order to facilitate the convergence of the numerical calculations, the brittleness of SL glass was reduced by adopting a rather low value of initial tensile strength for the glass-glass elements (short term tensile strength of float SL glass is about 45 MPa Schittich et al., 2007;Haldimann et al., 2008). The values E and m of the bulk SL glass were adopted according to EN 572-1:2004CEN (2004. The values of v o and G f i of the mortar-mortar interface elements and E of mortar were fixed on the basis of the experimental values reported in Liaudat et al. (2018) for the control concrete. Those values were corrected taking into consideration the different mixture proportions and the particle size distribution of the limestone aggregates.
The values of v o and G i f of the mortar-glass interface elements were assumed to be one half of their respective values adopted for the mortar-mortar interface elements, based on previous simulations of concrete at the meso-scale with similar mechanical models (López et al., 2008a;López et al., 2008b;Idiart et al., 2011). The parameters adopted for the interface elements representing the pre-existing cracks guarantee that they offer very low (practically zero) resistance to normal opening, very high stiffness upon closure, and a (practically) purely frictional behaviour under tangential relative displacements in compression.
For the diffusion problem, natural boundary conditions were assumed. The initial concentrations of primary aqueous species in pore solution were assumed to be the same as those adopted for the simulations of Interfacial Expansion Tests (Section 7).
The effective diffusivities (D b ) and the total porosity (/) considered for the continuum elements are given in Table 6. The SL glass and the loading plates were assumed to be non-porous and impervious. The total and gel porosity of the mortar was estimated on the basis of its mix composition, chemical composition of the cement, and particle size distribution of the non-reactive mortar aggregate. The effective diffusivities in the mortar were estimated with the procedure propose by Oh and Jang (2004). Further details on the steps followed for the estimation of mortar porosities and effective diffusivities can be found in Liaudat (2018). For the aqueous alkalis, the diffusivity of Na + was adopted.
For the interface elements representing frictionless constants the diffusivities were assumed constant, with small values in the longitudinal direction (D b l =D b w ¼ 1 Â 10 À10 ) and very high values in the transversal direction (D b n =D b w ¼ 1 Â 10 6 m À1 ). For the rest of the interface elements, the longitudinal and transversal diffusivities evolved according to the evolution law proposed in Section 3.6.2. In all cases, the parameters adopted were those calibrated for the Interfacial Expansion Tests (Section 7), with the exception of the width of the ITZ (w ITZ ) which was set null for the mortar-mortar and glass-glass interfaces. The diffusivities in bulk water (D b w ) at 60°C are given in Table 3. The initial volume of the interface elements (U o mp ) and the corresponding initial surface area of dissolution fronts (S o mp ) are given in Table 5. The values taken for the fitting parameter of the reduction factor f S 2 were A S ¼ 13 m 2 /mol for mortar-glass interfaces and A S ¼ 6:5 m 2 /mol for glass-glass interfaces.
The remaining parameters, namely g a o ; / a o ; k g ; r a th ; K a sp ; k a f ; k a d ; k g , the shape parameters k and p of functions f a / , and the exponent of the hydroxyl ions activity in the kinetic law of SL glass dissolution, were the same as those used for simulating Interface Expansion Tests (Section 7).
As mentioned above, three different boundary conditions were considered for the mechanical problem ( Fig. 20): (i) without applied loads (q x ¼ q y ¼ 0), (ii) compressive loading in the vertical direction (q x ¼ 0; q y ¼ 10 MPa), (iii) compressive loading in both vertical and horizontal direction (q x ¼ q y ¼ 10 MPa). Note that the combination of very stiff loading plates, frictionless contacts, n/a n/a n/a n/a and the imposed boundary conditions guarantees that the sides of the specimen do not bend nor rotate due to the applied loads and/ or the development of internal pressures. The time discretization of the staggered scheme consisted in time steps of 3 h. Each of these steps was further subdivided in 12 and 50 equal increments in DRACFLOW and DRAC, respectively. For the cases with external loading, the q x and/or q y were applied in one additional step at the beginning of the simulations, in which the time increments in DRACFLOW were null. In that way, the external loading of the specimens may be regarded as 'instantaneous'.

Free expansion
The first loading case considered is q x ¼ q y ¼ 0 (Fig. 20). The axial expansion curve obtained is given in Fig. 21 in terms of the homogenized strain u x =L, where u x [m] is the length change in the horizontal direction and L is the side length of the mortar square. Note that, in this case, the model has a symmetry axis in the direction of the pre-existing crack, and, consequently, the expansions in the horizontal and vertical directions were the same. The points marked on the expansion curve (Fig. 21) indicate significant times in the development of the ASR expansions, in which results will be depicted in subsequent figures, namely: For each of these moments, the deformed geometry, the spatial distribution of the reaction products (U A mp þ U B mp ), the Ca/Si molar ratio of the mixed reaction products, the internal pressure (r g ), and the work spent on fracture processes (W cr =G I f ) are given in Figs. 22-26. Fig. 22. Modelling of a single-aggregate specimen under free expansion conditions. Deformed geometry (Â500) at different times. Between the beginning of the simulation and day 0.875 (point A, Fig. 21), practically no expansion was developed. During this time, the reaction products formed in the mortar-glass interface elements progressively filled the available initial volume (U o mp ), eventually exerting a growing internal pressure (r g ) once the initial volume was exhausted (Figs. 23 and 25). Point A coincides with the moment in which r g reached a value slightly over v o and the cracking process started (Fig. 25). Simultaneously, reaction products were also formed within the pre-existing crack but without exhausting its initial volume, i.e. without exerting any internal pressure. The profiles of Ca/Si molar ratio indicate that the reaction products at the mortar-glass interface were mainly RPA and that the reaction products within the pre-existing crack were mainly RPB (Fig. 24).
After point A, reaction products continued to precipitate both at the mortar-glass interface and within the pre-existing crack (Fig. 23). As the normal cracking aperture (a cr n ) of the mortarglass interface elements grew, the reactive stress exerted by the  mortar matrix increased and the local rate of formation of reaction products decreased due to the reduction of the porosity (/ mp ). Eventually, the formation of reaction products at the mortarglass interface elements stopped altogether at day 5.5, when r g became equal to the threshold stress of the mixed reaction product (r th ). Meanwhile, the formation of reaction products within the pre-existing crack continued (Figs. 23 and 24), exhausting the initial volume at day 2, and from there on exerting a growing internal pressure that eventually propagated the pre-existing crack both into the glass aggregate and into the mortar (Fig. 25). The crack propagation into the mortar was initiated earlier (point B), because the previous formation of reaction products at the mortar-glass  interface had tensioned the mortar but compressed the glass. The propagation into the SL glass started five days later, at day 13.75 (point C). The resulting cracking apertures enabled the dissolution of portlandite in mortar-mortar interfaces, the dissolution of SL glass in glass-glass interfaces, and the formation of RPA and RPB in both types of interfaces. Note that if no pre-existing crack was considered in the glass aggregate the ASR would have completely stopped when r g at the glass-cement interface had reached the threshold stress of the reaction products. That means that no cracks would propagate into the mortar matrix nor into the glass particle, and that no expansion would develop beyond that point. This result is qualitatively similar to the experimental results obtained by Maraghechi et al. (2016) and Maraghechi et al. (2012) using crack-free SL glass aggregates in Portland cement mortars. For a sensitivity study of the effect of the length and orientation of the pre-existing crack on the overall expansions of the specimen, see Liaudat (2018). It can be seen in Fig. 25 that from day 13.75 (point C) until the end of the simulation, there was a progressive increase of r g at the mortar-glass interface elements at both sides of the pre-existing crack, reaching values over r th . This increase of r g occurred as a consequence of the compression of the previously formed reaction products (da cr n =dt < 0), and not to the formation of additional amounts. The reduction of a cr n at the mortar-glass interface was the result of the propagation of the pre-existing crack in the framework of the kinematic restrictions imposed by the loading plates. Because of the same reasons, the normal stress (r n ) of the mortar-mortar interfaces at the left-bottom quadrant of the specimen progressively increased as the pre-existing crack propagated, eventually reaching v o at day 30.75 (point D), giving rise to three new cracks growing radially from the SL glass aggregate (Fig. 26).
From day 30.75 until the end of the simulation, the propagation of the cracks previously initiated continued, in some cases reaching the borders of the specimen (points E and F). The sudden rise of the expansion rate after point E (Fig. 21) indicates the configuration of a failure mechanism (separation of the right-bottom part of the specimen) that eventually led to the end of the simulation due to lack of convergence of the numerical solution. Note that in a more realistic concrete specimen with multiple coarse aggregate particles this kind of 'brittle' failure mechanism would not have occurred, since the crack propagation would have had to go over more tortuous paths before configuring a failure mechanism. Therefore, the development of ASR strain curves after point E cannot be considered as representative of the behaviour of real SL glass concrete.
It is remarkable that the spatial distribution of Ca/Si molar ratios of the mixed ASR products formed within the specimen (Fig. 24) is qualitatively similar to the one observed experimentally (e.g. Rajabipour et al., 2010;Maraghechi et al., 2012;Liaudat, 2018), i.e. low-calcium products formed at locations in contact with the cement paste, and high-calcium products within glass cracks.

Expansion under uniaxial and biaxial compression
Let us now consider the uniaxial compression (q x ¼ 0; q y ¼ 10 MPa) and the biaxial compression (q x ¼ q y ¼ 10 MPa) cases, as depicted in Fig. 20.
The axial strain curves obtained are given in Fig. 27, together with the 'free expansion' curves (Section 8.2.1). The initial strains generated by the application of q y and/or q x were deducted from the strain curves, in order to facilitate the comparison of the ASR-induced strains obtained under different loading cases. The corresponding volumetric strain curves, considered as the sum of the horizontal and vertical strains, are given in Fig. 28. The deformed geometry and the spatial distributions of the reaction products (U A mp þ U B mp ), the internal pressure (r g ), and the work spent on fracture processes (W cr mp =G I f ), are given for one modelling time (day 34) in Figs. 29 and 30. The time chosen (day 34), corresponds to the moment just before a crack reached the upper border of the specimen in the uniaxial loading case.
Under uniaxial (vertical) loading, the specimen developed a progressive expansion in the lateral (horizontal) direction and an associated contraction in the vertical direction ( Fig. 27a and b). The rate of expansion in the horizontal direction was higher than in the free expansion case. However, the increase in the rate of hor- Fig. 27. Modelling of a single-aggregate specimen under horizontal (q x ) and vertical (q y ) uniform external pressure. Axial strains in the (a) horizontal and (b) vertical directions. The initial strains generated by the application of q y and/or q x were deducted from the strain curves. izontal expansion was not sufficient to compensate the contraction in the vertical direction, and, therefore, the volumetric expansion rate of the uniaxial case before the cracks reached the specimen border (day 34), was lower than in the free expansion case (Fig. 28). As in the free expansion case, the sudden rise of the expansion rate in the uniaxial case after day 34 was due to the configuration of a failure mechanism (separation of the right side part of the specimen) when cracks reached the borders of the specimen.
The effect of uniaxial loading on the development of the ASR expansions can be also appreciated in the cracking patterns. While in the free expansion case the cracking pattern was symmetric with the pre-existing crack propagating in the diagonal direction (Figs. 22 and 26), in the case of uniaxial loading the propagation of cracks was vertically oriented (Figs. 29a and d).
In the biaxial case, the applied loads practically inhibited the ASR expansions (Figs. 27a and b). The inhibition of the ASR expan-  sions is confirmed by the absence of cracks in the specimen after 34 days (Figs. 30a and d).
Same as for the Interfacial Expansion Tests, some time (0.33 days) was needed by the ASRs to fill the initial volume (U o mp ) of the mortar-glass interfaces, before r g started to be exerted and strains were developed (Figs. 27). From day 0.33 onwards, r g at the mortar-glass elements started to grow, inducing strains in the mortar specimen. In the free expansion and biaxial load cases, the expansion curves remained identical until day 0.875, when the expansion rate of the free expansion case suddenly rose as the mortar-glass elements developed cracks. In the case under biaxial loading, the isotropic compression, which was well over the threshold stress of the calcium-rich reaction product, prevented the cracking of the mortar-glass elements and, consequently, the rate of formation of reaction products asymptotically tended to zero as r g ! r th . In the case under uniaxial loading, the cracking of mortar-glass elements occurred slightly earlier than the free expansion case, with different consequences for the vertical and horizontal strains. In contrast to the free expansion case, in which the cracking of the mortar-glass elements occurred in pure tensile mode (mode I), under uniaxial loading the cracking occurred in mixed shear-compression mode, due to the initial stress state around the SL glass aggregate. This initial shear-compression stress state was determined by the higher stiffness of the glass aggregate in relation to the mortar matrix. Under those initial conditions, the growing r g at the mortar-glass elements had the effect of progressively reducing the effective compressive normal stress (r n ) and therefore also the shear strength determined by it (Section 4.2).
Eventually, if r g was sufficiently high, the hyperbolic failure surface was reached and the cracking process started (Section 4.2). As the mortar-glass elements cracked, the tangential relative displacement (a cr l ) associated to the mixed cracking mode allowed the mortar to 'slip' around the aggregate under the effect of the vertical load, with subsequent descent of the top-side loading plate (Fig. 27b). It can be expected that in a more realistic concrete specimen with multiple coarse aggregate particles this mechanism would be less significant due to the additional internal kinematic restrictions.
During the first five days, the horizontal expansion curve of the uniaxial load case was very similar to the horizontal expansion curve of the free expansion case. But approximately at day 5, the cracks in the uniaxial loading case started to propagate vertically both into the mortar matrix and inside the aggregate and, consequently, the horizontal expansion curve grew faster than in the free expansion case.
As cracks propagated vertically, the passive restraint exerted by the mortar on the expansion of the mortar-glass interfaces in the horizontal direction, was reduced, allowing further interfacial expansions with the subsequent degradation of the shear strength. This reduction of the shear strength, as explained above, caused further descent of the top-side plate under the effect of the vertical load.
It is remarkable that the effect of loading on the development of the ASR strains in the single-aggregate model is qualitatively similar to the effects observed in the experimental Confined Expansion Tests with SL glass concrete specimens reported in Liaudat et al. (2018), namely, (i) reduction of the volumetric expansion rate for increasing volumetric stress, and (ii) increase of the expansion rate in one direction when compressive stress is applied perpendicularly to this direction. It might be argued that the concrete mesostructure is intrinsically 3D and therefore that the 2D analysis presented in this section might not be representative of the real concrete behaviour. However, based on their experience in previous 2D/3D numerical studies performed with similar models (purely mechanical Caballero et al., 2006;Caballero et al., 2007;López et al., 2008a;López et al., 2008b andchemo-mechanical Idiart et al., 2011;Pérez et al., 2017), the authors do not expect that 3D simulations with the proposed ASR model would lead to significant differences in the qualitative results. Nonetheless, this presumption clearly need to be confirmed with further research.

Concluding remarks
A coupled C-M FE model for simulating ASR in SL glass concrete at the meso-scale has been formulated theoretically and implemented numerically. The development of the model has led to the following numerical achievements: A new FE formulation for diffusion-reaction processes within discontinuities represented via double-node zero-thickness interface elements has been proposed. The formulation is capable of reproducing both dissolution/precipitation of solid species occurring in the water filling the discontinuity as well as dissolution fronts occurring at both sides of the discontinuity. This diffusion-reaction formulation is coupled to the mechanical formulation through the normal relative displacement variable, in this way making it possible to introduce the effect of crack opening both in the evolution of the transport parameters and in the kinetics of the chemical reactions. A new constitutive law for interface elements has been proposed, which is able to reproduce the propagation of a crack induced by the development of an internal pressure exerted by solid reaction products formed within it. The numerical implementation of the diffusion-reaction formulation is done in analogy with traditional mechanical implementations. This is carried out by assuming that the concentrations of the primary aqueous species in pore solution (analogous to the displacements in the mechanical problem) at each integration point vary linearly within each time increment. Then, this linear variation of the concentration is introduced in separated routines called 'chemical constitutive law' and 'diffusivity law' (analogous to the mechanical constitutive laws) where the evolution of the solid species fractions and effective diffusion coefficients within the time increment are calculated.
The model considers several C-M coupling mechanisms. In particular, a reaction-expansion mechanism, qualitatively proposed by the authors elsewhere (Liaudat et al., 2019) on the basis of experimental results, is introduced in order to reproduce the effect of compressive stresses on the development of ASR expansions. It must be noted that SL glass has some peculiarities as concrete aggregate, and, thus, the model includes some assumptions that may be not directly applicable to the more usual case of concrete made with natural reactive aggregates. However, the changes needed in the formulation in order to include such type of reactive aggregates, should be relatively minor and, therefore, adaptation should be relatively straight forward. Moreover, although the model is presented and developed in 2D, the conclusions reached and the formulation developed should be valid and readily extensible to 3D.
The model has been used for analysing two application examples under different boundary conditions. In the first example, the model has been able to reproduce the experimental results of Interfacial Expansion Tests, both regarding the development of ASR expansions and regarding the chemical composition of the reaction products. Remarkably, the model has been able to reproduce, using a unique set of model parameters, the effect of three different lateral sealing schemes on the experimental expansion rates. In the second example, using the same model parameters calibrated in the first example, the effects of external loading on the development of ASR expansions and cracking in SL glass concrete have been qualitatively reproduced with the simulation of ideal 2D single-aggregate specimens. These effects are the reorientation of cracks, the reduction of the volumetric expansion rate for increasing volumetric stress, and the increase of the expansion rate in one direction when a compressive stress is applied in the perpendicular direction. Additional application examples with a single aggregate can be found in Liaudat (2018). Single-aggregate examples such as the ones presented in this paper are necessary steps in the process of validating the proposed model. Validation is still an ongoing task which must include additional more complex cases such as multiple SL glass particles and 3D analysis. However, the results obtained so far seem to indicate that the reactionexpansion mechanism proposed in Liaudat et al. (2019) and implemented in this model is capable of explaining the main aspects of the C-M mechanisms involved in the ASR expansions of SL glass concrete.
Moreover, the C-M formulation for zero-thickness interface elements with double nodes proposed in this paper has demonstrated its capabilities to simulate diffusion-reaction processes in discontinuities combined with mechanical actions. This formulation could be used in the future to model other C-M expansive processes in interfaces or discontinuities. Additionally, it is remarkable that the simplifying assumption of considering that chemical reaction can only occur in the interface elements, representing ITZ and open cracks, significantly reduces the computational cost without decreasing the ability of the model to reproduce the experimental behaviour.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.