Investigations of the intrinsic corrosion and hydrogen susceptibility of metals and alloys using density functional theory

1 Introduction

Physics-based simulation of materials aging and degradation through the mechanisms of corrosion and hydrogen embrittlement remains a challenging field of research (Barrera et al. 2018; Taylor 2015). Such research is motivated by the necessity for decision-makers in materials engineering and design to develop new materials that are more robust under a variety of environments, to establish fitness for service criteria, and to extend the lifetime of aging structures (Cordeiro et al. 2020; Taylor et al. 2018c). Credible models can help augment and interpret the often-limited set of relevant historical and experimental databases that typically inform such decision-making processes (Melchers 2019). Such models can fall into two generic categories: empirical models, with sometimes (but not always) physically motivated functional forms that are parameterized to laboratory and/or field data (Feliu et al. 1993a, 1993b; Melchers 2019; Mendoza and Corvo 2000; Oesch and Faller 1997); and first-principles models that are physically motivated and parameterized only by fundamental physical constants (Huang et al. 2019; Lin et al. 1981; Sharland 1987; Taylor 2015; Taylor et al. 2018c).

The development of first-principles models for corrosion and hydrogen embrittlement will require a careful re-construction of the key physicochemical processes that contribute to those phenomena (Gunasegaram et al. 2014). These processes include mass-transport from the environment to the near-surface region; chemical and electrochemical reactions at the electrified interface between the material and solution; compositional and microstructural changes in the material itself; interactions of solutes and surfaces with static and dynamic strain fields in the material; and formation and evolution of surface films and scales (Anderko et al. 2008, 2014, 2015; Barrera et al. 2018; Rice and Wang 1989; Shibata 2007; Tehranchi et al. 2020). These processes are dynamic, involve many different species, maybe close to or far from equilibrium, and heavily dependent upon environmental variables such as temperature, pH, electrochemical potential, and the presence of impurities in both the solution and the material itself. Clearly, deriving a completely first-principles approach to corrosion and hydrogen embrittlement phenomena is a challenge for theorists.

In this review, we describe the development of an approach that began with a relatively simple analysis from first principles of the electrode/electrolyte interface (Filhol and Neurock 2006) and moved up to consider the effect of different environmental variables on chemical transformations that take place on activated metal surfaces (Taylor et al. 2007c). The approach was then adapted incrementally to increase in complexity, such that electrochemical reactions on alloy materials in more complex solution chemistry environments can now be estimated entirely from first principles (Ke and Taylor 2020b; Ke et al. 2020b). The method’s assumptions and idealizations are reviewed in detail and potential improvements are suggested that could be made in future developments of the approach. Some new frontiers for the development and application of this methodology are also presented.

2 Density functional theory

Density functional theory (hereafter, DFT) is a method for predicting the structure and properties of matter based on the solution for the ground-state electron density distribution for a given potential field (Jones 2015). In the typical kinds of DFT calculations being performed in contemporary science and engineering, these potential fields often consist of the superposition of individual potential energy fields that are mapped to the nuclei and core electrons of the atoms that constitute solids, liquids, and the molecules and ions in their environment (also known as effective core potentials or pseudopotentials) (Norskov et al. 2010). Additional terms can be added into the electronic structure problems that are solved by DFT to simulate other effects, such as electric fields across materials or interfaces (Hansen et al. 2008; Jacob 2007; Jinnouchi and Anderson 2008b; Karlberg et al. 2007; Lozovoi and Alavi 2007; Ohwaki and Yamashita 2001), the physics of spin–orbit coupling (Boettger and Ray 2002), magnetization and charge-localization constraints (Chen et al. 2014; Mills et al. 2003; Souvi et al. 2013; Scivetti and Persson 2013; Wodtke et al. 2008), self-energy corrections (Mattsson and Jennison 2002), van der Waals interactions (Antony and Grimme 2006; Grimme 2006; Tkatchenko et al. 2010, 2012), and explicit or implicit methods to simulate the effect of solvation on system energies (Marenich et al. 2009; Mathew et al. 2014; Miertuš et al. 1981; Taylor 2009). DFT calculations are time-consuming, even with state-of-the-art computer hardware, and therefore also expensive. For this reason, efforts to reduce the computational burden of DFT calculations, or to approximate their outputs using machine learning and/or coarse-graining, have been popular areas of research for some time (Fellinger et al. 2010; Hegde and Bowen 2017; Martinez et al. 2013, 2020; Schleder et al. 2019).

Because DFT calculations provide information about hypothetical configurations of matter at the level of electronic granularity, a large number of properties can be predicted (Hafner et al. 2006; Ke and Taylor 2019; Parr and Yang 1984; Saal et al. 2013; Yang et al. 1984). Two of the most commonly reported outputs from DFT calculations include ‘optimized geometries’ for atoms and molecules in various states of interaction or compound formation, and the electronic internal energy, which is useful for comparing different states for their relative thermodynamic stabilities and to estimate heats of reaction from the difference in internal energies of product and reactant states (Alfe et al. 2000).

DFT calculations are solved by self-consistent solution of the Kohn–Sham equations (Kohn and Sham 1965), which mirror the concept of ‘molecular orbitals’ in computational quantum chemistry, although the ideas do not exactly correspond (Stowasser and Hoffmann 1999). In the Kohn–Sham equations below, the first term in parentheses is the kinetic energy functional and the second term is the effective potential energy field that is interacting with the occupied one-electron Kohn–Sham orbitals ϕi(r). ϵ_i is the Eigenergy for the corresponding Kohn–Sham orbital. The electron probability density function ρ(r) is given by summing over the squared magnitude of the occupied individual Kohn–Sham orbitals. The total system energy is the sum of occupied orbital energies with some additional terms for the Coulomb energy of electron–electron interactions and some corrections from electron exchange and correlation (see following). N is the number of electrons.

Whereas the physics of the interaction between charged particles that are fundamental to DFT is formally exact, determining an accurate solution of these equations is limited by (a) the present lack of knowledge of what is called the exchange-correlation functional (but good approximations exist), and (b) the burdensome computational requirements needed to solve this problem which scale as n³, where n is the number of electrons (Korth 2017; Mortensen et al. 2005; Su and Xu 2017). The two effects of exact exchange and electron correlation relate to the quantum mechanical properties of electrons and the difficulty in describing the multi-electron wavefunction in terms of separable electronic orbital functions. A hierarchy of approaches has been developed to incrementally approach the exact energy limit, and increasing computational complexity is required to ascend this hierarchy, known in quantum chemistry as “Jacob’s Ladder” (Car 2016). Whereas this is one source of error in the formalism, which itself is an active field of research (Mardirossian and Head-Gordon 2017), when the method is being applied to practical engineering problems, such as corrosion and hydrogen embrittlement, other sources of error are probably more important and result from the necessity of creating idealized, solvable representations that mimic problems of industrial significance. Even comparison to experiments conducted under well-controlled laboratory conditions often involves many steps of approximation and idealization. For this reason, it is the creation, evolution, and interpretation of DFT models in the context of corrosion science that forms the focus of this review.

A flowchart showing the workflow for DFT calculations is given in Figure 1. A list of useful outputs from first-principles DFT models that can be applied to help interpret corrosion and hydrogen embrittlement is shown in Table 1, adapted from Ke and Taylor (2019).

Figure 1:

Workflow of DFT modeling.

3 Using DFT to model the electrode-electrolyte interface

Aqueous corrosion is a redox process involving the dissolution of metal atoms as cations (at the anode), and a counter-electrode reaction (at the cathode) that involves the reduction of species in solution with consumption of the electrons that are produced during anodic dissolution (Tomashov 1961). Some of the fundamental elements that need to be modeled from first principles, therefore, are the structure and properties of the electrochemical interface (Jinnouchi and Anderson 2008b; Lozovoi et al. 2001), and the influence of electrochemical and materials parameters on the reactions that occur there (Anderson and Kang 1998; Jinnouchi and Anderson 2008a; Taylor et al. 2018b). To construct models for materials/environment interfaces, an appropriate function ν_eff(r) needs to be constructed that creates the potential energy field around which the electronic structure will be solved via the self-consistent field theory. The predominant model used for such scenarios is the “periodic slab” method in which a three-dimensional box is constructed which has a finite number of layers of atoms that minimally represent the material and its surface (i.e., one could use say five layers of metal atoms arranged in a specific crystallographic orientation of interest; often low-index surfaces are used because of their close-packed, stable structure) and either a void or solvated space adjacent and normal to the surface is populated with the environmental molecules of interest (Makov and Payne 1995; Taylor et al. 2006b). This method has been widely used to simulate the catalytic properties of metal surfaces used in industrial catalysis and then later adapted to study systems of interest to corrosion science (Norskov et al. 2010; Skúlason 2015).

Motivated to study the impact of changing electrochemical potentials on the functionality of electrocatalysts and corrosion reactions, Filhol and Neurock (2006), and Taylor et al. (2006b) expanded upon this framework by introducing a virtual potentiostat into the cell (Figure 2). Electrification of the interface to simulate anodic (or cathodic) potentials was effectuated in this method by changing the number of electrons, N in equation (1b), as the constraint to the Kohn–Sham DFT model (Lozovoi and Alavi 2007; Lozovoi et al. 2001). The corresponding variation in electrochemical potential was detected in the simulations by monitoring the difference in the electrostatic potential between the Fermi level of the metal slab system against the potential at the point of the virtual potentiostat (located at the mid-point of the cell, in the environment region); i.e., essentially the work function of the coupled electrode/electrolyte system. This method was used to simulate the electrochemical activation of water, and the impact of changes in potential on the water molecule arrangements near the interface (Desai and Neurock 2003; Filhol and Neurock 2006; Rossmeisl et al. 2006; Taylor et al. 2006a, 2006b, 2007b, 2007c). Of particular interest was to evaluate the effect of surrounding, explicitly modeled water molecules, on the predictions made in the DFT approach. Because the inclusion of water molecules in the slab model adds significant computational complexity, many studies of adsorption exclude “spectator” solvent molecules. The study by Rossmeisl et al. showed that these can be omitted from the calculation with little effect on the chemisorption energies and final results (Rossmeisl et al. 2006). For other calculations, such as metal dissolution, in which solvation of the product state (i.e., a dissolved cation) plays a key role in lowering the system energy, such idealizations would be expected to radically influence the final results.

Figure 2:

An example of a periodic slab model with virtual potentiostats shown to enable calibration to a vacuum reference state (top), and then electrified reference states (bottom) (Taylor et al. 2006b). Reproduced with permission from the American Physical Society (© 2006). The electrostatic potential profile across the periodic slab is plotted, with the vacuum potential φ(v), the Fermi level φ(m) marked for the top figure, and reference points for the ‘double electrode’ method, φ′₀(m), marked and the solution potential φ′₀(w). The vacuum potential is used to align the Fermi level of the system relative to vacuum as a reference point for simulations on charged, electrochemically polarized slab systems.

The electrochemical activation of water corresponds to either the cathodic reaction (upon application of a negative potential) or the anodic reaction (for more positive potentials) as follows:

The superscript * denotes that the associated molecule or fragment is adsorbed to the metal surface.

When cathodic potentials are applied, the DFT calculations produce stable geometries for water molecules with the dipole of the water molecule rotating so that the hydrogen atoms become oriented toward the metal electrode (Taylor et al. 2007b). As the potential becomes more cathodic, bonds form between the hydrogen atoms and the metal, breaking the O-H bond in water and producing adsorbed hydrogen and OH⁻: this is an important step in both the cathodic evolution of hydrogen (Conway and Bockris 1957), which can be a counter electrochemical reaction to anodic dissolution in aqueous corrosion and also in the activation of hydrogen before absorption and embrittlement of the material (DeLuccia 1988). Conversely, when anodic potentials are applied, the DFT model correctly reproduces the opposite behavior: water breaks up to produce OH adsorbed to the surface and liberating H⁺ which forms H₃O⁺ and H₅O₂⁺ moieties in the interslab space, a process that has been observed to occur on metallic surfaces upon polarization in an aqueous environment via in situ surface probes (Maurice and Marcus 2018; Maurice et al. 2000). Figure 3 shows how the DFT energy varies with the electrochemical potential applied by the computational potentiostat model (Figure 2). As the potential varies from cathodic to anodic, the stable phases of the surface, in order, the transition from H adsorbed, to water adsorbed, and hydroxide and oxide, because of the cathodic and then anodic electrochemical activation of water (equations 2a, 2b) and then the subsequent electrochemical activation of OH.

Figure 3:

Water activation on Cu(111) using the computational virtual potentiostat method (Taylor et al. 2007b). Reproduced with permission from Taylor and Francis (© 2007). The plot shows the free energy for the Cu(111) system with adsorbed water (filled circles), hydroxide (open circles), oxide (open squares), or hydrogen (open triangles) as obtained from direct DFT calculations of solvated metal slabs.

Of relevance to corrosion was modeling the impact of changing the electrochemical potential of the system not just on the reactivity of water but also on the changes that resulted in the metal slab itself. Using the virtual potentiostat model, a Cu(111) surface was polarized incrementally to progressively more positive surface charge densities, and geometry optimizations performed at each step (Taylor 2012a). After activation of water to OH occurred, a second activation step to adsorbed O resulted from the geometry optimization.

Following this step, a Cu surface atom in the (111) plane was lifted out of the surface and displaced by the adsorbed O, effectively forming a stoichiometric layer of Cu₂O on the surface. The results indicated that DFT calculations could provide a powerful method for understanding the electrochemical destabilization of a/the metal surface and could be used in the future to investigate a variety of scenarios relevant to corrosion science and engineering (Figure 4).

Figure 4:

Cu atom lift-out from a Cu(111) slab with oxygen adsorbed occurs in the DFT simulated system when the electrochemical potentiostat anodic potential is increased to potentials beyond water and hydroxide activation (Taylor 2012a). Reproduced with permission under Creative Commons license, Christopher D. Taylor (© 2012).

4 Modeling adsorption on metal surfaces

The ability of DFT calculations to provide information about reactions and adsorption of molecules and atoms from the environment on to metal surfaces has resulted in their application to the study and comparison of modes of adsorption for chemical corrosion inhibitors (Borck et al. 2007; Cho et al. 1995; Costa and Marcus 2015; Durnie et al. 2005; Jiang and Adams 2003; Jiang et al. 2004; Kokalj et al. 2011; Kovačević and Kokalj 2012, 2013; Kovacevic et al. 2017; Lebrini et al. 2008; Motta et al. 2012; Peljhan and Kokalj 2011a; Sun et al. 2012), to oxygen adsorption as a function of electrochemical potential to form pre-passive monolayer films (Francis and Taylor 2013; German and Efremenko 2004; Hyman and Medlin 2005; Karlberg et al. 2007; Ke and Taylor 2020b; Ke et al. 2020b; Samin and Taylor 2017b, 2018b, 2018c; Taylor 2013; Taylor et al. 2018b; Yeh et al. 2009), and to the influence of interactions between co-adsorbates or alloying elements on processes such as hydrogen uptake (Greeley and Nørskov 2007b; Limmer et al. 2017; Taylor et al. 2011).

Taylor, Kurapati, and Mondal developed a model for molecular adsorption of corrosion inhibitors on to defect sites on passivated metallic surfaces by using a cluster model that calculated the binding energies of the inhibitor molecule to an M(OH)_n molecular fragment, building on earlier work by Mondal and Taylor (Mondal and Taylor 2014; Taylor et al. 2018a). The binding energy is computed by the difference inbound (M + A, signifying Metal plus Adsorbate, i.e., the chemical inhibitor molecule) and un-bound (M, A) systems:

Following this equation, more negative binding energies imply stronger chemical bonding because the system M + A has lower potential energy than the isolated parts. Inhibitors can be compared based on their relative binding energy: a more negative binding energy can imply a more effective chemical inhibitor, although other factors such as packing on the surface, availability in mixed-media environments, interactions with corrosion products, etc. may also play a role (Hackerman 1952, 1962; Taylor et al. 2015b). The field of application of quantum chemical methods (including DFT) to the problem of effective corrosion inhibitor selection and design, is proliferous and could entail a review entirely of its own. One approach, developed in the late 60s, involves the regression of computable parameters for molecules based on their electronic structure, against experimentally determined activities-the so-called quantitative structure-activity relationships (Durnie et al. 1999, 2005; Gece 2008; Vosta and Eliasek 1971). This is in effect a semi-empirical way to indirectly link the properties of inhibitor molecules to performance, but does not involve a first-principles mechanistic treatment, and so will not be further described herein.

Oxygen adsorption on metal surfaces is the precursor to the formation of passive films (Francis and Holby 2015; Seyeux et al. 2005). Oxygen binding energies are typically strongly negative values on the order of a few eV, indicating strong binding to metallic surfaces: hence, the propensity for metals to oxidize in the natural environment. Oxygen adsorption is the final step in electrochemical water activation (Taylor et al. 2007b, 2007c), but can also occur directly because of dissociative adsorption of O₂ on metal surfaces, a process which is exothermic and therefore spontaneous (Calvo 2014; German and Efremenko 2004; Janik et al. 2009; López-Moreno and Romero 2015; Norskov et al. 2004; Panchenko et al. 2004; Wang et al. 2008). The study of oxygen adsorption on metal surfaces using DFT has seen considerable interest because of the push to develop more efficient oxygen reduction catalysts for fuel cell technologies in the 2000s (German and Efremenko 2004; Hansen et al. 2008; Hyman and Medlin 2005; Janik et al. 2009; Karlberg et al. 2007; Mills et al. 2003; Norskov et al. 2004; Panchenko et al. 2004; Wang et al. 2011, 2012; Yeh et al. 2009). The nature of oxygen binding to metals that are challenging to study experimentally, for example, radioactive elements such as Tc and Pu, has also been studied by Taylor (2013, 2014 to understand more about the corrosion properties of these materials. As oxygen accumulates on close-packed metallic surfaces the adsorption energy changes because of lateral interactions between the adsorbates (Cheng et al. 2014; Tian et al. 2009). When atomic oxygen adsorbs on closest-packed metal surfaces, it typically resides in hollow-sites, bonded to three or four metal atoms, depending on the configuration (Francis and Holby 2015). The adsorption is accompanied by charge transfer. The surface metal atoms that are ‘bonded to’ or first-nearest neighbors to the oxygen atom each develop a partial charge, indicating that the metal surface, even at sub-monolayer coverage of oxygen, already has an oxidized character (Francis and Taylor 2013). As the coverage of oxygen increases, reconstructions can occur in the top layer of metal atoms because of this polarization of the metal atoms, and the repulsions that accumulate between the oxygen atoms.

A study was made of the accumulation of oxygen on the Ni(111) surface to produce monolayer, bilayer, and trilayers of NiO(111) for the purpose of characterizing the changes in valence, coordination, and reactivity that accompany the initial stages of growth of the oxide film (Olatunji-Ojo and Taylor 2013). Table 2 summarizes some of the findings of the paper, showing that, when an appropriately higher level of theory (PBE + U) is used (Dudarev et al. 1998a, 1998b), the so-called “octopolar” reconstruction of the Ni(111) surface to form the first partial monolayer of NiO(111) results from the model, which is a finding that is in agreement with the experimental observations (Barbier et al. 1999). As subsequent layers of the oxide phase form, the energy of oxygen adsorption approaches the energy of formation for bulk NiO (Table 2). The same convergence in energy occurs with the lower level of theory (PBE), but in that case, there is a reduced driving force to form the oxide because PBE underestimates the energy of formation for the oxide phase. Further such studies of these processes will help to elucidate the mechanisms that lead to passive film formation and growth on these metal systems, including identification of the relevant charge transfer processes (Yu et al. 2012), how the reaction energies change with the time evolution of the film structure and thickness, and identification of the role of alloy elements (Castle and Asami 2004; Han et al. 2015; Ke and Taylor 2020b; Ke et al. 2020b; Liu et al. 2015b; Marcus 1994; Seyeux et al. 2013; Yu et al. 2018), solution chemistry (Anderko et al. 2008; Bouzoubaa et al. 2009a, 2009b, 2010; Ke et al. 2020b; Liu et al. 2015a; Samin and Taylor 2018b; Taylor 2012b; Yu et al. 2019), the adhesion at metal/oxide interfaces (Jiang et al. 2008; Wang et al. 2006), and defect states at the metal/oxide and oxide/environment interface that are most relevant for limited and controlled oxide growth as a function of environment and materials composition (Krishnamurthy et al. 2002; Lin et al. 1981; Todorova and Neugebauer 2014, 2015; Samin and Taylor 2018a, 2018d; Seyeux et al. 2013).

The ability to study the interactions of adsorbates with one another on metal surfaces directly using DFT methods was used to investigate a long-standing question regarding the ability of certain chemical species such as As and S to influence hydrogen uptake kinetics on structural materials (McCright 1977; Protopopoff and Marcus 1987, 1988; Taylor et al. 2011). A Ni(111) surface was studied by Taylor et al. and hydrogen adsorption and absorption energies calculated using DFT (Taylor et al. 2011). Sulfur, a known hydrogen recombination poison, was then adsorbed on the simulated surface and its influence on hydrogen adsorption and absorption energies of nearby sites was determined. Using the energies extracted from the DFT calculations, the hydrogen recombination versus uptake (i.e., absorption into the metal) kinetics were estimated. It was shown that up to certain quantities of sulfur adsorbed on the surface, hydrogen absorption kinetics could be enhanced while recombination was quenched (i.e., the poisoning effect). The kinetics were estimated by raising or lowering the activation energy according to the shift in chemical potential of the adsorbed hydrogen because of its proximity to the sulfur atom. Raising the surface chemical potential by an amount δE was anticipated to lower the activation energy by an equal amount. Beyond a certain sulfur coverage, however, the ability for hydrogen to either adsorb or absorb was entirely quenched as full coverage of sulfur effectively sealed off the surface to further hydrogen adsorption (Figure 5).

Figure 5:

DFT calculations of sulfur and hydrogen co-adsorbed show two zones of influence of S on nearby H adsorption energies. Zone A increases the surface chemical potential mildly, lowering the activation energy enough to promote an increase in H absorption into the metal by the uptake enhancement factor, which is dependent upon electrochemical overpotential, h. Zone B increases the surface chemical potential to a much higher extent, creating a field that repulses hydrogen too strongly to allow for either adsorption or absorption to occur. At higher S coverages, there is more zone B than zone A on the surface, thus the uptake effect is coverage dependent (Taylor et al. 2011). Reproduced with permission from IOP Publishing (© 2011).

5 Electrochemisorption and surface Pourbaix diagrams from first-principles

One of the most fundamental tools in the electrochemical corrosion scientist’s toolkit is the Pourbaix diagram (Deltombe et al. 1974). Pourbaix diagrams express in a visual way the most stable phases of an element in the two-dimensional electrochemical phase space of electrochemical potential (usually in units of V NHE) and pH. The Pourbaix diagram is calculated using Gibb’s free energy of the various possible phases and the Nernst equation to adjust the free energies for the conditions of potential and pH. For example, the free energy for the reaction:

is calculated using the Nernst equation as follows:

Where, in this case, m = 2 for the number of H⁺ exchanged in the reaction and n = 2 for the number of electrons exchanged. RT is the ideal gas constant multiplied by temperature. F is Faraday’s constant, and U is the potential; U₀ is the equilibrium potential for the reaction. Owing to the typical linear form of equation (5), Pourbaix diagrams can be solved by linear programming methods and have the appearance of linear divisions of the U-pH two-dimensional space. The Pourbaix diagram, therefore, is a kind of map, where lines mark the transition from one stable phase of the metal/environment system to another. The various regions of the map are sometimes referred to as “corrosion”, “immunity”, and “passivity” depending on the properties of each phase.

Using chemisorption energies one can compute something akin to a Pourbaix diagram for the most stable chemical species that form on a metallic surface: this approach has been extensively developed since the 1980s by Protopopoff and Marcus (Marcus and Protopopoff 1993, 1997; Protopopoff and Marcus 2003, 2005, 2012). For example, water may exist in the undissociated molecular form of H₂O on a metallic surface under some conditions of potential and pH, but as the potential becomes anodic or cathodic it may transform according to the reaction sequences shown above (2a–c)a–c(2a–c). The chemisorption energies can be shifted for different conditions of pH and potential using the same rules as the Nernst equation:

Adsorption reactions that include electrochemical transfer can be referred to as electrochemisorption v(Conway et al. 1974).

The electrochemisorption windows of potential for water on a variety of metals were calculated by Taylor, Neurock, and Kelly, as shown in Figure 6 (Taylor et al. 2007c). More noble elements typically have a broad range of electrochemical stability for water, and, in particular Au. Pd and Pt, good hydrogen reduction catalysts, have underpotential regions of stability for hydrogen adsorption from water on the surface. Ni and Mo are electrochemically active toward water under all potentials studied. Mo, Ni, and Cu are also ‘earlier’ oxide formers in the sense that they allow the surface in contact with water to oxidize at lower electrochemical potentials. These first-principles predictions are in broad agreement with the known activities and behaviors of these metals (Henderson 2002; Taylor and Neurock 2005).

Figure 6:

Electrochemisorption windows for H (gray), H₂O (white), and OH (black) on a series of close-packed metallic surfaces. Standard conditions were assumed (i.e., pH 0, 298 K). A full Pourbaix style treatment would also include an additional axis of pH variation to this kind of phase diagram for each of the surfaces studied, and additional phases (Taylor et al. 2007c). Reproduced with permission from IOP Publishing (© 2007).

Protopopoff and Marcus have created surface phase diagrams for several materials/environment scenarios, including sulfur, hydroxide, thiosulfate, and lead adsorption on various metals and alloys. Their technique primarily used adsorption energies taken from surface science experiments (Marcus and Protopopoff 1993, 1997; Protopopoff and Marcus 1987, 1988, 2003, 2005, 2012). Surface Pourbaix diagrams were created based on DFT calculated adsorption energies for nickel in high purity water under both ambient and relevant temperatures, relevant to the nuclear industry (Taylor et al. 2006a, 2007a). The method has continued to be adopted by others, including Williams et al. (2015) and Yuwono et al. (2019) for Mg, and Kim et al. (2014) for Tc. An example of the Pourbaix diagram for Mg(0001) computed by Yuwono et al. is shown in Figure 7.

Figure 7:

A Pourbaix diagram constructed from DFT for the Mg(0001) surface by Yuwono et al. (2019). Reproduced with permission from Elsevier (© 2019). The various θ values correspond to surface coverages of the adsorbates (Mg*O, or Mg*OH, or Mg*H) used for the various colored curves in the plot.

6 Competitive electrochemisorption on metals

Going beyond the adsorption of water and its dissociation products, competitive chemisorption of other species can occur when ions such as chloride compete with water, oxygen, etc. for adsorption sites on the metallic surface (Conway et al. 1974). Chloride ions can adsorb electrochemically to the surface to form a chemisorbed chloride:

This process can compete with the adsorption of oxygen, hydrogen, hydroxide, or molecular H₂O. Analysis of the electron density maps produced using DFT shows that the Cl maintains its negative charge on adsorption through polarization of the metal atoms it contacts on adsorption (so the complex Cl–M(surf) could also be written as Cl⁻–M⁺ (surf) or even Cl⁻ (M₃)⁺, although Cl adsorption may not always be at threefold hollow sites depending on the surface topology and/or presence of defects) (Andryushechkin et al. 2015; Ke and Taylor 2020b). The change in free energy for each of these reactions can be computed with DFT, and then by using suitable thermodynamic models based on Born-Haber cycles and the knowledge of the reference equilibrium potentials tabulated for hydrogen evolution, chloride reduction potential, etc.(Taylor 2012b) This approach translates the zero Kelvin chemisorption energy (ΔE_ads, using ads more generally for adsorption-chemical (chemisorption) or physical (physisorption)) into ΔG at realistic temperatures for the reactions relevant to problems in alloy corrosion, crack growth, etc.

The Born–Haber cycle that could be used to compute the free energy for chloride adsorption onto a metallic surface starting with Cl⁻ in the aqueous phase (corresponding to Eq. (9)), ∆G^target, is shown in Figure 8. Similar cycles can be derived for hydrogen adsorption, hydrogen sulfide adsorption, water adsorption, etc.

Figure 8:

Example Born–Haber cycle for chloride adsorption onto a metallic surface from Cl⁻ in the aqueous phase. ∆E^abinitio is the adsorption energy from the quantum chemical calculation, and the various steps in the cycle refer to the thermodynamic parameters required to get the overall free energy of interest ∆G^target (Taylor 2012b). Reproduced with permission from NACE International (© 2012).

Using this data, Taylor computed phase diagrams for the competitive chemisorption of O, OH, H, H₂O, N, NH, NH₂, NH₃, and Cl on Fe(110) surface under conditions relevant to ammonium chloride corrosion, a problem of concern in the refining industry (Akpanyung et al. 2019; Alvisi and de Freitas Cunha Lins 2008; Taylor 2012b). In this failure mode, ammonia and hydrochloric acid condense onto the inner walls of refinery equipment to form the NH₄Cl salt, which, being hygroscopic, absorbs water thereby resulting in highly concentrated solutions of Cl⁻ at the saturation concentration around and underneath the salt deposit. Using the DFT adsorption energies, the constructed phase diagram shows that chloride can preferentially chemisorb onto the Fe(110) surface over a certain range of environmental conditions, thereby inhibiting the process of passivation on the metallic surface. Inhibiting the passive film formation can be expected to have implications for the ability of metals in the early stages of localized corrosion to repassivate and hence the transition from metastable to stable pit growth. From experimental investigations, Pletnev et al. (2000), found that chloride accelerates the corrosion of mild steel during active dissolution. Accordingly, Pletnev attributed this to the formation of a surface complex involved adsorbed chloride. Alvarez studied the influence of chloride on pitting potentials and speculated that the specific adsorption of chloride hindered pit repassivation and increased dissolution rate (Alvarez and Galvele 1984). These experimental investigations and interpretations are consistent with the results from these phase diagram studies based on the DFT computations.

7 Quantifying surface coverage: using DFT to generate adsorption isotherms

Whereas the surface Pourbaix and phase diagram approaches allow prediction of the transitions from one dominant surface phase to another, the coverage of a surface can also be mixed with multiple adsorbates having negative free energies of adsorption under any particular condition (Saraby-Reintjes 1985). In such cases, it is of interest to know quantitatively how much of each species is present and in what ratio. The simplest model to predict the surface coverage of a species in the environment onto the surface of a solid is the Langmuir isotherm (Gileadi 1987; Masel 1996). This model assumes that there is no significant interaction between adsorbed atoms or molecules on the solid surface, other than site-blocking (i.e., one molecular species is allowed per site, which can simulate repulsive interactions based upon how one selects the effective size of a ‘site’). For a reaction of the type:

where S represents the solid surface and A the adsorbate, then the equilibrium between surface and environment is given by the Langmuir equation:

here K_eq is the equilibrium constant, ΔG_ad is the free energy of adsorption, R is the ideal gas constant, and T is the absolute temperature. The activities are given by the terms, aSA, a_S, and a_A. In the Langmuir approximation the activities are given by the fraction of available surface sites populated by the species, θ, and the concentration of species A in the environment (partial pressure in the gas phase or the solution concentration, C_A, i.e., ideal solution behavior).

The noninteraction assumption inherent in the Langmuir approach may be considered reasonable if one considers that, at the limit where interactions start to be important, the coverage is ‘capped’ by the (usually) repulsive interactions: that is, the fractional surface site coverage term takes this effect implicitly into account through an effective site-blocking, where the ‘size’ of a site is determined by the limit of repulsive interactions between adsorbates (Conway et al. 1974). Although the Langmuir approach may not get the exact shape of the adsorption isotherm, it may be considered approximate enough to start making inferences about corrosion and hydrogen susceptibility.

To quantitatively evaluate this assumption for the case of Cl and O co-adsorption, Samin and Taylor performed some more accurate isotherm calculations using the numerical method of Grand Canonical Monte Carlo (GCMC) (Samin and Taylor 2018b). GCMC is a technique that allows the number of species in an atomistic/molecular simulation, usually on some kind of lattice that represents a crystalline solid or a surface, to vary to meet a given chemical potential (Blum et al. 2005; Byrd et al. 2000; Ceder 1993; Sanchez et al. 1984). Using this technique one can construct an adsorption isotherm without requiring an analytical expression, such as the Langmuir isotherm (Samin and Taylor 2017a; Schmidt et al. 2012). Instead, the surface coverage and arrangement of co-adsorbed species are optimized such that the chemical potential of the surface system obtained by the use of a model Hamiltonian is equal to that of the environment (Blum et al. 2005; Hart et al. 2005; Hegde and Bowen 2017). In the case of chemisorption, this can be achieved by using DFT to develop a database of interaction energies between adsorbed species (known as figures of merit). These interactions can be one-body (i.e., adsorption energy in the limit of no-interaction), two-body (pair-wise interactions), three-body, etc. (see Figure 10). The more interactions that are taken into account in the cluster expansion, the more accurate the isotherm will be.

Samin and Taylor (2018b) demonstrated this technique for Cl and O chemisorption on a Ni(111) surface. Once parameterized, the GCMC algorithm samples from a large ensemble of such configurations that are used to generate the adsorption isotherm as a function of the chemical potential. The many-body interactions produce plateaus in the total surface coverage versus chemical potential plot, indicating that significant chemical potential pressure must be exerted to go beyond a certain nominal coverage (i.e., repulsive interactions between the atoms on the surface effectively prevent them from becoming close neighbors until a high chemical potential is reached through extremely high partial pressures or high concentration solutions): example isotherms are shown in Figure 11. Ultimately the results from Samin and Taylor (2018b) indicated that when chemical potentials remain modest (i.e., typical solution chemistry and ambient temperatures) then the Langmuir approximation can be applied. In such cases, i.e., the left-hand side of Figure 11, the effective surface coverage of ∼0.3 ML (calculated in terms of the ratio of surface adsorbates to the total number of fcc binding sites) is equivalent to 1 ML in the Langmuir definition (100% occupancy of the total accessible surface sites).

8 DFT-informed Langmuir isotherms for Co-adsorption from aqueous solution

When there are multiple species that can interact with the surface, then there may exist a competition among those species for adsorption onto surface sites. In this case, the general form of equation (4) still applies but the total fractional surface coverage must be constrained by an equation of the form:

where j = M (empty site), Cl, O, OH, H₂O, etc.

When the adsorbates can vary considerably in size a scale factor may need to be taken into account since one adsorbate may occur >1 equivalent surface sites of one of the smaller species (Conway et al. 1974). Applying this constraint, and a series of equations of form Eq. (11) for each possible adsorption reaction, one can then simultaneously solve for the coverage of each adsorbate with a matrix equation. Alternatively, a Flory–Huggins isotherm may be used (Conway et al. 1974).

Taylor et al. (2018b) applied this Langmuir technique to consider scenarios for a series of commonly used structural base metals. The chemical potentials were determined from the pH and chloride concentration, and the electrochemical potential, and the adsorption equations were solved simultaneously to determine the coverage of each coadsorbate based on the free energy expression. The coverage of Cl, O, and H₂O on Ni(111) at various pH, potential coordinates is shown in Figure 12. The coverage is plotted against the electrochemical potential spectrum. The results are of a similar nature to what was seen in Figure 9 for Fe(110): at cathodic potentials water adsorbs (and in fact cathodically dissociates to also form adsorbed hydrogen, although that reaction was not simulated in Figure 12) (Taylor 2012b). As the potential is made more anodic, first the one-electron adsorption of chloride occurs, followed by displacement of the Cl by O. This behavior becomes more pronounced at acidic pH and higher chloride concentrations. The adsorption of chloride becomes less significant or is omitted altogether at higher pH and low chloride concentrations. These observations were interpreted against the scenario that occurs at the bottom of corrosion pits: low pH environments produced by hydrolysis of cations, and high chloride concentrations because of electromigration can lead to conditions which favor chloride adsorption on the exposed metallic surface, thus preventing repassivation (Anderko et al. 2004; Frankel et al. 2017; Sharland 1987). In the next section, the extension of this modeling approach to consider how alloying affects this process will be reviewed.

Figure 9:

Gibbs’ free energy of adsorbates on Fe(110) at pH 4.45 in a saturated solution of NH₄Cl and 298 K (Taylor 2012b). Reproduced with permission from NACE International (© 2012).

Figure 10:

Example of the configurations used to train the cluster expansion model for chloride-oxygen adsorption isotherms on Ni(111). For alloys, the number of configurations needed to train the alloy would need to be significantly expanded (Samin and Taylor 2018b). Reproduced with permission IOP publishing (© 2018).

Figure 11:

Isotherms for O and Cl adsorption as a function of the chemical potential of chloride. Typical chemical potentials for chloride ions in solution are from −7 to 0 eV (higher chloride concentration = higher chemical potential). Chemical potentials for oxygen range from −7 to −2 eV, depending on pH (higher pH = higher chemical potential) (Samin and Taylor 2018b). Reproduced with permission from IOP Publishing (© 2018).

Figure 12:

Ni(111) surface coverage versus the potential for different Cl molarities (1 mM, 10 mM, 0.1 M, and 1 M, respectively, from left to right) and pH (4, 7, and 10 from top to bottom, respectively) (Taylor et al. 2018b). W stands for water, O for adsorbed oxygen, and Cl for adsorbed chloride. Reproduced with permission from Elsevier (© 2018).

Isotherm models are frequently used in the interpretation of chemical corrosion inhibitor experiments, particularly those that measure inhibitor efficiency versus the concentration of the inhibitor (Christov and Popova 2004; Mondal and Taylor 2014). At the same time, there is some debate about the best way to interpret inhibitor efficiency data and linking it to adsorption phenomena (Kokalj 2010; Kovacevic et al. 2017; Walczak et al. 2019). Taylor et al. published a model for how a first-principles approach to inhibitors could be developed that couples adsorption behavior along with other phenomena that are important such as interaction with corrosion products, partitioning between oil and water phases, speciation, and ion-pair and multilayer formation (Taylor et al. 2015a, 2015c). This remains an area of much interest and such models have yet to see a full development (Kovačević and Kokalj 2013; Kovacevic et al. 2017; Peljhan and Kokalj 2011b). Coupling DFT with other atomistic simulation approaches for inhibitor film formation has seen renewed effort in recent years (Du et al. 2014; Khaled 2003, 2006, 2008, 2009, 2010; Khaled and Sherik 2013; Khaled et al. 2005).

9 Competitive electrochemisorption on alloy surfaces

The development presented so far only treats adsorption on a system where the adsorption energy is homogeneous across the metal surface. When an alloy material is used, however, there will be inhomogeneities across the surface, both at the atomistic level (i.e., in the statistical distribution of elements of the solid-solution alloy matrix) and at the microstructural level (grain-boundaries, secondary phases, intermetallics, etc.). For solid solution alloys, one can use purely statistically distributed (i.e., probabilities of occurrence based on composition) or otherwise informed distributions of the alloying elements and weigh the adsorption isotherm according to the expected distribution of surface sites with distinct adsorption energies. For the case of microstructurally heterogeneous materials, one could also construct adsorption isotherms for each phase and study them in their statistical distribution collectively or individually to understand the unique properties of each microstructural feature and their differentiating surface chemistries that may influence corrosion. For scenarios such as work-hardened zones which have local strain/dislocations, the adsorption energies will shift according to the extent of local strain; these shifts can also be determined from DFT calculations on strained surfaces (Greeley and Mavrikakis 2005; Lei et al. 2013). Francis and Curtin applied DFT to study the impact of localized surface stresses on the adsorption of environmental molecules for their relevance to electrocatalysis (Francis and Curtin 2015). An analogous effort could be applied to infer the impact of local stresses and strains on corrosion-relevant processes at the electrochemical interface.

As an example, consider the case of chloride adsorption on an alloy with unique sites expressed by the term M_i. If there are N unique sites each with their own corresponding adsorption energy and local environment, the possible reactions that could occur on the surfaces can be generally represented as:

For each unique reaction there is a corresponding equilibrium equation:

where there is a coverage term θxi for species x at the type of surface site i. The fractional surface coverage at each site i should sum to unity as before, and the total coverage across each type of site should be consistent with the statistical likelihood of each site i existing on the surface. This likelihood may be determined from a statistical assumption for an ideal solid solution or from external experimental information. The set of matrix equations is solved simultaneously as before and summed over the N sites of type i.

For most alloy materials, the number of unique surface adsorption sites that would need to be studied with DFT can get extremely large, making it computationally costly to run all combinations. To expedite the determination of adsorption energy at a site i the authors leveraged a machine learning approach by taking the results of ∼150 DFT calculations and then used logistic regression to estimate the adsorption energy from a sum of two-body terms between the metal atoms making up the adsorption site i (Ke and Taylor 2020b; Ke et al. 2020a, 2020b). This approach can be varied based on the geometry of the surface plane and the atomic make-up of the adsorption site. This approach was applied to O and Cl adsorption on close-packed alloy surfaces based on a stainless steel system (Ke et al. 2020a) and a Ni–22Cr system (Ke and Taylor 2020b; Ke et al. 2020b).

Figure 13 shows the performance of the local bonding model, which assumes that the adsorption energy can be obtained only as a sum of nearest-neighbor interactions. The scatter in the parity plot indicates that second-nearest neighbors and the overall matrix composition play a role in influencing the binding energy of both O and Cl to metal alloy surfaces. The order of significance of the elements can be assessed from Figures 13(c) and (d). For oxygen binding, the order from strongest oxygen binding to weakest is: Cr > Mo > W ∼ Mn > Fe > Co > Ru ∼ Ni > Cu. For chloride binding, the order is similar but Mn ∼ Fe > W. Thus, only W looks like it clearly differentiates between Cl and O in terms of binding affinity. Owing to the scatter in the fits, however, it was decided that more detailed analysis was necessary.

Figure 13:

Parity plots for the regression against (a) oxygen and (b) chlorine adsorption energies for the series of Ni–Cr alloys and stainless steel alloys. Contributions of alloy elements to (c) oxygen and (d) chloride to adsorption energies (eV) relative to the mean introduced into the Ni–22Cr alloy system.

The first focus of this approach was on studying the role of alloying in influencing the localized corrosion susceptibility of corrosion-resistant alloys (Ke et al. 2020b). The term “corrosion-resistant alloys” refers to the resistance of the alloys to uniform (or general) corrosion because of the presence of robust passive films across the bulk of the material surface. However, microstructural heterogeneities or discrete probabilistic events can create the risk of localized corrosion even in “corrosion-resistant alloys”: hence, many corrosion scientists have favored taking a stochastic approach to modeling localized corrosion (Organ et al. 2005; Stewart and Williams 1992; Riley et al. 1991; Williams et al. 1991). Factors known to enhance the corrosion resistance of an alloy are the adoption of alloying elements such as chromium and molybdenum, and/or the use of nickel-based rather than iron-based alloys (Marcus 1994; Sedriks 1996; Taylor et al. 2018c).

Following the general concept introduced in Figure 12, Ke and Taylor (2020b) sought to quantify the competition between chloride and oxygen for surface adsorption sites as a means to assess a material’s susceptibility versus resistance to localized corrosion within a givenvironment. The adsorption energies of Cl and O were computed for a series of model Ni–22Cr alloy slabs generated using the special quasirandom structure model (Zunger et al. 1990), which is a method for approximating in a periodic model the atomic correlations that exist in truly random solid solutions. The series consisted of the special quasirandom Ni–22Cr slab, generated earlier by Samin and Taylor (Samin and Taylor 2018c), with various elements replacing either Ni or Cr atoms on the surface layer of the slab, at which site the O or Cl atoms were chemisorbed. The chemisorption interaction of Cl and O with the surface is a mixture of covalent and ionic chemical bonds. The possible alloying elements considered were taken from the set: Cr, Mn, Fe, Co, Ni, Cu, Mo, Ru, W.

Because O and Cl both adsorb into three-fold hollow sites on the metallic surface, there are numerous possible combinations of those three binding metal atoms that could coordinate to the O or Cl atom. i.e., for the nickel-chrome binary system, chemically distinct adsorption sites can include Ni–Ni–Ni–O, Ni–Ni–Cr–O, Ni–Cr–Cr–O, etc. It is also likely that second nearest-neighbor effects may play a role in the adsorption environment. The plot of total chlorine versus oxygen adsorption energies for the alloys considered is shown in Figure 14. On the left are all the data points, and on the right are the data points plotted for the case in which the adsorption site has all three of the adsorption sites in common (e.g., Fe–Fe–Fe or Cu–Cu–Cu).

Figure 14:

Adsorption energies for Cl versus O on a series of Ni–22Cr alloys with dopant elements added into the surface adsorption sites (Ke and Taylor 2020b; Ke et al. 2020b). The solid squares represent those surface sites which consist of all three elements bonded to Cl or O having the identity as labeled. The remaining data points (blue circles) are composed of mixtures of those elements. Reproduced with permission from IOP Publishing (© 2020).

As seen in Figure 14, generally the chloride and oxide adsorption energies scale with one another. Therefore, from this basis alone it is not clear to determine scenarios in which certain alloying elements or their concentration in the alloy may contribute to corrosion resistance. Instead, the competitive adsorption isotherm model was extended to the alloy surfaces to create surface coverage plots. The coverage plots for a Ni–Cr binary system with water, chloride, and oxygen adsorption as a function of the electrochemical potential (the critical controlling parameter in aqueous corrosion) are shown in Figure 15. The data are shown for an assumed condition of pH 4 (mildly acidic) and high chloride concentration (1 M). Apparent in these plots is the presence of three electrochemically controlled phases on the Ni–Cr binary surface, along with regions of co-adsorption. On the cathodic side, water adsorption is the dominant surface phase (potentially hydrogen adsorption from proton or water reduction as well, although that was not included in this analysis, see the next section). As the potential becomes more anodic, a chloride adsorption phase displaces the water, which we can associate with more active corrosion. At more anodic potentials the oxygen adsorption begins to dominate, as this is a two-electron adsorption process, and will ultimately lead to passivation as multilayer adsorption and formation of the oxide occurs. The multilayer formation of oxide is not able to be addressed within this competitive adsorption framework. Further investigations with multilayer isotherms developed using DFT studies of oxide/metal interfaces like that of Choudhury et al. (2014) would be required in the future. Modeling solute capture in the formation of passive oxide films is also a new frontier that is starting to see interest from surface science measurements and theoretical studies (Sherman et al. 2019; Yu et al. 2018).

Figure 15:

Ni–Cr binary alloys of up to 25% Cr composition with the coverage of unique surface phases plotted as a function of the electrochemical potential, which is the critical factor controlling aqueous corrosion. Because of the thermodynamic equilibria, no OH adsorbed phases are present in appreciable quantities for the range of potentials and conditions plotted. The conditions used to solve the Langmuir isotherm model are 300 K, pH 4, and 1 M Cl⁻ concentration (Ke et al. 2020b). Reproduced with permission from IOP Publishing (© 2020).

As seen in Figure 15, the effect of mixing in Cr adds additional structure to the electrochemical phase diagrams. From one perspective, the scale of electrochemical potential can be considered as providing an energy spectrum for the corrosion activity of the alloy surface. Chromium decreases the ability of chloride to totally cover the surface during the electrochemical potentials associated with active corrosion, thus decreasing the corrosion currents possible during localized corrosion. For example, at 25% Cr, there is only a very small range of potentials where Cl is the dominant surface phase, and even in this case the total coverage is relatively small.

For the purposes of understanding the role of additional elements in changing the corrosion susceptibility of an alloy to corrosion, Ke et al. constructed an index that aggregates the information in Figure 15 by integrating the total chloride adsorption across the electrochemical energy spectrum. This quantity was called the chloride susceptibility index (CSI) (Ke et al. 2020b):

Formally, since the surface coverage is a unitless area fraction quantity, the CSI will have units of V.

In Figure 16, Ke et al. (2020b) illustrated the use of the CSI to infer the localized corrosion susceptibility of Ni–Cr–X ternary alloys with substitution of Ni by different elements X. The model shows that Fe and Cu make the alloy more susceptible to chloride interference with passivity (and hence, decreasing the tendency for repassivation), whereas Co, Mn, Ru, Mo, and W, in increasing order of significance, all have positive effects. Of course, there are many other influences these variables can play in affecting corrosion mechanisms; such as the ability to grow sufficiently protective oxide films, their resistance to breakdown, and metallurgical effects that might influence surface segregation, microstructure, and diffusion (Clayton and Olefjord 2011; Grabke et al. 2004; Greeley and Nørskov 2007a; Han et al. 2015; Kirchheim et al. 1989; Liu et al. 2015b; Ma et al. 2017; McCafferty 2000; Punckt et al. 2004; Scully 2009; Shibata 2007; Tang et al. 2014; Qiu et al. 2017).

Figure 16:

Chloride susceptibility index plotted for a Ni_90−xCr₁₀A_x for various potential alloy elements A = [Cu, Fe, Co, Mn, Ru, Mo, W] that could be selected to change the alloy’s net resistance to localized corrosion by integrating the surface coverage of chloride over the relevant electrochemical energy spectrum (Ke et al. 2020b). Reproduced with permission from IOP Publishing (© 2020).

The CSI was applied to the Ni–Cr and stainless-steel systems and correlated against the repassivation potential obtained from electrochemical measurements (Figure 17) (Ke and Taylor 2020b; Ke et al. 2020a, 2020b). The correlation bears out the general theory proposed by Anderko and Sridhar in which a thermodynamic-kinetic model for repassivation versus salt film formation was proposed as a criterion for localized corrosion pit growth (Anderko et al. 2004, 2006, 2008, 2014). Their model contained adsorption constants and mass-transport data that was fitted to experimental repassivation potential data and polarization curves. The model used in Figure 17, on the other hand, uses only first-principles information but is seen to lead to a similar interpretation and validity relative to the experimental findings.

Figure 17:

CSI plotted against the experimentally determined repassivation potentials, E_rp: (a) Ni–Cr–Mo alloys (Ke et al. 2020b); and (b) stainless steels (Ke et al. 2020a). Reproduced with permission from IOP Publishing (© 2020).

Figure 18:

First-principles Evans diagrams from Luo et al. (Left) Schematic showing the principle of superposition of anodic and cathodic branches, and (right) First-principles polarization curves of pure Mg, Mg–0.3 wt%Ge (labeled as Mg–0.3Ge), and Mg–0.1 wt%Cu (labeled as Mg–0.1Cu) in saturated Mg(OH)₂ solution (Luo et al. 2019). Reproduced with permission from IOP Publishing (© 2018).

The chloride susceptibility index can potentially be generalized as a corrosion susceptibility index. On-going work is assessing other potential adsorption phenomena of interest, and the atomistic intermediate states involved in the dissolution of metal atoms as cations to come up with additional features for indexing the susceptibility of a surface to corrosion processes (Ke and Taylor 2020a).

10 Simulation of anodic dissolution and cathodic reduction kinetics from DFT

Anodic dissolution is one of the quintessential processes that come to mind when thinking of aqueous corrosion: the dissipation of metal atoms as ions into solution, along with the formation of oxide scales and rust (Waldman 2015). Discovering opportunities to manipulate materials’ features and composition to mitigate the rate of anodic dissolution using a first-principles approach could have impacts on the integrated computational materials engineering (ICME) of a new generation of corrosion-resistant materials (Taylor et al. 2018c). Following the idea of Marcus, alloying elements can be plotted in a two-dimensional space where the axes are indicative of the metal-metal bond strength (Marcus used cohesive energy) and the metal-oxygen bond strength (Marcus used adsorption energy of oxygen to the surface) (Marcus 1994). The concept is that metals with high metal-metal bond strength will be dissolution moderators, as they will have higher resistance to the reaction which leads to dissolution and/or migration of the metal atoms from the metallic phase into a corrosion product oxide film. At the same time, metals with a high affinity for oxygen will promote the formation of strong robust oxide films. Windl and co-workers recently used DFT data to expand upon this idea by studying the effect of dissimilar elements on metal-metal bond energies (Oberdorfer and Windl 2019). Their results showed that elements with high cohesive energy can, effectively, propagate that high cohesive energy to dissimilar elements when they are alloyed together. In soon-to-be-published work, they have made similar connections for the impact of dissimilar elements on the metal-oxide bonding energy (Chen and Windl 2020). This approach will allow expansion of the Marcus approach, which considers only elements in isolation, to alloy environments.

Simulating the dissolution process using DFT is a challenging undertaking. Greeley et al. evaluated adatom adsorption energies for a series of alloys as a means for estimating dissolution and deposition potentials (Greeley and Nørskov 2007a; Greeley 2010). Taylor examined the metal-metal bond-breaking process that occurs during dissolution by simulating the transition of a metal adatom from a Cu surface into its final state as a solvated ion (Taylor 2009). This work, which is currently being extended to Ni and Ni–Cr alloy systems by Ke and Taylor (2020a), should allow the first-principles estimation of the activation energy for dissolution (i.e., something akin to a Dissolution Resistance Index), and hence allow theorists to go beyond the thermodynamic predictions of the preceding efforts to first-principles models for electrochemical kinetics.

A semi-empirical approach for electrochemical kinetics, which borrows from first-principles calculations but also draws heavily from empirical parameters to scale the model predictions to experiment, has been applied in recent years to the efforts to understand corrosion kinetics on magnesium and magnesium alloys (Ma et al. 2017, 2019, 2020; Liu et al. 2011; Lü et al. 2011). Ma et al. (2017) propose that anodic dissolution kinetics can be derived using the following set of equations:

These equations provide a first-principles basis for deriving the anodic dissolution exchange current density, given a knowledge of the work function of the material, and the energy to create defect states; namely the free surface, vacancies at the surface, and adatoms. In a way, this equation is an analytical, derivable form of the numerical off-lattice kinetic Monte Carlo approach developed by Taylor using the embedded atom method to directly calculate the activation energy of metal atoms at different defect states on a metallic surface (Taylor and Liu 2013), which was in turn inspired from the on-lattice version assumed by Erlebacher and Newman (Artymowicz et al. 2009; Erlebacher 2004). The work function and surface energy may be directly determined from DFT calculations in Ma et al.’s approach (Ma et al. 2017, 2019, 2020), whereas other terms remain fitting parameters because of the challenge in performing DFT calculations sophisticated enough for their accurate determination.

The cathodic reduction kinetics have also been derived based on oxygen reduction or hydrogen reduction, using first-principles data. Taylor investigated the anomalous hydrogen evolution kinetics on Mg (Atrens and Dietzel 2007; Bender et al. 2011; Frankel et al. 2013; Williams et al. 2013) using a variety of DFT calculations performed in the literature and a proposed mechanism that involved Tafel adsorption, site-blocking by OH, and elimination of MgOH⁺ by dissolution at anodic potentials (Taylor 2016). This last reaction was proposed to account for continued HER at positive potentials by creating new sites for availability for hydrogen adsorption and recombination by removing the OH site “poisons” through an anodic mechanism. Yuwono extended this idea by adding in the Heyrovsky mechanism and parameterizing a full kinetic model using DFT that accounts for much of the observed anomalous hydrogen evolution behavior in Mg and Mg alloy materials (Yuwono et al. 2019).

Corrosion kinetics are traditionally expressed through the mixed potential theory, in which the corrosion current density (i.e., corrosion kinetics, usually on a log scale) is plotted against the electrochemical potential (i.e., the thermodynamic driving force) (Muralidharan 2002). This approach, developed by Evans and Hoar in 1932 and subsequent works (Evans and Hoar 1932, 1934), provides a useful way to consider the contribution of relative cathodic and anodic reactions to the overall current density and corrosion performance of a metal or alloy in an environment, and the effect of various inhibition and corrosion mitigation strategies (Frankel 2016). Reconstruction of these “Evans diagrams” for real and hypothetical Mg alloys using the first-principles basis has been published in recent work (Luo et al. 2019). The approach, similar to the above uses a combination of first-principles DFT evaluated data for work functions and surface adsorption properties, alongside empirically fitted scale factors to try to forecast the effect of different alloying choices on the final electrochemical kinetics of these materials (Luo et al. 2018, 2019). Luo et al. proposed that a semi-empirical implementation of the first-principles electrochemical kinetics using the Butler–Volmer kinetics can be applied to evaluate the relative corrosion behavior of Mg–Ge and Mg–Cu alloys (Luo et al. 2018, 2019). The anodic branch was fitted from the experimental data and first-principles predictions of the anodic equilibrium potential, and the cathodic branch was a combination of fitting parameters and activation energies for HER estimated by evaluating the energy of the intermediate adsorbed hydrogen state for the different alloy combinations and relevant intermetallic phases (Luo et al. 2019).

11 Electrochemisorption effects on hydrogen assisted crack growth

When a crack develops in a material that is exposed to an electrochemical environment, processes similar to those occurring during initiation of localized corrosion can activate (Figure 19): freshly exposed surfaces can either repassivate or remain bare for metal dissolution and hydrogen uptake depending on pH, chemistry, and electrochemical potential (Figure 18). The surface adsorption of oxygen can be a precursor to passive film formation. On the other hand, water or hydrogen adsorption can create a source for surface hydrogen that may then diffuse into the material, thereby potentially accelerating crack growth because of hydrogen modifying the plasticity in the fracture process zone or weakening the atomic planes within the lattice or at grain boundaries, thus lowering the resistance to decohesion (Barrera et al. 2018; Van der Ven and Ceder 2004). Other mechanisms for hydrogen-assisted fracture might include enhanced vacancy formation and coalescence because of H-vacancy complex formation (Fukai 2002), or the formation of brittle hydride phases, depending on the material.

Figure 19:

Schematic showing processes that can occur during stress corrosion cracking and/or corrosion fatigue; lowering electrochemical potential from passive to active region, metal dissolution, and water reduction leading to high pH and hydrogen generation which can cause embrittlement. Reproduced with permission from Ramgopal Thodla (2020).

Scully and Moran (1984) and Frankenthal and Milner (Quinlan 1986) conjectured that water reduction reaction on the exposed bare metal surface created on a surface under strain may produce a steady-state concentration of atomic hydrogen on the surface when the reduction reaction occurs via the pathway:

The driving force to create an equilibrium between H_ads and H_abs will create a boundary condition for the flux of hydrogen into the material, thus providing one impetus for further hydrogen embrittlement of the material via its entry to the fracture process zone.

Environmentally assisted cracking, including stress corrosion cracking and corrosion fatigue, is commonly described as a failure mode that involves the intersection of (a) a susceptible material, with (b) a corrosive environment, and (c) a mechanical driving force (Ford 1984). The flow chart in Figure 20 outlines some of the various kinds of environmentally assisted cracking and the processes and sub-processes that lead to crack initiation and growth in the susceptible material. Although this is not the place for a detailed review of environmentally assisted cracking (see a recent review of hydrogen embrittlement of steels, for example, Barrera et al. (2018)), examination of this flow chart indicates that surface processes may play decisive roles in limiting or enabling environmentally assisted cracking because of the potential significance for mechanisms such as dissolution at the crack tip, water and hydrogen/hydronium adsorption, repassivation, transport, and complexation.

Figure 20:

Process and sub-processes that may be active and rate-controlling in corrosion fatigue and stress corrosion. Adapted from Ford (1984).

In a combined theoretical/experimental study, crack growth rate measurements were made under different testing conditions on pipeline steels as a function of the applied cathodic protection potential by Thodla et al. (Thodla et al. 2020b). The crack growth rate trend as a function of potential was compared to the DFT predictions of the potential-dependent surface chemistry. The authors studied the crack growth rate in three grades of pipeline steel that differed in strength and under both static and dynamic loadings. In all three cases, when the potential was scanned from −1050 mV SCE to −800 mV SCE, a drop of between 1–3x orders of magnitude was observed in the crack growth rate, corresponding directly with the potentials at which the coverage of water and atomic hydrogen on the surface drops from 1.0 to 0 monolayers (Figure 21). According to Wei et al. (Lu et al. 1981; Wei 1980; Wei and Gangloff 1989; Wei et al. 1980), there are three possible rate-determining steps for hydrogen-assisted crack growth:

• – Transport of species to the crack tip

• – Reaction of species at the surface

• – Diffusion of hydrogen to the fracture process zone

Figure 21:

Comparison of the crack growth rate measured at pH 8.2 and water adsorption on Fe110 surface (at pH = 11, which are the expected conditions at the crack tip over a range of cathodic potentials) as a function of applied potential: (a) a 904 MPa strength low alloy steel (sample ID 2683) at 0.1 mHz, ΔK = 22 MPa, and R = 0.6; (b) the 904 MPa strength low alloy steel at various values of K (i.e., limiting static crack growth rate); and (c) an 1166 MPa strength low alloy steel for various values of K (Thodla et al. 2020b).

However, other possible rate-limiting processes may be considered, such as (Thodla et al. 2020b):

• – Development of local strains/strain rates high enough to interact with the hydrogen

• – Hydrogen interaction with local strain to cause crack propagation

Given the correlation between the predicted change in surface chemistry as a function of electrochemical potential, the reduction in crack growth rate (CGR) with decreasing magnitude of the applied cathodic potential on low alloy steels seems to indicate that surface processes become rate-limiting around −900 mV SCE, and below those potentials, microstructural factors appear to play a more significant role. Figure 21 shows the comparison.

The development of the DFT-based surface coverage model for alloyed materials follows the same approach to be taken for studying the electrochemical potential-dependent behavior of the hydrogen-assisted crack growth for more complex materials, such as corrosion-resistant alloys required for extreme conditions in the oil and gas industry. A study of the hydrogen embrittlement of Inconel alloy 718 under cathodic polarization showed that, such as pipeline steels, hydrogen-assisted crack growth rate can significantly increase once a threshold in the applied cathodic electrochemical potential is crossed (Thodla et al. 2020a). Potentiostatic scratch measurements indicated that the transition from cathodic current transients to anodic transients occurred at the potential associated with the transition from predominantly water-covered metal surface to oxide-covered surface, as predicted by the alloy adsorption isotherm model (Figure 22(a)). Likewise, around this same transition, the crack growth rates associated with hydrogen assisted crack growth increased by several orders of magnitude (Figure 22(b)).

Figure 22:

(a) The transitional range for the coverage of water versus oxygen as a function of applied potential to the bare metal surface of alloy Inconel 718 aligns with the transition region between predominantly cathodic transients to anodic transients produced from scratching experiments. (b) Likewise, the transition to higher water coverage/lower oxygen coverage on the Inconel 718 alloy surface induced by cathodic polarization correlates with the increase in crack growth rates on the two heats (120 and 150 ksi) studied for Inconel 718 (Thodla et al. 2020a).

The crack growth characteristics of Inconel alloy 725 under cathodic polarization were studied by Ramgopal et al. (Thodla et al. 2021) and the surface chemisorption model applied to this alloy as a third case study for testing this computational modeling approach. A series of constant K tests on different heats and for different K values showed that crack growth rate associate with hydrogen embrittlement sharply increased by 1–2 orders of magnitude from −850 to −950 mV SCE and then more gradually increased from −900 to −1050 mV SCE. As in the former studies, it was observed in this joint computational/experimental work that the decrease in crack growth rates at less cathodic potentials correlated with the decrease in water coverage that occurs when oxide instead begins to cover the surface sites on the alloy. Based on these insights, a crack-tip strain rate model was developed to account for the changes in the crack growth rate as a function of both the environmental and loading conditions. Figure 23 shows the correlations between CGR and surface coverage for the various heats of 725 studied.

Figure 23:

A comparison between the crack growth rate of two heats of 725 (H2 and H4) as a function of water coverage with that of 718 from comparable values of K (∼80 and 100 MPa√m) (Thodla et al. 2021).

In addition to allowing simulations of the potential surface interactions that might occur between a material and the environment that are relevant to hydrogen embrittlement, DFT simulations have also been performed to study some of the details of how hydrogen affects the microstructure. These include the interactions between hydrogen and the metallic lattice (Teus et al. 2007), including diffusion (Chohan et al. 2018), hydrogen, and the decohesion energy of crystallographic planes or grain boundaries (Geng et al. 2005; Tehranchi and Curtin 2017; Tehranchi et al. 2020; Van der Ven and Ceder 2004), hydrogen–vacancy interactions (Lu and Kaxiras 2005; Nazarov et al. 2010; Olsson et al. 2015; Tateyama and Ohno 2003; Taylor et al. 2013, 2014), hydrogen effects at the interface between cementite and the iron matrix (McEniry et al. 2018), formation of hydride phases (Blomqvist et al. 2011; Lillard et al. 2014; Udagawa et al. 2010), and hydrogen interacting with stacking faults and localized plasticity (He et al. 2019; Juan et al. 2002; Lu et al. 2001, 2002; Zamora et al. 2012). These kinds of studies continue to fill out the fundamental knowledge base regarding the influence of hydrogen on materials behavior and will, in time, provide the basis needed for developing fully first-principles informed models of hydrogen-assisted cracking.

12 Passive film breakdown studies

Adsorption reactions on low-index metal surfaces are relatively straightforward calculations compared to those required to study reactions that might take place on the surfaces of oxide materials. This challenge notwithstanding, there remains interested in understanding how oxide films become compromised because of mechanical or chemical interactions, and, in particular, the role of halides in initiating localized corrosion. Efforts that applied DFT in this area in the recent literature include studying the role of chloride in changing the surface stability and vacancy formation mechanisms on passive aluminum oxide surfaces (Marks 2018; Yu et al. 2019), and the displacement of surface OH in hydroxylated nickel oxide surfaces with chloride (Bouzoubaa et al. 2009a, 2010). There still remains to be created, however, an interpretative mechanistic framework that connects these atomistic steps with predictive thermodynamic or kinetic models for passive film breakdown and connection to the electrochemically measured parameters such as the breakdown or repassivation potentials, although one such framework has been initiated that breaks down the various phenomenological contributing factors into kinetic triggers and thermodynamic driving forces (Yu and Marks 2018). We forego a more detailed review of the structure, properties, and reactivity of oxide surfaces surveyed using DFT at this time, and anticipate that advances over the coming years will mature the field toward this direction of integration of DFT insights into phenomenological, electrochemical models for passive film breakdown.

13 Perspectives

As this review indicates, DFT calculations have come a significant way in describing elements of the electrochemical interface relevant to corrosion and hydrogen embrittlement phenomena. There still remains much to be carried out, however, to improve the accuracy and credibility of these models, to expand into new domains of materials corrosion and degradation phenomena, and to develop both larger datasets of relevant parameters computed by first-principles and publicly available algorithms for making those datasets and models usable to the larger corrosion science and engineering community.