A new theoretical phase eld-based formulation for predicting electro- chemo-mechanical corrosion in metals is presented. The model com- bines electrolyte and interface electrochemical behaviour with a phase eld description of mechanically-assisted corrosion accounting for lm rupture, dissolution and repassivation. The theoretical framework is numerically implemented in the nite element package COMSOL MUL- TIPHYSICS and the resulting model is made freely available. Several numerical experiments are conducted showing that the corrosion predic- tions by the model naturally capture the in uence of varying electro- static potential and electrolyte concentrations, as well as predicting the sensitivity to the pit geometry and the strength of the passivation lm. Keywords: Phase eld, localized corrosion, mechanically-assisted corrosion, electro-chemo-mechanical modeling, nite element method, COMSOL 1 Introduction Corrosion is widely recognized as one of the most common yet destructive fail- ure mechanisms of engineering components and structures, carrying a cost of 1 2 about 3.1% of the gross domestic product (GDP) in industrialized countries [1]. Corrosion can be classi ed into two categories: general corrosion and local- ized corrosion. The rst one refers to uniform metal dissolution and can be prevented by physical and electrochemical measures, e.g., through the use of coatings and cathodic protection. However, localized corrosion, such as pitting corrosion, is more dicult to prevent and detect in engineering practice. More- over, by interacting with mechanical loads, localized corrosion can often lead to catastrophic structural failures through phenomena such as mechanically- assisted corrosion and stress corrosion cracking (SCC). Thus, there is a need for developing physically-based models capable of predicting localized corrosion failures. Mechanistic predictions of localized corrosion require resolving the under- lying physical processes, which has long been considered a remarkably complex task [2]. Two main challenges exist. First, predicting localized corrosion requires resolving a strongly coupled electro-chemo-mechanical problem involv- ing, at the very least, the transport of ionic species, interface reactions, changes in electrolyte conductivity and mechanical behaviour of the metal, as well as their interactions. Secondly, localized corrosion is a complex interfacial prob- lem; the electrolyte-metal interface exhibits a highly irregular morphology and its evolution depends on the local chemistry and mechanics, which are themselves dependent on the interface morphology [3]. Promising progress has been achieved, independently, in the modeling of these two key challenges. New mathematical models describing mass transport, homogeneous chemical reactions and charge balance have been proposed to better understand the electro-di usive mass transport of multiple species and its e ects on corrosion evolution (see, e.g., [4{6] and Refs. therein). However, these works are limited by the assumption of a stationary electrolyte-metal interface. On the other hand, a number of computational schemes have been recently presented to track the evolution of the corrosion front [7{14]. Among these, the phase eld model stands out for its thermodynamics roots and its ability to eciently capture corrosion morphologies of arbitrary complexity. In addition, it can be readily integrated into existing nite element schemes. Instead of explic- itly tracking a moving boundary (and the associated boundary conditions), the phase eld paradigm describes the evolution of the interface by means of an auxiliary order parameter that varies smoothly between the two phases (metal and electrolyte, in the case of corrosion) and evolves based on a suitable governing equation [15{17]. The application of the phase eld paradigm to cor- rosion was rst shown by Stahle and Hansen [18] and Abubakar et al. [19]. Mai and co-workers extended these works to present a phase eld corrosion frame- work that accounts for the di usion of dissolved species and reactants in the electrolyte [8, 20]. Later, Cui et al. [12] incorporated the role of mechanics in the corrosion process, capturing both the enhancement of corrosion rates due to mechanical elds [21] and the process of passive lm formation, rupture and subsequent repassivation [22]. More recently, mechanics-enhanced phase eld corrosion models have been extended to include mechanical fracture [23, 24], 3 which could be assisted by the ingress of aggressive species such as hydrogen [24]. However, these models do not fully resolve the electrochemistry of the system; there is a need to develop electro-chemo-mechanical phase eld formu- lations for localized corrosion that can capture the sensitivity to the applied potential and the concentration of species in the environment. In this work, we present a new electro-chemo-mechanical phase eld-based formulation for predicting localized corrosion in elastic-plastic solids. The model combines, for the rst time, an electrochemical description of ionic species transport and electrostatic potential distribution in the electrolyte with a mechanics-dependent interface kinetics law built upon mechanochemical theory [21] and the lm rupture{dissolution{repassivation (FRDR) mech- anism [22]. The theoretical framework is numerically implemented in the nite element package COMSOL MULTIPHYSICS and numerical experiments are conducted to showcase the ability of the model in capturing the main experimental trends. Focus is also placed on the details of the COMSOL imple- mentation, making the model freely available and providing details that enable reproducibility and facilitate usage. The remainder of this paper is organ- ised as follows. In Section 2 we present our theoretical framework. Then, the COMSOL implementation is described in Section 3. The results of the numer- ical experiments conducted are presented and discussed in Section 4. Finally, concluding remarks end the manuscript in Section 5. 2 Theory 2.1 General considerations The process of localized corrosion is usually initiated by the rupture of the passivation lm that protects the metal from corroding. Following the local failure of the passive lm, the metallic surface is exposed to the corrosive environment, triggering the dissolution process and releasing cations into the electrolyte, i.e. M!Mn+ M+nMe(1) where M is the corroded metal and nMis the charge number of metal M. The present theory aims at encapsulating the electrochemical and mechanical mechanisms involved in this process. An overview of the elements of our theory is provided in Figure 1. There, it can be seen that the evolution of the corrosion front is described by a so-called phase eld order parameter , which varies from 0 (electrolyte) to 1 (intact metal) within a di use region. The interplay between activation-controlled corrosion and di usion-controlled corrosion is captured by simulating the transport of metal ions. To this end, a normalized concentration cM=c0 M=csolid is de ned, where c0 Mis the metal ion concentration and csolidis the concentra- tion of atoms in the metal. Accordingly, cM= 1 inside the metal while cM= 0 in regions of the electrolyte that are far away from the electrolyte-electrode interface. Following our previous work [12], we characterize the mechanical behaviour of the metal using an elastic-plastic constitutive model and capture 4 Fig. 1 Schematic illustration of the main variables involved in the electro-chemo-mechanical phase eld-based formulation presented for predicting localized corrosion as a function of the environment, material and loading conditions. The illustration showcases a domain that encompasses the metal ( = 1) and electrolyte ( = 0) phases, and the di use interface in-between. as well the interplay between mechanics and corrosion (mechanochemistry, FRDR mechanism). Unlike our previous work [12], the current model also resolves the electrochemistry of the electrolyte. This is achieved by solv- ing for the concentration of multiple ionic species ci(whereiequals, e.g., Na+;Cl;H+;OH, etc.), and the electrolyte potential 'l. Accordingly, the primary variables of the theory are the phase eld order parameter , the normalized metal ion concentration cM, the electrostatic potential 'l, the concentration of multiple ionic species ciand the metal displacement vec- toru. Based on these preliminaries, we shall rst de ne the energy functions (Section 2.2), then derive the balance equations of the primal elds (Section 2.3), introduce the main corrosion-mechanics coupling relations in the electro- chemo-mechanical system (Section 2.4), particularize the electrochemistry of the electrolyte (Section 2.5) and, nally, summarize the governing equations of our theory (Section 2.6). 2.2 Free energy Consistent with our prior research [12], we consider two distinct systems, namely the electrochemical system and the mechanical system. These two sys- tems possess their respective free energies, denoted by Eand M, which will be formulated in Sections 2.2.1 and 2.2.2, respectively. The coupling of these two systems will be addressed in Section 2.4. 5 2.2.1 Electrochemical system In this section, we introduce a general expression for the energy function in corrosion problems. Consider rst the free energy of the electrochemical sys- tem, denoted by E, which is an integral of the electrochemical energy density Eover the domain , such that E=Z EdV=Z  ch+ i+ el dV (2) where ch; i, and eldenote the chemical, interfacial, and electrostatic free energy densities, respectively. These are de ned below. Chemical free energy density ch. The chemical free energy density ch can be further decomposed into the energy for metal dissolution ch;dand the energy stored in the dilute solution ch;s, ch= ch;d+ ch;s(3) We follow the KKS model [25] to de ne ch;d. Accordingly, each material point is a mixture of both solid and liquid phases with di erent concentrations but equal chemical potentials. Thus, ch;dis given by ch;d=h( ) ch;d S+ [1h( )] ch;d L (4) whereh( ) =2 3+3 2is a degradation function intrinsic to the phase eld model. Here we de ne h( = 1) = 1 as intact metal and h( = 0) = 0 as fully corroded region to describe the phase transition due to metal dissolution. Also,h0( = 0) =h0( = 1) = 0 must be satis ed to ensure that the energy converges to a nite value where the domain is locally intact/fully damaged. ch;d Sand ch;d Lare the chemical free energy density terms respectively asso- ciated with the concentrations of the solid phase cSand the liquid phase cL. Accordingly, @ ch;d S(cS) @cS=@ ch;d L(cL) @cL(5) and, cM=h( )cS+ [1h( )]cL (6) Consistent with the KKS model, ch;d Sand ch;d Lare given by ch;d S=A(cScSe)2and ch;d L=A(cLcLe)2(7) wherecSe=csolid=csolid= 1 andcLe=csat=csolidare the normalized equilib- rium concentrations for the solid and liquid phases. Also, Ais the free energy density parameter, which is assumed to be equal for the solid and liquid phases and can be determined by benchmarking the chemical free energies obtained from KKS model with those obtained from thermodynamic databases [19, 23]. 6 Combining Eqs. (4)-(7), one reaches ch;d=A[cMh( )(cSecLe)cLe]2(8) The second element of ch, the energy stored in the dilute solution, ch;s, is de ned as follows: ch;s=nX i=1ciRgT(lnci1) +ci0 i (9) whereiis the associated ionic species (e.g., Na+;Cl;H+;OH),Rgis the gas constant, Tis the temperature, and 0 iis the reference chemical potential. Here, one should note that the chemical energy associated with the metallic ions Mn+ Mis already accounted for through the KKS model. Interfacial free energy density i. The interfacial strain energy density iis de ned as the sum of double-well potential energy and the energy corresponding to the phase eld gradient, reading i=g( )w+ 2jr j2(10) where is the gradient energy coecient and wis the height of the double-well potentialg( ), which is chosen here as g( ) = 2(1 )2. As discussed and derived in Appendix A, the height wand the gradient energy coecient in Eq. (10) can be related to the interface energy per area and its thickness `as: =r w 18and`=r 8 w(11) Electrostatic free energy density el. Finally, the electrostatic energy density elis de ned as a function of the charge density as follows el=F'l nMcMcsolid+nX i=1nici! (12) whereFis Faraday's constant, and niis the charge number of the ithionic species. 2.2.2 Mechanical system We proceed to de ne the mechanical strain energy M. We consider an elastic-plastic solid with a strain energy density Mthat can be additively 7 decomposed into its elastic ( e) and plastic ( p) parts, such that M=Z MdV=Z h( ) ( e+ p) dV (13) where the corrosion degradation function h( ) captures the reduction of mate- rial sti ness due to metal dissolution. The elastic strain energy density is given as e=1 2"e:Ce:"e(14) whereCeis the linear elastic sti ness matrix and "eis the elastic strain tensor. The plastic behaviour of the solid is characterized in this work by means of von Mises J2 plasticity theory. Accordingly, for a plastic strain tensor "p, the plastic component of the strain energy density reads, p=Zt 00:_"pdt (15) where0is the so-called undamaged or e ective Cauchy stress tensor. 2.3 Balance equations Consistent with the free energy de ned in Section 2.2, we shall now formulate the governing equations of the electro-chemo-mechanical phase eld corrosion theory. 2.3.1 Phase eld evolution As in previous phase eld corrosion models [8, 12], the phase eld evolves based on the following Allen{Cahn equation: @ @t=L@ E @  r2  (16) with @ E @ =2A[cMh( ) (cSecLe)cLe] (cSecLe)h0( ) +wg0( ) (17) Here,Lis the interface kinetics coecient, which can be related to the corro- sion current density [24, 26] or the overpotential [27]. Following Ref. [12], we shall enrich the de nition of Lto incorporate the role of mechanics in enhanc- ing corrosion rates and breaking the passivation lm - see Section 2.4. For the moment, let us assume that Lis constant and denote LLafor the case where mechanical elds are absent. Then, as discussed in Ref. [24], Lais proportional 8 to the corrosion current density iawhen corrosion is activation-controlled. Accordingly, =La ia=L0 i0(18) wherei0is the exchange current density and L0is the interface kinetics coecient when the overpotential is zero. Thus, as shown in Section 4.2, the proportionality constant can be determined by conducting a numerical experiment under activation corrosion conditions (small La). The corrosion current density iacan usually be estimated by using the so-called Butler{Volmer equation: ia=i0 expaanMF RgT exp (1aa)nMF RgT (19) whereaais the anodic charge transfer coecient. Consistent with Eqs. (18)- (19),Lacan also be estimated by a Butler{Volmer-type equation: La=L0 expaanMF RgT exp (1aa)nMF RgT (20) And Eq. (20) can be further simpli ed as a Tafel-type equation if only the anodic reaction is relevant at the dissolution surface; i.e., La=L0expaanMF RgT (21) Finally, the overpotential is expressed as: ='s'lEeq (22) where'sis the solid (applied) potential and Eeqis the equilibrium potential. In this work, we assume Eeq= 0 for convenience. 2.3.2 Transport of species The transport of ionic species is governed by mass conservation. Accordingly, the rate of change in time of any of the species must be equal to the sum of its concentration ux through the boundary @ and the reactants/products due to chemical reactions in , such that Z @c @tdV=Z @ JndS+Z RdV (23) where Jis the concentration ux and Ris the chemical reaction term. Here, one should note that Eq. (23) is valid for both c0 Mandci. Since Eq. (23) must 9 hold for any arbitrary volume, recalling that c0 M=cMcsolidand using Gauss' divergence theorem one reaches @cMcsolid @t+rJM=RM (24) and@ci @t+rJi=Ri (25) The transport process is driven by the chemical potential . For the metal ions,Mcan be further decomposed into two terms, one associated with the KKS-based phase eld formulation, M1, and another one related to the migration process, M2, such that M=@ E @cMcsolid= 2A[cMh( ) (cSecLe)cLe]=csolid| {z } M1+F'lnM|{z} M2(26) Where we emphasise that csolidis a constant used for normalising the metal ion concentration. Accordingly, the ux JMcan be calculated by a Fick law-type relation, JM=JM1+JM2=DM 2Ac2 solidrM1[1h( )]DMcMcsolid RgTrM2 =DMcsolidr[cMh( ) (cSecLe)cLe][1h( )]DMcMcsolid RgTFnMr'l (27) Note that the term [1 h( )] is present in JM2to ensure the transport of species is only valid in the electrolyte and along the interface ( <1). However, this term is not necessary for JM1, given that cMand are naturally coupled by the KKS model. Also, note that we use the real metal ion concentration c0 M=cMcsolidin (27) to maintain the dimensional consistency of JMandJi. Now, inserting (27) into the mass conservation equation (24), the transport of metal ions is formulated as @cM @trDMr[cMh( ) (cSecLe)cLe] r[1h( )]DMcM RgTFnMr'l =RM csolid(28) Similarly, the driving force for other species, i, is given by, i=@ E @ci=RgTlnci+0 i+F'lni (29) 10 and the ux Jiis de ned as, Ji=[1h( )]Dici RgTri=[1h( )]Dirci[1h( )]Dici RgTFnir'l(30) The resulting mass transport equation can be obtained by inserting Eq. (30) into the balance equation (25), rendering a phase eld-dependent form of the Nernst-Planck equation: @ci @trf [1h( )]Dircigr[1h( )]Dici RgTFnir'l =Ri(31) Finally, for both RMandRi, we introduce a generalized form of the reaction term, which is given by Rj=mtX m=18 < :kmpY vjm>0cvjm jkmrY vjm<0cvjm j9 = ;(32) wheremtis the total number of chemical equations, vjmis the stoichiometric coecient for species jin the chemical reaction m, andkmpandkmrrespec- tively denote the rate constant of products and reactants in reaction m. Note that we de ne vjm>0 for products and vjm<0 for reactants. Also, we empha- size that when the metal ion is involved, the actual concentration c0 M=cMcsolid must be adopted in Eq. (32). 2.3.3 Electrostatic potential The distribution of electrostatic potential 'lcan be estimated by the following Poisson-type equation [28, 29], r(r'l) =nMFcsolid@ @t(33) In Eq. (33), the variation of charge density due to the chemical reaction shown in Eq. (1) is accounted for by de ning an additional term on the right- hand side, with the term csolid@ =@t capturing the creation of electrons due to the dissolution of the metal electrode. Also, is the electric conductivity, which is de ned as =h( )s+ [1h( )]l (34) wheresandlare the conductivity in solid and liquid phases, respectively. The magnitude of the solid conductivity is chosen to be a suciently large value (s= 1107S/m), so as to ensure a uniform distribution of 'lin the solid phase. Thus, 'l, which is indistinctly referred to as electrolyte or electrostratic potential, is solved for in the entire domain but its magnitude is only relevant within the electrolyte and at the electrolyte-electrode interface. 11 Then, upon the assumption of a dilute solution [30], lis given by the following concentration-dependent function l=F2 RgT cMcsolidDMn2 M+X iciDin2 i! (35) 2.3.4 Mechanical deformation The governing equation of the mechanical system is derived by minimizing the strain energy density M= 0. Let us neglect for simplicity body loads and external tractions. Accordingly, see Eqs. (13)-(15), one reaches: Z h( )0:"dV=0 (36) with"="e+"pbeing the total strain tensor. By application of the Gauss divergence theorem and considering that Eq. (36) must hold for any arbitrary variations, we obtain the following balance: r[h( )0] =0 (37) where0is the undamaged or e ective stress tensor, which is given by 0=Cep: ("e+"p) =Ce:"e(38) withCepbeing the elastic{plastic consistent material Jacobian. It follows that the homogenized or damaged Cauchy stress tensor is given by =@" M= h( )0. The solid is assumed to exhibit isotropic strain hardening, which is char- acterized by means of the following power law relationship between the ow stressand the equivalent plastic strain "p: =y 1 +E"p yN (39) whereEis the Young's modulus, yis the initial yield stress and Nis the strain hardening exponent (0 6N61). 2.4 Dissolution-mechanics interactions Two important physical couplings are relevant to our theory. Firstly, the evolution of localized corrosion will result in material damage and redistribu- tion of mechanical elds, see Eq. (37). Secondly, the mechanical deformation of the solid will impact metallic dissolution by enhancing corrosion rates (mechanochemical theory [21]) and by fracturing the passivation layer (FRDR mechanism [22]). The latter is captured by enhancing the phase eld mobility 12 coecient, as proposed by Cui et al. [12]. Thus, the mechanical work required to fracture the passivation lm is characterised using the equivalent plastic strain"p, such that lm rupture will occur when "preaches a critical value "f. After a lm rupture event, passivation will result in the deposition of an initially unstressed oxide layer on the newly exposed metallic surface. Thus, corrosion rates are the result of a competition between the kinetics of passi- vation and straining. Accordingly, the relevant time scale over which strains accumulate is the time that it takes for the new passivation layer to rupture, since its deposition. This rupture-dissolution-repassivation cycle time interval is here denoted tiand accordingly, "p i="fwith"p i=Zti 0_"pdt (40) Once lm rupture occurs, the bare metal interface kinetics coecient Lais immediately recovered. Bare metal interface kinetics are sustained for a certain period t0, as it takes time for the passive lm to be suciently stable to impact corrosion kinetics. Once the lm is suciently stable, disso- lution rates are gradually reduced, with the decay process being dependent on the environment-material system, and captured here by means of a stabil- ity parameter k. After a time tfsince the decay starts, lm rupture occurs again,"p i="f, and a new FRDR cycle begins. Thus, the time interval for each rupture-dissolution-repassivation cycle equals ti=t0+tf, withtfbeing deter- mined by the mechanical elds. Accordingly, in terms of the corrosion current density, each rupture-dissolution-repassivation cycle is given by, if(ti) =( ia;if 0<ti6t0 iaexp [k(tit0)];ift0<ti6t0+tf(41) In addition, and independently of the FRDR process, corrosion kinetics are accelerated by mechanical elds [21]. Following Gutman's mechanochemical theory [31], we introduce an additional term kmto describe this phenomenon. Thus the corrosion current density reads, i(t) =km("p;h)if(ti) ="p "y+ 1 exphVm RgT if(ti) (42) wherehis the hydrostatic stress and Vmis the molar volume. The latter is de ned asVm=Mm=m, such that for a stainless steel with density m= 7930 kg=m3and molar mass Mm= 0:056 kg/mol, the molar volume equals Vm= 7:1106m3=mol. Accordingly, building upon the connection between the mobility coecient Land the corrosion current density i, a generalized Lcan be de ned that incorporates: (i) the FRDR mechanism, via(41); (ii) the sensivity of corrosion kinetics to mechanical elds, via(42); and (iii) the impact of the overpotential 13 on the corrosion current, via(19)-(21). Hence, L=L("p;h;) =( km("p;h)La();if 0<ti6t0 km("p;h)La() exp [k(tit0)];ift0<ti6t0+tf (43) 2.5 Species and reactions in the electrolyte The transport of ionic species and the homogeneous chemical reactions in the electrolyte have an impact on localized corrosion. Here, we assume that the electrolyte is a NaCl-based solution containing the following six ionic species: Mn+ M, M(OH)(nM1)+ , H+, OH, Na+and Cl. These result in the following chemical reactions: Mn+ M+ H 2Okf1*)kb1M(OH)(nM1)+ + H+(44) H2Okf2*)kb2H++ OH(45) Thus, the reaction term for each ionic species in Eq. (32) can be re-written as: RM=kf1cMcsolid+kb1cHcMOH RM(OH) =kf1cMcsolidkb1cHcMOH RH=kf1cMcsolidkb1cHcMOH +kf2kb2cHcOH ROH=kf2kb2cHcOH(46) wherekfandkbrespectively denote the rate constants for the forward and backward reactions, with the subscripts 1 and 2 being employed to distinguish between the reactions (44) and (45), respectively. Chemical reactions typically occur over much shorter times scales than mass transport and, as a result, the reactions are typically assumed to be in equilibrium. Under equilibrium condi- tions, the concentrations involved must remain proportional to each other and consequently an equilibrium constant can be de ned for each of the reactions being considered; i.e., here one nds1 K1=kf1 kb1=cHcMOH cMcsolid; K 2=kf2 kb2=cHcOH (47) Two approaches are typically followed to introduce the equilibrium assump- tion in the Nernst-Planck equations (31) [32]. One can solve for some of the ionic species assuming Ri= 0 and then estimate the remaining concentrations viatheir equilibrium relationships. For example, focusing on reaction (45), a 1Note that, for dimensional consistency, some authors choose to de ne a unit activity concen- trationca0, such that K2=cHcOH=c2 a0, withca0usually taken to be 1 mol/L. Here, we choose to drop this term for simplicity. 14 numerical solution for cHcan be obtained by solving Eq. (31) with RH= 0, and then Eq. (47b) can be used to estimate cOH, as the magnitude of K2is known. Alternatively, the magnitude of Ri, the source term in Eq. (31), can be adequately chosen so as to ensure that the equilibrium conditions are ful lled. I.e., for the case of reaction (45), the magnitudes of RHandROHare chosen such that the numerical estimates for cHandcOHalways satisfy cHcOH=K2. Here, the latter approach is adopted. Thus, following Eq. (47), the expressions for the reaction terms (46) can be re-formulated as RM=kb1(K1cMcsolid+cHcMOH) RMOH =kb1(K1cMcsolidcHcMOH) RH=kb1(K1cMcsolidcHcMOH) +kb2(K2cHcOH) ROH=kb2(K2cHcOH)(48) and by inserting (48) into (28) and (31), the governing equations for all concerned species are obtained. Here, it is important to note that the local equilibrium assumption for the chemical reactions implies that one does not need to know the magnitude of the backward and forward reaction rates ( kfi andkbi), it suces to know the value of the equilibrium constants Ki. Con- sider for example Eq. (48)d, the term kb2acts as a penalty term as increasing its magnitude will constraint the solution to ensure that cHcOH=K2is met. Accordingly, the choice of kbiis purely numerical, with the equilibrium condi- tion been enforced for kbi!1 . As discussed in Section 4.1, the magnitude ofkbiis chosen to be suciently large to approximate equilibrium conditions but not so large so as to induce convergence problems. 2.6 Summary of governing equations The balance equations can be particularized upon the consideration of the constitutive choices made in Sections 2.3, 2.4 and 2.5. A summary of the governing equations is provided in Table 1. This overview emphasises the cou- plings between the di erent elements of our theory. First, it can be observed that mechanics plays a role in the evolution of the corrosion front, viathe termL("p;h;) in the phase eld evolution equation (T.1). Secondly, the evolution of the corrosion front leads in turn to a degradation of the material sti ness and a re-distribution of the mechanical elds, see (T.5). Thirdly, the phase eld evolution equation is impacted by the electrostatic potential viathe dependency of the mobility coecient on the overpotential ( L("p;h;)), with the overpotential being related to the electrostatic potential 'lthrough Eq. (22). Due to electromigration, the electrostatic potential also has an impact on the transport of solid phase ions, see (T.2), and on the transport of the electrolyte ionic species, see (T.3). These transport equations for ionic species also contain a phase eld dependent term, to ensure that transport is lim- ited to the electrolyte. Finally, the calculation of the electrostatic potential is in uenced by both the phase eld, as (T.4) accounts for the creation of 15 Table 1 Summary of the governing equations of the electro-chemo-mechanical phase eld corrosion model presented. Phase eld: @ @t=L("p;h;) r2 @ E @  (T.1) Transport of solid phase ions: @cM @trDMr[cMh( ) (cSecLe)cLe]rn [1h( )]DMcM RgTFnMr'lo =RM csolid(T.2) Transport of ionic species: @ci @trf [1h( )]Dircigrn [1h( )]Dici RgTFnir'lo =Ri (T.3) Electrostatic potential distribution: r(r'l) =nMFcsolid@ @t(T.4) Mechanical balance: r[h( )0] =0 (T.5) with: @ E @ =2A[cMh( ) (cSecLe)cLe] (cSecLe)h0( ) +wg0( ) (T.6) L("p;h;) =8 >>< >>:km("p;h)L0expaanMF RgT ;if 0<ti6t0 km("p;h)L0expaanMF RgT exp [k(tit0)];ift0<ti6t0+tf(T.7) =h( )s+ [1h( )]F2 RgT cMcsolidDMn2 M+P iciDin2 i (T.8) 0=Cep: ("e+"p) (T.9) electrons, and the concentration of ionic species (T.2)-(T.3), due to the in u- ence of those on the electrolyte conductivity, as shown in (T.8). Thus, the electro-chemo-mechanical system is fully coupled through these interactions. 3 COMSOL implementation The electro-chemo-mechanical phase eld formulation presented in Section 2 is implemented in the nite element package COMSOL MULTIPHYSICS. The primal elds and nodal degrees-of-freedom (DOFs) are the phase eld order parameter , the displacement components u, the concentration of metal ions 16 cM, the concentrations of the ionic species involved ci, and the electrostatic potential'l. In COMSOL, as a result of its symbolic di erentiation capa- bilities, the governing equations of the model can be formulated in either the weak or the strong form, without the need to provide explicit expres- sions for the residuals and tangent sti ness matrices. As described below, ve COMSOL physics interfaces are used in the implementation, two in-built ones ( Solid Mechanics and Transport of Diluted Species ), and three user-de ned interfaces that exploit the Mathematics module to simu- late the evolution of the phase eld, the transport of metal ions, and the distribution of electrostatic potential. So-called State Variables are used to capture the lm rupture-dissolution-repassivation mechanism. The COMSOL implementation is made freely available at www.empaneda.com/codes. 3.1 Module setup Phase eld (T.1) . We use the Coefficient Form PDE interface to de ne the evolution of the phase eld order parameter . When using this interface, COMSOL provides the following generic form to de ne PDEs that contain derivatives up to second order in both time and space, ea@2 @t2+da@ @t+r(cr  + ) + r +a =f (49) To mimic (T.1), the PDE coecients are chosen as da= 1=L,c= , f=@ E @ ,ea=a= 0, and = = =0. Transport of metal ions (T.2) . We use the General Form PDE interface to de ne the transport of solid phase ions. This interface enables de ning di erential equations of the following form: ea@2cM @t2+da@cM @t+r=f (50) Accordingly, Eq. (50) can be particularised to (T.2) upon the following choices:da= 1,ea= 0,f=RM/csolidand xi=DM@cM @x+DM(cSecLe)@h( ) @x[1h( )]DMcM RgTFnM@'l @x(51) The source term of Eq. (50) contains the reaction variable RM, whose evolution is described by the algebraic equation (48), which is de ned as a COMSOL Variable . Transport of ionic species (T.3) . The Nernst-Planck equations are implemented using the Transport of Diluted Species interface, which is part of the Chemical Species Transport module. Di usion and migration 17 in electric eld are considered to be the main transport mechanisms (i.e., neglecting convection) and we particularise the implementation to ve species (Dependent Variables ) with concentrations cH,cCl,cOH,cMOH, andcNa. As per Eqs. (48)b-d, three associated reactions are de ned: RH,ROH, and RM(OH) . The coupling term [1 h( )] is incorporated viathe de nition of the di usion coecients and the electric potential is de ned to be the electrolyte potential'l. Electrostratic potential distribution (T.4) . We use the Poisson's Equation interface to de ne the distribution of electrostratic potential 'l. The Poisson's equation is written as: r(cr'l) =f. Thus, we de ne c= andf=nMFcsolid@ =@t to mimic (T.4). Mechanical balance (T.5) . The Solid Mechanics interface is used to implement the governing equation of the mechanical problem (T.5). The degradation of the material sti ness associated with the evolution of the corrosion front ( ) is captured at the weak form level by editing the Weak expression . If the body is assumed to be linear elastic, then one can attain the same e ect by degrading only the Young's modulus. However, in the case of an elastic-plastic material, such a strategy would imply that the degree of plastic ow is estimated using the homogenized (degraded) stress, whereas it is more common to use the e ective (undegraded) stress [33, 34], as it corresponds to the actual stress acting on the undegraded area. Finally, two State Variables are used to implement the FRDR mech- anism, (T.7). To this end, two history-dependent variables are de ned that respectively store the last values of equivalent plastic strain and time at which the last lm rupture event took place: "c=( "p;if ("p"c)>"f "c;otherwise(52) and tc=( t; if ("p"c)>"f tc;otherwise(53) Hence,"p"cis a measure of the accumulated equivalent plastic strain since the last lm rupture occurrence, and the time tiin (T.7) equals to ttcin this de nition. Note that tcis updated at the end of each FRDR cycle, such thattibecomes zero at the beginning time each FRDR cycle. Also, at the beginning of the calculation, we de ne "c(t= 0) =tc(t= 0) = 0. Thus, storing "cenables determining when the lm rupture event occurs ( "p"c="f) and, usingtc, this information is used to identify what stage of the FRDR cycle corresponds to the current time t. Accordingly, the interface mobility is de ned 18 in COMSOL as L("p;h;) =km("p;h)L0expaanF RgT minf1;exp [k(ttct0)]g (54) such thatLcorresponds to the mobility coecient of the bare metal if the time since the last rupture event ( ttc) is smaller than the time that it takes for the lm to stabilise ( t0), and is degraded otherwise. 3.2 Solution strategies The time-dependent solution step is used. COMSOL provides both monolithic and staggered solution algorithms, respectively termed Fully Coupled and Segregated . While monolithic approaches are appealing due to their uncon- ditional stability, they generally lead to a poor performance convergence-wise, unless quasi-Newton methods such as BFGS are used [35{37], which are cur- rently not available in COMSOL. Thus, a staggered approach is adopted here. First, the charge balance problem (T.4) is solved to obtain the distribution of electrostatic potential 'l. Then, the Nernst-Planck equations (T.3) are evalu- ated and the concentrations of M(OH)(nM1)+ , H+, OH, Na+, and Clare obtained. Then, the deformation of the solid is estimated by solving the balance of linear momentum, (T.5). Finally, solutions are obtained for the phase eld evolution equation (T.1) and the transport of solid phase ions (T.2). Converge is assessed for all elds in what is usually referred to as a Jacobi-type multi- pass solution approach. This is de ned by specifying a number of iterations in theSegregated node settings. 4 Results 4.1 Pencil electrode test with electrochemistry We shall rst employ the proposed electro-chemo-mechanical phase eld formu- lation to simulate the so-called pencil electrode test, a paradigmatic benchmark in corrosion science. As shown in Figure 2, a stainless steel wire with a height of Hs= 150m and a diameter of d= 25m, is mounted into an epoxy coating, leaving only the top edge of the sample exposed to the corrosive solution. This boundary value problem has been investigated using multiple computational techniques, including phase eld [8, 12, 27, 38], but without incorporating the electrochemical behaviour of the electrolyte. Here, we explore the impact of the electrochemical process on the pencil electrode test by considering the electrochemical parameters listed in Table 2 (unless otherwise stated). In previous phase eld corrosion analyses of this benchmark test (see, e.g., Refs. [12, 38, 39]), only the stainless steel wire was modeled, with the Dirichlet boundary conditions =cM= 0 being applied at the top edge. This scenario corresponds to assuming that the environment near the corrosion front is com- pletely isolated [40]. It is unclear whether this simpli cation is sensible, even 19 d epoxymetalElectrolyte Hl Hs Fig. 2 Schematic description of the pencil electrode test, including the boundary conditions adopted at the top edge. Table 2 Electrochemical parameters for the pencil electrode test, values taken from Refs. [4, 12, 28]. Parameter Value Unit Interface energy 10 J =m2 Interface thickness ` 0.005 mm Temperature T 293:15 K Di usion coecient DM;Di 8:5104mm2=s Average charge number nM 2:1 | Molar volume Vm 7:12106m3=mol Free energy density curvature A 53.5 N =mm2 Average concentration of metal csolid 143 mol =L Average saturation concentration csat 5.1 mol =L Anodic charge transfer coecient aa 0.5 | Chemical equilibrium constant K1 1:625104mol=m3 Chemical equilibrium constant K2 1:0108mol2=m6 Initial electrolyte pH 7 :0 | Initial NaCl concentration 1 :0 mol =m3 Applied potential 's 0:2 V for the case of occluded environments such as cracks. A less restrictive and more realistic approach is to simulate the electrolyte and enforce the boundary conditions =cM= 0 at the edge of the electrolyte. We shall explore the out- puts from both approaches here and in the context of a corrosion pit (Section 4.2). In this case study, we take advantage of our electro-chemo-mechanical model and simulate the electrochemical behaviour of the electrolyte, applying the Dirichlet boundary conditions =cM= 0 at the top edge of the elec- trolyte domain and varying its height Hlto assess its in uence. Regarding the initial conditions, we set (t= 0) =cM(t= 0) = 1 for the metal and (t= 0) =cM(t= 0) = 0 for the electrolyte, such that the dissolution will naturally initiate from the interface. In the absence of electrolyte (the case 20 Hl= 0) the Dirichlet boundary conditions =cM= 0 are prescribed at the top edge, and the initial conditions are de ned as (t= 0) =cM(t= 0) = 1 for the entire domain. Following Ref. [12], the interface mobility coecient is chosen to be L0= 2106mm2=(Ns), and one should note that neither lm rupture nor mechanical e ects are relevant in this case study. The element length in the expected corroded region is made to be at least 3 times smaller the interface thickness `, to ensure mesh objective results. 0 40 80 120 160 200020406080100120 Fig. 3 Pencil electrode test with electrochemistry: predictions of pit depth versus time tas a function of the electrolyte height H l(for a xed metal height H s= 150m). First, we investigate the interplay between the transport of solid phase ions and corrosion by neglecting the role of the electrolyte potential and the transport of the remaining ionic species ( 'l=ci= 08x). The results obtained are shown in Figure 3, in terms of the pith depth (in m) versus time (in s). The pit depth is de ned as the distance between the initial metal/electrolyte interface to the current corrosion interface. It can be observed that corrosion rates decrease when the initial electrolyte height H lincreases. In the limit case, Hl= 0, corrosion is clearly in the di usion-controlled regime, and thus the rate- limiting step is the transport of metal ions away from the pit boundary. These results can be rationalized based on the impact of H lon the gradient of the normalized concentration of metal ions cM. Decreasing Hlleads to a reduction in the magnitude of rcM, as well as in the corrosion rates. This is shown in Figure 4, where the distribution of the normalized concentration of metal ions cMis given along the vertical axis for two H lchoices. The concentration of 21 metal ions equals the concentration of atoms in the metal csolidfor lowy=Hs values, until the interface is reached ( 0:5), where the concentration of metal ions drops until reaching the interface saturation value csat. Then,cM progressively diminishes until becoming zero at the far end of the electrolyte. Hence, the larger the electrolyte, the smaller the gradient of cM. This is readily seen in the gure, with the slope of the last stage being more pronounced for the case of H l= 0, relative to the H l= 200m result. This leads to a reduction of the corrosion rates through Eqs. (T.2), (16) and (17). This can also be readily seen in the gure, as the corrosion front is closer to the bottom edge of the metal ( y= 0) for the case of H l= 0. 0.0 0.5 1.0 1.5 2.0 2.5 0.0 0.2 0.4 0.6 0.8 1.0 1 2 0.00 0.02 0.04 Fig. 4 Pencil electrode test with electrochemistry: normalized metal concentration cMdis- tribution along the vertical direction ( y=Hs). Results are obtained for a time of t= 100 s for two selected values of the electrolyte height (H l= 0 and H l= 200m). We shall now turn our attention to the role of the applied potential and the concentration of ionic species in governing corrosion kinetics. As shown in Table 2, we consider a solution with pH=7 and molarity 0.001 M NaCl; i.e., a bulk NaCl concentration of 0.001 mol/L (1 mol/m3). Thus, the initial values ofciatt= 0 can be determined as follows. First, the concentrations cMand cMOH are chosen to be zero, given that no dissolution occurs at t= 0. Secondly, the concentration cHis determined from the initial pH using the standard rela- tion pH=log10cH. FromcH, one can readily estimate the initial value of cOH using the equilibrium condition (47)b. Finally, cClandcNaare same as the 22 initial NaCl concentration (i.e., 1 mol =m3in this case). These choices automat- ically satisfy the electroneutrality condition: nMcMcsolid+P inici= 0. These initial values of ciare also prescribed at the top edge of the electrolyte domain through Dirichlet boundary conditions. In what follows, the results are com- puted using a reference interface kinetics coecient of L0= 0:001 mm2=(Ns) and an electrolyte height of H l= 200m. As discussed in Section 2.5, the equilibrium conditions (44)-(45) can be automatically satis ed when kb1and kb2are large enough. However, one shall note that high values of kb1andkb2 can also worsen numerical convergence. After a sensitivity study, we nd that kb1= 5102m3=(mols) andkb2= 5103m3=(mols) are suciently large to e ectively satisfy the equilibrium conditions without hindering convergence. 0 40 80 120 160 200020406080100 Fig. 5 Pencil electrode test with electrochemistry: predictions of pit depth versus time tas a function of applied potential 's. The results obtained for a varying applied potential 'sare shown in Figure 5, where'sspans the range 0.1 to 0.3 V. In agreement with expectations, corrosion occurs faster with increasing 's, as a result of the enhancement of the overpotential , see Eqs. (21) and (22). For all cases, corrosion is activation- controlled with a linear relation between pit depth and time, revealing that the electrostatic potential 'ldoes not change signi cantly in time. Figure 6 further compares the distributions of 'lwith two typical applied potentials, 0.1 V and 0.3 V. Results show that the electrostatic potential 'lremains constant in the solid phase and then drops as we move away from the interface 23 0.020 0.04 0.06 0.08 0.10 0.12 0.14 0.16 0.18 0.0040 0.008 0.012 0.016 0.020 0.024 0.028 0.032 0.036Interface Interface Fig. 6 Pencil electrode test with electrochemistry: contours of electrostatic potential 'l aftert= 200 s for two selected values of the applied potential 's. in the electrolyte region. We emphasize that all electrolyte-related variables are solved for in the entire domain, with the phase eld being used to track the location of the electrode-electrolyte interface. Thus, the distribution of 'lin the solid phase ( = 1) has no practical interest. Another interesting observation from Figure 6 is the increase in magnitude of 'lwith increasing 's, as readily observed when comparing legends. This can be understood by considering the dual role that the applied potential plays. On the one hand, corrosion kinetics are accelerated for larger values of the applied potential, as shown in Figure 5 and discussed above. On the other hand, see (T.4), the increased dissolution rate ( d =dt ) resulting from higher applied potentials translates into a higher gradient of 'land, given that l= 0 at the top edge, in a higher value of the electrolyte potential. We conclude this benchmark test by assessing the in uence of initial NaCl concentration. As evident from (T.4) and (T.8), changing the initial concen- tration of solvent species would a ect the distribution of electrostatic potential 'lby changing the electrolyte conductivity, and this has an impact on the overpotential, viaEq. (22), and on the interface mobility coecient (T.7). To examine the role of solute concentration, we x the applied potential to be 's= 0:2 V and consider four selected values of the NaCl content: 0.001 M, 0.05 M, 0.1 M, and 0.5 M. The results obtained are shown in Figure 7, in terms of pit depth versus time. In agreement with experimental observations [41], sim- ulations predict an increase in corrosion rate with bulk NaCl concentration. The results shown in Figs. 5 and 7 demonstrate that the present electro-chemo- mechanical phase eld formulation can capture the sensitivity of corrosion rates to the applied potential and the solute concentration, respectively. 24 0 40 80 120 160 200 0 15 30 45 60 75 Fig. 7 Pencil electrode test with electrochemistry: predictions of pit depth versus time tas a function of the initial NaCl concentration. 4.2 Mechanically-assisted corrosion from a semi-circular pit: electro-chemo-mechanical analysis We now turn our attention to the interplay between electrochemistry and mechanics. A rectangular plate of width W= 0:3 mm and height Hs= 0:15 mm is considered. Plane strain conditions are assumed. The plate contains a circular pit of radius 10 m and is exposed to a NaCl-based electrolyte solution that extends over a domain with height Hl= 0:3 mm. The geometry of the metal-electrolyte system is shown in Figure 8, together with the initial and boundary conditions. As shown in Figure 8, we assume that the pit has nucleated as a result of the localised rupture of the passive lm, and therefore consider the presence of a passive lm along all regions of the metal-electrolyte interface outside the pit region. Within the void region ( R= 10m), we assume the passive lm is less stable, allowing the FRDR process described by (T.7) to be applicable. Conversely, the passive lm outside of the void region is assumed to be stable enough to preclude corrosion and the transport of ionic species. This is implemented by separating the electrolyte and the metal over those protected regions with an impenetrable and non-corroding 1 m layer, as indicated by the bold line in Figure 8. For simplicity, the same electrochemical parameters used in the pencil elec- trode test case study are adopted (Table 2), unless otherwise stated. The 25 RPassive ﬁlm Stainless steelElectrolyteBCs: xy Fig. 8 Mechanically-assisted corrosion from a semi-circular pit: geometric setup, initial and boundary conditions. applied potential is chosen to be 's= 0:4 V. The Dirichlet boundary con- ditions =cM= 0 and'l= 0 are prescribed at the top edge. The initial values ofciand the associated Dirichlet boundary conditions at the top edge of the electrolyte are same as those in the pencil electrode test. As in Ref. [24], the constitutive behaviour of the stainless steel is characterized by a Young's modulus of E= 190 GPa, a Poisson's ratio of = 0:3, a yield stress of y= 400 MPa, and a strain hardening exponent of N= 0:1. Film rupture, dissolution and repassivation is assumed to take place along the pit surface, and is characterized by the following FRDR parameters k= 0:0001,t0= 50 s, and"f= 0:003 [24]. A constant mechanical load is applied by constrain- ing the left side of the metallic sample and prescribing a constant horizontal displacement of u1= 0:8m in its right edge. The entire metal-electrolyte system is discretized by means of approximately 20,000 triangular linear nite elements, with the mesh being particularly re ned in the expected corrosion region, where the characteristic element length is at least three times smaller than the interface thickness `. Let us now estimate the reference mobility coecient L0by exploiting Eq. (18) to establish a connection with the experimentally-measured value of exchange current density i0. To this end, we conduct a simulation considering 26 only phase eld corrosion (T.1) and the transport of metal ions (T.2), such that Lis a constant independent of the mechanical elds and the overpotential. By choosing a suciently small value of Lwe make sure to be under activation- controlled corrosion conditions and proceed to measure the velocity at which the corrosion front evolves, vn. For the choice of L= 1103mm2=(Ns), the nite element model predicts a corrosion velocity of vn= 0:0032m/s. The corresponding corrosion current density ican then be estimated using Faraday's second law: vn=vn=i nMFcsolid(55) and this leads to a proportionality constant of = 1:081011in (18). Using the experimentally measured exchange current density i0= 1:422 mA=cm2 reported in Ref. [42], the reference interface mobility parameter is found to be L0= 1:51010mm2=(Ns). (a) (b) (c)1.0 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0.0 Fig. 9 Mechanically-assisted corrosion from a semi-circular pit: phase eld contours illus- trating the pit evolution over three selected time stages: (a) t= 1800 s, (b) t= 3600 s and (c)t= 5400 s. The results obtained are shown in Figure 9, in terms of the phase eld contours, which illustrate the evolution of the corrosion front. A more quan- titative description of the corrosion pit is shown in Figure 10, where the pit growth along the depth and width directions is shown. Here, pit depth refers to the growth along the vertical direction (or = 90for a polar coordinate sys- tem centred at the pit mouth), while pit width refers to the growth along the horizontal direction ( = 0). The results show that pitting is not symmetric, with corrosion rates being faster along the depth direction due to the role of mechanical elds in enhancing corrosion kinetics - see Eqs. (42) and (43). Dif- ferences are notable, with the pit extending along the vertical direction more than twice the pit width after 5400 s. Moreover, results reveal that corrosion rates reduce with time, deviating from the linear behaviour reported in the results obtained for the pencil elec- trode test (Figures 5 and 7). This is despite L0being much smaller in the present case study; L= 1:51010mm2=(Ns) vsL= 0:001 mm2=(Ns) for Figures 5 and 7. It is of interest to investigate the source of this drop in cor- rosion rates, which can be potentially caused by three factors: (i) the role of lm passivation, (ii) a shift from activation-controlled corrosion to di usion- controlled corrosion, and (iii) a raise in electrolyte potential 'l. Inspection of thecMdistribution shows that the concentration of metal ions at the inter- face is below the saturation concentration csat; hence, we proceed to disregard 27 0 900 1800 2700 3600 4500 540001530456075 Fig. 10 Mechanically-assisted corrosion from a semi-circular pit: predictions of pit exten- sion along the depth (vertical) and width (horizontal) directions. For a polar coordinate system (r;) centred at the pit mouth, the depth direction corresponds to r; = 90while the width direction corresponds to r; = 0. (a) (b) (c)0.0030 0.0032 0.0034 0.0036 0.0038 0.0040 0.0042 0.0044 0.0046 0.0048 0.0050 Fig. 11 Mechanically-assisted corrosion from a semi-circular pit: contours of electrostatic potential'lat times: (a) 1800 s, (b) 3600 s, and (c) 5400 s. (ii). Then, the role of the electrostatic potential is investigated by plotting the contours of 'lfor relevant time intervals. The results are shown in Figure 11 for timest= 1800 s,t= 3600 s, and t= 5400 s, with the solid domain ( = 1) marked with a grey colour. It can be seen that, while 'levolves with time, the change in magnitude is very small, of roughly 0.001 V from 1800 s to 5400 s. Moreover, the electrolyte potential decreases with time, which should lead to a higher dissolution rate, as per (T.7) and (22). It remains to assess the role of lm repassivation. To this end, we vary the magnitude of the stability parameter k, considering three selected values: k= 0:0001,k= 0:0002 andk= 0:0003. The results obtained are shown in 28 Figure 12 in terms of pit depth (i.e., growth along the vertical direction) versus time. In agreement with expectations, the more stable the lm (larger k), the smaller the corrosion rates. Results reveal signi cant sensitivity to changes in k, suggesting that lm repassivation is the largest contributor to the observed reduction in corrosion kinetics with time. 0 900 1800 2700 3600 4500 5400 0 10 20 30 40 50 60 Fig. 12 Mechanically-assisted corrosion from a semi-circular pit: pit depth versus time t predictions as a function of the lm stability parameter k. Next, we examine the distribution of the ionic species in the electrolyte. Figure 13 shows the contours of cMcsolid, pH, andcClfor a time of t= 5400 s. The contours show how metal ions accumulate close the corrosion front, how the pH remains relatively constant but changes signi cantly when approach- ing the electrolyte edge, and how the concentration of chloride ions decreases progressively from the pitting interface to the bulk solution. It is also worth noting that the concentration of metal ions cMcsolidremains well below the saturation concentration ( csat, see Table 2), indicative of activation-controlled corrosion conditions. Overall, it is clear that the concentrations of ionic species near the pit mouth are signi cantly di erent from those at the top edge of the electrolyte (the bulk solution). A more quantitative insight into the dis- tribution of the various ionic species is shown in Figure 14, where the natural logarithm of the concentration of relevant ionic species is shown as a function of their position along the vertical axis, going from the bottom of the pit to the top edge of the electrolyte. The results show how the concentration of metal 29 ions is highest near the corrosion front, due to metal dissolution, and then decreases as we move deep into the electrolyte. The accumulation of metal ions results in lower pH values (increasing cH) inside the pit due to the hydrolysis reaction (44). Consequently, as per reaction (45), the magnitude of cOHinside the pit becomes relatively small and increases as we approach the edge of the electrolyte. In addition, changes in cClare shown to be of secondary nature, but a higher Clconcentration is attained near the corrosion front to main- tain the solution charge balance. It is important to note that the accumulation of H+at the pit surface may favour hydrogen uptake into the metal and thus hydrogen embrittlement damage mechanisms [43, 44]. (a) (b) (c) 4.50 4.05 3.60 3.15 2.70 2.25 1.80 1.35 0.90 0.45 0.00 7.0 6.6 6.2 5.8 5.4 5.0 4.6 4.2 3.8 3.4 3.0 1.20 1.18 1.16 1.14 1.12 1.10 1.08 1.06 1.04 1.02 1.00 Fig. 13 Mechanically-assisted corrosion from a semi-circular pit: Contours of (a) cMcsolid, (b) pH, and (c) cCl. The results are shown for a time t= 5400 s. Finally, we investigate the di erences between the present modeling strat- egy, where the electrolyte domain is resolved, and the simpler, commonly used approach of neglecting the electrolyte and prescribing on the pit mouth the concentrations of ionic species in the bulk solution (see, e.g., [12, 20]). The dif- ferences between the predicted values of ciand'lat the pit surface and the top of the electrolyte (see Figures 11 and 13) suggest that results will be sensitive to this choice of modeling approach. Figure 15 shows the extension of pit depth along the vertical direction predicted both by accounting for the electrolyte domain (\electrolyte simulation" curve) and by replacing the electrolyte by Dirichlet boundary conditions for the concentrations and electrolyte potential at the pit mouth (\Pit mouth Dirichlet BCs" curve). The results show that accounting for the role of the electrolyte domain is important not only to rig- orously predict the pit chemistry but also to accurately estimate the resulting pitting. In particular, the simplistic approach of prescribing electrochemical boundary conditions at the pit mouth results in a sharper pit morphology that translates into an earlier pit-to-crack transition at the pit base. The slower corrosion rates predicted for the model resolving the electrolyte are likely asso- ciated with the drop in ionic concentration (as transport is enhanced), which results in a reduction in conductivity through (T.8) and, accordingly, to a slower corrosion rate. 30 100 200 300 400 -8 -6 -4 -2 0 2 4 100 200 300 400 0.00 0.03 0.06 0.09 Fig. 14 Mechanically-assisted corrosion from a semi-circular pit: predictions of ionic con- centrations along the y-axis, starting from the pit base and going all the way until the top edge of the electrolyte (as characterized by the distance to the bottom edge of the metal domain). The results are shown for a time t= 5400 s. 5 Conclusions We have presented an electro-chemo-mechanical phase eld framework for pre- dicting localized corrosion. Its main ingredients are: (i) a phase eld description of the evolution of the corrosion front, (ii) the modeling of mechanics phenom- ena such as the enhancement of corrosion kinetics with mechanical elds and the lm rupture-dissolution-repassivation (FRDR) process, and (iii) the char- acterization of electrochemical electrolyte behaviour, including the transport of multiple ionic species and the distribution of electrostatic potential. These three components are fully coupled and impact corrosion kinetics by means of a mechanics- and electrochemistry-dependent interface mobility parameter L("p;h;). The model is implemented in COMSOL MULTIPHYSICS, with the numerical implementation being described and made freely available to the reader at www.empaneda.com/codes. Numerical experiments were conducted on two case studies of particular interest to investigate the interplay between electrolyte behaviour and localized corrosion; key ndings include: The size of the electrolyte domain plays an important role on the corrosion process. If the electrolyte domain is chosen to be very small, then the transport of metal ions is hindered and, as a result, corrosion shifts from 31 0 900 1800 2700 3600 4500 5400 0 15 30 45 60 75 Fig. 15 Mechanically-assisted corrosion from a semi-circular pit: Comparison of pit growth along the vertical direction considering two modelling strategies: (i) resolving the electrolyte domain (\electrolyte simulation" curve), and (ii) neglecting the electrolyte and prescribing the concentrations of the bulk solution on the pit mouth (\Pit mouth Dirichlet BCs" curve). activation-controlled to di usion-controlled. Increasing the applied potential and the initial bulk solution concentration leads to a rise in dissolution rates. On the other hand, faster corrosion kinetics can potentially increase the electrolyte potential, reducing the overpotential and the magnitude of the mobility coecient. There are signi cant di erences between the bulk and local pit electrochem- ical behaviours. Among other, due to the hydrolysis reaction, an increase in the concentration of hydrogen ions is observed inside the pit, which can enhance hydrogen ingress into the metal and embrittlement. The simplistic yet widely used modeling approach of replacing the simula- tion of the electrolyte by prescribing the bulk concentration and electrolyte potential at the pit mouth overpredicts corrosion rates and underpredicts the time required for the pit-to-crack transition to occur. 32 Appendix A. Interface energy and thickness in phase eld corrosion Here, we follow a similar approach to that of Ref. [45] to derive the interface energy per area and the interface thickness `. In the absence of mechanical and chemical contributions, can be de ned as: =interface Ainterface=Z+1 1" wg( ) +1 2 d dy2# dy (56) where interfaceandAinterfaceare the total energy and the area of the interface, respectively. Note that we use yas the direction normal to the interface to be consistent with Sections 4.1 and 4.2. In equilibrium we have  = 0, such that Id dy@I @(d =dy ) = const (57) whereIis the integrand of (56). Combining (56)-(57), one reaches w 2(1 )21 2 d dy2 = const (58) Assuming that = 0 wheny!1 and that = 1 wheny!+1, from (58) we obtain d dy=r 2w (1 ) (59) De ning the location of the interface at = 0:5 (y=y0), the solution for reads =1 exph q 2w (yy0)i + 1(60) such that the interface thickness `is derived as, `=1 d =dy y=y0=r 8 w(61) Finally, by combining (56) and (59), the interface energy can be written as =Z+1 1" wg( ) +1 2 d dy2# dy= 2wZ1 0 2(1 )2r 2w1 (1 )d =r w 18 (62) 33 Disclosure statement No potential con ict of interest was reported by the authors. Acknowledgements C. Cui acknowledges helpful discussions with S. Kovacevic and M. Makuch (Imperial College London). C. Cui and R. Ma acknowledge nancial support from the National Natural Science Foundation of China (grants 52178153 and 51878493). E. Mart nez-Pa~ neda acknowledges nancial support from UKRI's Future Leaders Fellowship programme [grant MR/V024124/1]. C. Cui addi- tionally acknowledges nancial support from the China Scholarship Council (grant 202006260917). Notes on contributor(s) C. Cui: Conceptualization, Investigation, Methodology, Software, Writing { original draft, Writing { review & editing. R. Ma: Funding acquisition, Supervision. E. Mart nez-Pa~ neda: Conceptualization, Funding acquisition, Investigation, Methodology, Software, Supervision, Writing { review & editing.