Phase field modeling for the morphological and microstructural evolution of metallic materials under environmental attack

The complex degradation of metallic materials in aggressive environments can result in morphological and microstructural changes. The phase-field (PF) method is an effective computational approach to understanding and predicting the morphology, phase change and/or transformation of materials. PF models are based on conserved and non-conserved field variables that represent each phase as a function of space and time coupled with time-dependent equations that describe the mechanisms. This report summarizes progress in the PF modeling of degradation of metallic materials in aqueous corrosion, hydrogen-assisted cracking, high-temperature metal oxidation in the gas phase and porous structure evolution with insights to future applications.


INTRODUCTION
Corrosion amounts to a staggering 3.4% 1 and 3.34% 2 of the gross domestic product in the USA and China, respectively; in the USA, this cost greatly exceeds the costs arising from all-natural disasters combined 1,3 . The degradation of metallic materials in natural environments is difficult to avoid due to the practical applications of these materials in corrosive environments. A metal surface is usually covered by a passive film that was formed during service. When this passive film is damaged or in poor quality, the metal is exposed to the environment, resulting in a localized or nonlocalized attack. Of these types, the localized attack is the most dangerous because it is difficult to predict and detect which can result in sudden structural failure 4 .
Well-designed experiments can be used to evaluate key parameters affecting and controlling the degradation of a material of specific composition in a given environment. Theoretical work can improve understanding of the observed phenomena, which is important in the design and development of tools for service life prediction and material protection. There are two key areas of theoretical research. Thermodynamics addresses whether a degradation process is feasible; it also provides possible stable phases and equilibrium properties of the system. Kinetics explores how fast a degradation process can proceed under given environmental conditions.
Theoretical modeling of degradation kinetics can be classified into three levels by a length scale. Macroscale theories predict the average behavior of systems over a millimeter to meter length. Mesoscale models predict the effects of the local microstructure and environment on charged particle transport, interface stability, phase transition, topological evolution, and kinetics at the material surface on a scale of micrometers to millimeters, although some can even resolve nanometer-scale phenomena. Atomistic models track the interactions and motion of individual atoms/molecules or the mechanisms of bonding and/or charge transfer at interfaces, enabling an understanding of degradation mechanisms at the atomic level and the accurate estimation of some thermodynamic and kinetic properties (e.g., formation energy, binding energy, and energy barriers to diffusion and reaction). Atomistic modeling cannot encompass the length and time scales over which degradation takes place in engineering applications. On the other hand, the macroscale approach makes it difficult to capture key features, such as the effects of microstructural features and defects on the initiation of degradation and the kinetics of morphological and microstructural evolution during the localized attacks. The majority of theoretical modeling research are in mesoscale.
The challenges encountered by most kinetic models proposed to date include interface movement, the electrical double layer at the metal-electrolyte interface and the evolution and effect of insoluble corrosion products (ICPs). Several early steady-5-10 and transient-state 11-20 models do not allow for changes in the shapes and dimensions of the metal as corrosion proceeds, and most attempt to solve a series of transport equations (modified Fick's law, also called the Nernst-Planck equation) in an electrolyte under concentration and electrostatic potential gradients while considering the chemical reaction rates in the system. In an early effort, Sharland and Taskeer 5 proposed a one-dimensional (1-D) steady-state model that solves a set of Nernst-Planck equations coupled with a neutrality condition in the electrolyte to obtain molar concentration and electrostatic potential profiles. This model considers the metal-electrolyte interface as a fixed boundary and assumes a Neumann boundary condition that follows Tafel kinetics. It can predict the metal dissolution rate and determine the relevant concentration and electric potential profiles. Galvele 21 , Turnbul and Thomas 22 , Walton 10 , and Krawiec et al. 9 proposed similar models that also solve a set of mass (or ion) transport equations with a steady-state assumption by ignoring interface movement. lifespans of materials. Several numerical models of metallic materials degradation have been developed based on these techniques. In the following, we provide a brief review of the important features and limitations described in prominent publications about the above-mentioned methods.
Xiao and Chaudhuri 23 developed a FEM-based model that solves a set of mass transport equations coupled with the Poisson equation or neutrality condition. This model considers the damaging surface as a sharp interface and requires a matching mesh at each time step. Although this strategy is promising, the sharp-interface assumption makes it difficult to implement in complex 2-D or 3-D problems. Several other numerical sharpinterface models that use FEM to solve the coupled equations include those proposed by Oltra et al. 24 and Sarkar et al. 25 for uniform metal corrosion and those proposed by Kota et al. and Brewick et al. for microstructure-sensitive corrosion 26,27 .
Several attempts have been made to model the degradation of steel and aluminum alloys using the FVM. In FVM models, a phasechange strategy is introduced by evolving the matching mesh based on the concentration variation of the ionic species [28][29][30] . These models address the matching mesh problem to some extent but continue to assume a sharp metal-electrolyte interface, which may reduce the accuracy of simulation results for complex 3-D problems.
CA is an advanced technique that has been used to model multiphase problems. In this approach, the transition of one phase to another depends on the predefined probability factor, which is related to the chemistry of the corrosion problems. CA models do not solve any partial differential equations (PDEs) and therefore have a low computational cost. Additionally, CA corrosion models can capture large-scale effects. These models have been used to study 2-D and 3-D geometries. Stafiej et al. used each cell to represent solid (metal, passive film, and reactive metal) and electrolyte (acidic, basic, and neutral) 31 . Several other CA models available in the literature for modeling corrosion include those proposed by Di Caprio et al. 32,33 , Pérez-Brokate et al. 34 , and Van der Weeën et al. 35 . CA modeling has also been used to address intergranular corrosion (IGC). For example, the models proposed by Lishchuk et al. 36 and Di Caprio et al. 33 focus on IGC in sensitized alloys. Although CA models have a low computational cost, as mentioned above, the lack of time-dependent dynamics makes it difficult to predict a corrosion event quantitatively.
ALE is a moving-boundary method that has been used to model corrosion problems. This approach also assumes a sharp metal-electrolyte interface. The models proposed by Yin et al. 37 and Wang et al. 38 include the formation of corrosion products. Yin et al. 37 and Wang et al. 38 also predicted a key role of the porosity of these products in corrosion kinetics and reported that the corrosion rate decreasesd significantly in the presence of corrosion products under the same environmental conditions. ALE is based on a moving mesh method that significantly increases the computational cost due to matching mesh generation at each time step when modeling non-uniform or non-regular geometry problems. According to the authors 39 , the mesh quality decreases over time even when modeling a 2-D geometry which requires specific mesh quality-based criteria to stop the simulation and redraw the geometry. PD has been used recently to model corrosion problems. PD models discretize the domain using a simple structured grid and track the interface location implicitly based on the concentration profiles. These models use the concept of effective diffusivity to control the rate of solid-phase dissolution. Bobaru and colleagues 40,41 initially proposed PD models for localized corrosion processes and recently developed modified PD models for stress corrosion cracking and IGC cracking in metallic materials [42][43][44][45] .
Unlike other computational techniques, the PF method is built by defining the diffuse interface (with finite thickness) between different phases, and the interface evolves due to the minimization of the total free energy of the system. FVM, PD and PF employ PDEs that govern the interface movement with time. However, FVM and PD focus only on solving concentrationdependent diffusion equations while PF employs an additional Ginzburg-Landau equation that governs the morphology evolution. CA, PD, and PF have shown good capability and promise in studying complex evolving morphologies where other methods have shown limitations. CA and PD are efficient methods to model the phase change process with far less computational cost as compared to PF method. CA method does not solve PDEs while PD method uses a meshfree technique that reduces the computational cost. However, PF is the only method that is built on the total free energy minimization of the system (which is the most fundamental way to address the phase change or microstructure or morphology evolution problem). Recently, Jafarzadeh et al. 42 compared the advantages and disadvantages of the FVM, PD, CA, and PF methods in a recent review article. In contrast, our review article focuses on the applications, benefits, and limitations of the PF method. Therefore, we refer readers to the review article by Jafarzadeh et al. 42 for more details and comparative analysis.
PF models are built on the principle of treating the evolving interface as a diffusive interface with finite thickness instead of using a sharp one by introducing smooth and continuous mathematical functions. Dimensionless phase-field variables (or order parameters) are introduced that represent the diffusive interface. These variables evolve due to the minimization of the total free energy of the system 46 and do not require explicit tracking. Consequently, PF models are more suitable for solving complex 3-D geometry problems due to their inherent numerical stability. Over the years, The PF method has been applied to model several material processes. The applications of this method range from solidification, dendrite growth, phase transformation, electrochemical corrosion & deposition, dislocation dynamics, crack nucleation and propagation, void nucleation and migration, gas bubble evolution, and intercalation process in battery materials 46,47 .
Thermodynamically, phase transformations, chemical reactions, and microstructure evolution are driven by free energy minimization in the system. In a pioneering study in the early 1970s, Asaro and Tiller 48 considered the contributions of strain energy and surface distortion to free energy in the context of stress corrosion cracking (SCC) via surface diffusion. That study focused on the criteria of fracture initiation rather than on the evolution of morphology. PF models make full use of fundamental principles that govern the phenomena and derive phase-field equations coupled with relevant transport equations to predict the system morphology. These models describe multiphase systems using a set of field variables (conserved and non-conserved) that have smooth transition across the evolving interfaces. The Cahn-Hilliard 49 equation is a non-linear diffusion equation while the Ginzburg-Landau equation is the relaxation equation (also known as Allen-Cahn equation) 50 governs the evolution (temporal and spatial) of the phase-field variables. Furthermore, by utilizing the kinetics and thermodynamics of the system under study as the input, the PF method can show great potential in modeling the complex microstructural and morphological evolution of the material without the need to explicitly track the interface 46 . As a mesoscale method, PF modeling has made significant advances over the last two decades 46,47,51,52 and has been used to predict many material processes, as mentioned above. This review focuses on progress in the use of the PF method to model morphological evolution in metallic materials under aggressive environmental and electrochemical conditions.

TYPICAL CHEMICAL REACTIONS ON THE METAL SURFACE
The mechanisms that lead to the morphology and microstructure evolution when metals are under environmental attack can vary T.Q. Ansari et al. significantly. The aqueous attack usually involves electrochemical reactions, while the gas phase attack usually only involves chemical reactions on the metal surface. However, irradiation will involve damage not only on the surface but also in the bulk of the metallic materials. HAC involves hydrogen ingress from the surface, usually initiated by metal corrosion. Hydrogen diffuses to defects inside the metal such as microcracks or voids and weakens the metallic bonds to cause crack growth at low-stress levels. Stress corrosion cracking is a brittle failure that requires three conditions present simultaneously: a corrosive environment, an SCC susceptible metal, and some component of tensile or shear stresses. The porous structure formation during dealloying is a sub-category of metal corrosion where less noble metal will be removed by corrosion from an alloy surface in an aqueous environment. It is evident that for most forms of environmental attack, except for irradiation, chemical reaction at the metal surface is involved. In the following, we discuss some typical reactions related to these topics in more detail.
In the case of metal corrosion, as noted earlier, metal surfaces are usually covered with a passive film that protects against oxidation in the presence of a corrosive electrolyte. Areas of damage to this film can act as nucleation sites for localized corrosion. The corroding metal surface is the anode. The chemical reactions on anode are oxidation reactions, i.e., the metal loses electrons to become metal ions. A key difference between metal aqueous corrosion and gas-phase oxidation is that, in the corrosion process, there must be a cathode to collect electrons lost by metal corrosion through an electric connection. The cathode can be the same metal or another conductor connected to the corroding metal. The chemical reactions on the cathode surface are reduction reactions in which some cations obtain electrons from the cathode surface. Several attempts have been made to understand the surface science and corrosion initiation mechanism using advanced experimental techniques at the atomic and nanometric scale, as recently reviewed by Maurice and Marcus 53 . Considerable advances in knowledge of reaction kinetics at the metal-electrolyte interface (at the atomic scale) and the understanding of early-stage corrosion mechanisms have been made due to ultra-high vacuum (UHV) studies based on scanning tunneling microscopy (STM) coupled with density functional theory (DFT). Readers are referred to the review article by Maurice and Marcus for a detailed overview 53 . Once the corrosion process is initiated, the metal oxidation reaction at the metal-electrolyte interface can be expressed as In an early example, Bockris et al. 54 proposed a mechanism for iron corrosion in aqueous environment that has been widely accepted in the corrosion science community. This mechanism includes the effect of pH (or H + ion concentration) on the corrosion rate and pit morphology in iron (Fe), as follows: where subscripts (s), (aq), and (ads) represent the solid, aqueous, and adsorbed phases, respectively. In this example, Fe atoms react with hydroxide ions to form adsorbed iron hydroxide, which then oxidizes to form iron ions (Fe 2+ ) through a one-electron transfer step. During pitting corrosion in salt water, several corrosion products, such as FeCl 2 , Fe(OH) 2 and Fe(OH) 3 , may form as follows: where the subscript (ppt) indicates a precipitate. In a corrosion cell, the cathode is the place where reduction reactions take place. The most common reduction reactions occur at the cathode are: Pourbaix diagrams are used to map the electrochemical phases of materials in various domains (e.g., metals, oxides, hydroxides, oxyhydroxides, and aqueous ions) over specific ranges of pH and electrode potential. Details related to the applications and uses of Pourbaix diagrams are omitted from this review, and readers are referred to a recent review article 55 .
In case of chemical dealloying, a strong acid is usually used as an electrolyte to form porous structure in an alloy. For example, in case of Au-Ag binary alloy, nitric acid (HNO 3 ) is used to remove Ag from the alloy to form nano-porous structure 56 .
To describe the core reaction of high-temperature metal oxidation, we take the example of zirconium oxidation at elevated temperature in air. The reaction of zirconia formation on zirconium surface in the presence of oxygen at high temperature can be given by Zirconia formation is a multistep process: at first, it is a common believe that O 2 molecules are attracted to the metal surface and dissociate into atomic oxygen, which induces Zr to oxidize to cations Zr 4+ and at the same time atomic oxygen become anions O 2− . Finally, Zr 4+ reacts with O 2− to form ZrO 2 . Zirconia is a n-type of semiconductor with a certain amount of anion vacancies in the oxide. Therefore, the oxide growth is controlled mainly by oxygen diffusion through the oxide layer to form new oxide at the metal-oxide interface.
HAC is a form of fracture that initiates and propagates due to the presence of atomic hydrogen inside the metal. The atomic hydrogen is often a result of hydrogen ingress from metal surface during corrosion (the hydrogen molecule in reaction (9) can further dissociate into atomic hydrogen), coating and plating operations, usually referred to as internal embrittlement. More details on its mechanism are given in sub-section "Hydrogenassisted cracking".
High energy particle radiations can result in changes of the morphology, microstructure, and thermo-mechanical properties. These changes are caused by atomic displacements during irradiation. More details on the mechanism are given in the section "Phase-field simulations of morphology evolution in metallic materials under environmental attack". We will also highlight the major differences and similarities related to these environmental attacks in PF modeling in the conclusion section.

PHASE-FIELD METHOD
The composition and microstructure of a metallic material are inherently inhomogeneous. In the PF method, two sets of variables (conserved and non-conserved) are used to describe phase evolution as a function of space and time. A set of conserved variables, c i ðr; tÞði ¼ 1; 2; :::; nÞ, governs the field of concentration variables; here, r and t represent the space coordinates and time, respectively. The non-conserved variables (i.e., phase-field or order parameters), η i ðr; tÞði ¼ 1; 2; :::; mÞ, are used to describe m different phases that may coexist in the system. All of the field (concentration and order) parameters in a PF model have a smooth transition across the interface and thus yield diffuse (finite width) interfaces. A PF model does not require explicit tracking of the diffuse interfaces. This makes the PF models less vulnerable to the singularity (such as stress and/or electric field) at the interface. This feature gives an advantage to the PF method over its counterparts that use sharp-interface techniques. The total free energy of the inhomogeneous system can be written as a function of all of the field variables, including chemical free energy, gradient energy (often referred as interface energy), elastic energy and electrostatic energy: where V is the volume of the system considered; f chem is the bulk chemical free energy density; f grad is the gradient energy density; f elas is the elastic energy density; and f elec is the electric energy density. The gradient energy density can be expressed as where κ c and κ η are the gradient energy coefficients related to the surface energy at the evolving interfaces. The diffusion of solutes and charged species and the evolution of morphology due to the corrosion and phase transition are governed by a set of PDEs under the appropriate initial values and boundary conditions. The evolution of the concentration fields, c i ðr; tÞði ¼ 1; 2; :::; nÞ, and order parameter fields, η i ðr; tÞði ¼ 1; 2; :::; mÞ, is governed by Cahn-Hilliard equations (which describe the transport of solute and charged species), and Ginzburg-Landau equations (which describe structural evolutions), whereŨ represents the velocity field in the fluid (electrolyte), and M i is the mobility of species i. The phase-field equations are coupled with mechanical equilibrium equations describing local stress distribution, electrostatic equilibrium equations describing electric potential distribution and classical Nernst-Planck equations describing ionic species transport in the electrolyte. Detailed descriptions of the chemical, electrostatic and elastic free energy density are provided in the following sections wherever necessary. These free energy terms are treated differently by different models, and readers are referred to two early review articles on the PF method for further details 46,57 . The total free energy of the system is minimized by solving Eqs. (14) and (15) together with other relevant equations.

PHASE-FIELD SIMULATIONS OF MORPHOLOGY EVOLUTION IN METALLIC MATERIALS UNDER ENVIRONMENTAL ATTACK
In this section, we discuss progress in the area of PF-based simulations of the morphological evolution of metallic materials under environmental attack. The relevant literature is divided into subsections according to the mechanisms considered while formulating the models and the simplifications made during simulations. In the first subsection on simulating pitting corrosion without passivation kinetics, we discuss the progress made in this area using early PF models. These models make two major assumptions: (1) that part of the metal is exposed to the electrolyte environment due to passive film damage, and (2) that no ICP formation (or passivation kinetics) occurs. The literature has largely ignored pit nucleation (i.e., passive film damage mechanism), a complex phenomenon. Pit nucleation and growth follow time and length scales that can differ by several orders of magnitude: nucleation usually occurs over approximately nanometers to micrometers, whereas growth occurs over micrometers to millimeters. Both processes largely depend on the material properties and the environmental. Therefore, it is very difficult (if not impossible) for a mesh-dependent methodology like PF modeling to combine these length scales using available computational resources.
In the second subsection, we discuss the methodology used and progress made in pitting corrosion studies and incorporate the effect of passivation kinetics. In the third subsection, we summarize and review progress in the prediction of corrosion kinetics related to stress (or mechanical loading) in metallic materials. We also highlight key assumptions applied to the simulations described in these publications. In the fourth subsection, we discuss microstructure-dependent corrosion kinetics, such as IGC, crystal-orientation-dependent corrosion and corrosion in reinforced metal matrix composites. In the fifth subsection, we discuss the application of the PF method to the formation of porous structures during liquid or chemical metal dealloying. In the sixth subsection, we discuss hydrogen-assisted cracking (HAC). In the final subsection, we discuss hightemperature metal oxidation in the gas phase.
Furthermore, microstructure evolution under irradiation represents an important class of metal degradation. Complex microstructural changes can occur in structural materials, such as those used in nuclear power reactors. Over the last two decades, the PF method has been used extensively to model degradation in irradiated metallic materials, and our group has published several research articles on this topic [58][59][60][61][62][63] . In this journal, Li et al. 47 published a comprehensive review on the use of the PF method to model microstructural or property evolution in irradiated nuclear materials, and reported that the PF literature on this topic accurately describes important phenomena, such as the migration and recombination of defects, spatial-dependent generation and dissolution induced by irradiation, the Soret effect, strong interfacial energy isotropy and elastic interactions between defects. Li and colleagues further suggested that the PF method can be used to predict 2-D and 3-D microstructural evolution qualitatively and quantitatively while considering complex microstructural effects, such as radiation-induced segregation, secondphase nucleation, void and gas bubble superlattice formation, void migration, hydrate formation and grain growth 47 . In their review, they also discussed the strengths and limitations of PF method for studying irradiation effects in nuclear materials. The readers are referred to the original review article 47 for further information.
Pitting corrosion without passivation Mai et al. 64 proposed a PF model of pitting corrosion in steel, in which the chemical free energy of the system was defined according to the model proposed by Kim et al. (i.e., the KKS model) 65 . Although this 64 work used the Tafel equation to link the kinetic interface parameter with overpotential (and thus ignored cathodic reactions), the kinetic interface parameter was assumed to be constant in the numerical simulations. In this early effort, a simplified PF model for pitting corrosion that ignored some important parameters was proposed. In short, this work ignored the effects of cathodic reactions, the electrolyte, ionic transport, overpotential dependence on ionic concentrations, passivation kinetics and mechanical loading on the pitting corrosion process. In 2018, researchers proposed several PF models for pitting corrosion that incorporated the electromigration effect and has ability to show a smooth transition from T.Q. Ansari et al. reaction-controlled to migration-controlled followed by diffusioncontrolled process. Xiao et al. 66 developed a PF formulation and benchmarked their simulation results with a sharp-interface model. The implementation limitation of their model to 1-D geometry can be considered as the biggest drawback. To overcome this problem, Ansari et al. 67 developed a PF model and demonstrated with their simulations that the predicted pit depth is quite close to the experimental results 68,69 . Ansari and colleagues also performed several case studies which showed the applications of their proposed formulation to capture the evolving surface morphology of complex microstructures. In 2018, Mai and Soghrati 70 , Chadwick et al. 71 , and Tsuyuki et al. 72 proposed PF models that would incorporate the effects of ionic species concentration in the electrolyte and the electric field strength and could also handle complex microstructures. PF models differ according to the definition of free energy density used in the system. Here, we briefly summarize the model proposed by Ansari et al. 67 , which applies the definition of chemical free energy used in the KKS model. Importantly, several other methods can be used to describe the chemical free energy density of a system. For example, the formulation proposed by Wheeler-Boettinger-McFadden (WBM model) 73,74 has been used widely to model phase change or transformation processes. For this model, the major drawback is that it places more limitations on interface thickness than the KKS model. The mesh element size used for numerical scheme depends on the interface thickness; therefore, the ability to use a relatively coarse mesh in the KKS model can significantly reduce the computational time. Similarly, formulations proposed by Khachaturyan 75 and Folch-Plapp 76 have been used extensively in the PF modeling literature. More details about these formulations can be found in the original articles.
The model 67 introduces an order parameter to describe two phases (solid and electrolyte) in the system. Evolution of this order parameter is governed by the Ginzburg-Landau Eq. (15), whereas evolution of the concentration field is governed by the Cahn-Hilliard Eq. (14). The total free energy of the system is expressed in terms of the chemical, gradient and electric free energy densities. The chemical free energy density, f chem (c 1 ,η), of the system is adopted from the model proposed by Kim et al. 65 for binary alloys. In the KKS model, each material point is considered as a mixture of two co-existing phases that are assumed to have local equilibrium: where c 1 represents the normalized concentration in the total system wherein c s and c l are the normalized concentrations in the solid and liquid phases, respectively, and h(η) is a monotonous function of the form η 2 ðÀ2η þ 3Þ, which takes the value 0 at h(0) and 1 at h (1). Therefore, f chem (c 1 ,η) can be given by the expression where w is the height of the double-well potential, and g(η) is the double-well function given by η 2 1 À η ð Þ 2 . The free energy densities, f s c s ð Þ and f l c l ð Þ, are assumed to be parabolic functions, consistent with the literature 77 . Transport in the electrolyte is governed by the classical Nernst-Planck equation: where i represents the ionic species present in the electrolyte; D i eff is the effective diffusivity of the ionic species in the system; and R i is the rate of the production or consumption of the species as a result of electrochemical reactions happening in the system. The conservation of charge in the system is governed by the Laplace equation. The kinetic interface parameter, L j , is related to the Buttler-Volmer equation. The PF model is solved using the finite element method, in which the 1-D and 2-D results are validated against 1D pencil electrode and 2D foil experiments, respectively, as shown in Fig. 1a, b. The corrosion process appears to be diffusion-controlled because the relation between the square root of time and the pit depth is approximately linear. The authors 67 further describe the evolution of all anions and cations available in the electrolyte and the electric potential distribution. Furthermore, Xiao et al. 66 and Mai and Soghrati 70 compared the distributions of anions and cations in an electrolyte using a sharp-interface model that was proposed by Sharland a few decades earlier 5,6 . Figure 1c shows the pitting corrosion kinetics for two individually growing pits. Few experimental 78 and numerical 79 studies have been devoted to study the interactions of multiple pits growing in a metallic material. Such interaction may result in an increase or decrease in the total corrosion loss, depending on the material and corrosive environment. Ansari et al. 67 reported that an interaction between two closely located pits leads to an increase in corrosion loss relative to that observed when the two pits never interact. A future PF study that investigates this phenomenon in detail by varying the distance between the initial pits, the electrolyte solution and the applied potential conditions may yield interesting results. Figure 1d, e present the pitting corrosion kinetics results reported by Mai and Soghrati 70 and Chadwick et al. 71 . In these simulations, the bulk electrolyte and reference electrode were placed directly at the metal surface. Such an approximation is valid when the process is controlled by an applied electric field but may not hold for an aggressive electrolyte. Through a 1-D analysis (verified by a sharpinterface model), Xiao et al. revealed huge differences in pH, metal ion molarity, chloride ion molarity and electric potential between the crevice bottom and the bulk electrolyte 66 . Tsuyuki et al. 72 were the first to report the corrosion morphology of a 3-D case determined using the PF method, as shown in Fig. 1f; however, their model lacked numerical and experimental validation (for all 1-, 2-and 3-D cases), as identified by the authors. Importantly, in all of the above PF models except that proposed by Xiao et al. 66 , the formulation can be extended easily to study a 3-D case. However, the formation/ dissolution of corrosion products was not explicitly included in studies of passivation kinetics based on these models. In the following section, we will discuss progress in the modeling of pitting corrosion kinetics in terms of incorporating the growth of soluble corrosion products and ICPs and their effects.

Pitting corrosion with passivation
The oxidation of metallic metals in an electrolytic environment often yields corrosion products that can be soluble or insoluble. However, the numerical modeling literature has largely ignored the effect of this phenomenon, and only a handful of theoretical modeling studies have studied the effects of ICP nucleation and growth during pitting corrosion. For instance, studies on pitting corrosion that have used ALE, a sharp-moving-interface method, have demonstrated that the corrosion rate and pit morphology are highly dependent on the new ICP phase [37][38][39] , which can become the rate-controlling parameter if it partially or fully covers the corroding pit. As mentioned in the "Introduction" section, however, the usefulness of the ALE is limited in the context of complex evolving morphologies.
Only a few notable publications in the area of PF modeling have explicitly discussed the effects of the electrolyte, electric field strength and newly formed porous ICP phase on the corrosion rate and pit morphology [80][81][82] . In the first such publication, Ansari et al. 80 proposed a multi-PF model (MPF) of the evolution of a new ICP, along with other microstructural effects of the corroding metal, for investigating the passivation kinetics of pit stability and growth. In the following year, Chen and Ruan 82 extended this MPF model by incorporating strain energy to include the effect of a constant mechanical load. The non-conserved PF equations in these models are derived from classical rate theory. Bazant summarized this theory of chemical kinetics and charge transfer, which is based on non-equilibrium thermodynamics, in his lecture notes 83 and a review article 84 . Since then, his research group has used the PF formulation based on this approach to study the intercalation kinetics of nanoparticles in Li-ion batteries [85][86][87] , while a team lead by Chen has used this approach to study dendrite growth in battery materials [88][89][90] . The PF equations for pHdependent pitting corrosion modeling, which were reported by Tsuyuki et al. 72 , are also based on this formulation but do not form an MPF model. Here, we briefly summarize the MPF model proposed by Ansari et al. 80 . We will address the formulation related to mechanical loading as proposed by Chen and Ruan 82 in the next section on SCC.
The net rate (R net ) of any chemical reaction can be expressed by the difference between the forward and backward reactions. The net rate, which describes the reaction kinetics of a process and satisfies the detailed balance of chemical potentials in a system, is expressed as where the excess chemical potential at the transition state is given by μ ex TS , and the chemical potentials in the initial and final state are given by μ 1 and μ 2 , respectively. The equilibrium constants for the forward and backward reactions are given by r ! 0 and r 0 , respectively, and can be assumed to be equal for a properly defined chemical potential (as discussed by Bazant 84 ). According to Chen et al. 90 , the PF parameters in this model are considered to evolve according to Eq. (18) as a function of thermodynamic force (μ k ¼ μ 2 À μ 1 ), and it is assumed that r ! 0 ¼ r 0 . Therefore, the PF equations can be re-written as where α is the charge transfer coefficient. μ k is expressed as a variational derivative of the total free energy of the system. The readers are referred to the original article 80 for the detailed derivation. Here, we give the final form of the PF equations for metal (η m ), the ICPs (η icp ) and the electrolyte (η e ): where L σ , L ψ1 and L ψ3 are the kinetic interface constants; a i represents the activity of the ionic species I, which is given by a i ¼ χ i c i ; and χ i is the activity coefficient. A saturation factor, S c ¼ 1 À c Fe þ2 =c sat;Fe þ2 , is also introduced in Eq. (22) to enable the inclusion of the effect of the metal ion concentration at the corroding surface in the PF equation for metal dissolution. The second terms on the right sides of Eqs. (22) and (23) are multiplied by the variables λ M and λ ICP , respectively, which are introduced to localize the electrochemical reactions at the evolving interfaces and expressed as λ M = 4(H m H e + H m H icp ) and λ ICP = 4 (H e H icp ). In a classical PF model with no more than two phases, a polynomial function with a non-zero value, 6η À 6η 2 , is usually applied at the interface. Moelans 91 proposed a formulation to express the phase fractions (H i ) for an MPF system, which was adopted in their model as where N is the total number of phases. The time-dependent transport of ionic species in the system is governed by T.Q. Ansari et al.
Nernst-Planck Eq. (19). The electric field distribution in the electrolyte and porous ICP phases is governed by a simple Poisson equation, ∇ σ eff ∇φ Â Ã ¼ I R , where I R is the current density, which is related to the reaction rate (please see the original article for more details 80 ). The effective diffusion coefficient of the ionic species inside the porous ICP phase is assumed to follow the Bruggeman relation 92,93 .
The schematics of the problem studied in their work 80 is shown in Fig. 2a. The corrosion rate was validated using an earlier experimental study conducted on mild steel in Swansea Public Wharf (SPW), as shown in Fig. 2b. Ansari et al. 80 simulated a 1-D geometry that included the growth of ICPs and varied the applied potential to demonstrate that their model could capture reaction-, migration-controlled and diffusion-controlled pitting corrosion kinetics. Using a 2-D geometry, they 80 performed a similar investigation to demonstrate the evolution of metal corrosion and ICP growth and the corresponding variations in the pH, metal ion concentration and electric potential of the electrolyte. In this scenario, metal corrosion occurs during the porous ICP phase when the pit is completely covered with ICP (see Fig. 2c). The maximum metal ion concentration and electric potential occur at the ICP/metal surface (pit bottom), whereas the corresponding values at the ICP/electrolyte surface are significantly lower. This gradient suggests that the corrosion rate is limited by slow diffusion through the ICP phase. Ansari et al. 80 also used the model to simulate under-deposit corrosion, and the results indicated a significantly lower pit growth rate when the pit surface was fully covered with the ICP phase. The results of this under-deposit study are shown in Fig. 3a. For additional details about this comparison, the readers are referred to the original article 80 .
Importantly, they assumed a constant porosity value for the ICP phase in most case studies and later discussed the effect of variation in this porosity value on pitting corrosion kinetics. Figure   3b 1 -b 4 shows the results of pit corrosion at one applied potential value under the effect of a simulated change in porosity. The results suggest that the corrosion process can be reactioncontrolled or diffusion-controlled at the same applied potential, depending on the porosity of the new ICP phase. Notably, the authors skipped the theory discussions related to the reaction kinetics that result in ICP with different porosity. For more detail, some early reports on the experimental investigation of the ICP might provide a good background on this process [94][95][96] . Chen and Ruan investigated the effects of variations in electrolyte (i.e., more and less aggressive cases) and constant mechanical loading on the corrosion rate. As shown in Fig. 3c, they found that a more aggressive electrolyte led to the formation of a deeper pit, consistent with the findings of earlier experimental studies. Their simulation results show that the inclusion of a mechanical load may cause the transition of a smooth pitting corrosion process to SCC. Furthermore, their PF model ignored the effect of mechanical loading on the ICP phase, although this variable may significantly affect passivation kinetics by damaging the new ICP. In the following section, the mechanism related to SCC and the relevant PF studies are discussed in more detail.
Stress-corrosion cracking SCC is an aggressive form of localized corrosion that accelerates the damage caused by existing cracks in a corrosive environment. In an early effort, Stahle and Hansen 97 developed a PF model based on the elastic free energy density that could capture the growth of an anisotropic crack from the metal surface. However, this phenomenological model lacked reaction kinetics. Mai and Soghrati 98 extended their earlier work on pitting corrosion 64 to study SCC initiated from pits in steel and obtained a total free energy definition of the system that was essentially the same as that reported by Mai et al. 64 . In this model, a spatially varying kinetic interface parameter incorporates the effective stress distribution due to mechanical loading via a linear relationship and enables inclusion of the effect of overpotential. Nguyen et al. 99 improved the model developed by Mai et al. 64 to study SCC in steel by including the elastic strain free energy term in the total free energy. Subsequently, Lin et al. 100 developed a more comprehensive model that incorporated both strain and the electrostatic free energy density. Recently, Cui et al. 101 proposed a PF model for SCC driven by dissolution, following the formulations described in two earlier works on PF modeling of SCC 64,99 . In the model by Cui and colleagues, the kinetic interface parameter is considered as a time-dependent spatially varying variable that is related to the film rupture-dissolution-repassivation (FRDR) mechanism. Although this work represents an advance relative to its counterparts 64,99 , it ignores the effect of the electrolyte (or electromigration process) on evolution of the corroding surface. In terms of the PF formulation, the model proposed by Lin et al. 100 is more comprehensive because it incorporates the effect of the electrolyte while retaining the other free energy terms (chemical, gradient, and strain). Nguyen et al. 99 , Mai and Soghrati 98 and Cui et al. 101 described the chemical free energy densities in their proposed models with reference to the KKS model 65 , whereas Lin et al. 100 followed the strategy proposed by Bazant 84 to derive PF equations based on rate theory. For the definition of chemical free energy density based on the KKS method, please see the section "Pitting corrosion without passivation" and Eqs. (14)- (17). For the derivation of the PF equations based on Bazant's work, please see the section "Pitting corrosion with passivation" and Eqs. (20)- (24). Next, we briefly discuss the strain energy formulation, proposed by Cui et al. 101 , which contains both elastic and plastic energy density terms, expressed as where ε e is the elastic strain, ε0 p is the plastic strain rate, C 0 is the stiffness matrix, and σ 0 is the Cauchy stress tensor for the undissolved solid. The isotropic power-law hardening behavior is assumed by defining the flow stress as a function of plastic strain: where E is Young's modulus, σ y is the yield strength, and N is the hardening exponent. In the chemo-mechanical term proposed by Gutman, the accelerating effects of mechanical stresses and strain on corrosion are related to the corrosion current density 102 . This corrosion current density is further related to the kinetic interface parameter (L), as proposed in an earlier work by Mai et al. 64 . This work 101 updates the definition of corrosion current density by including the effects of plastic and yield strain in the acceleration of corrosion. Nguyen et al. 99 investigated the transition from pitting corrosion to SCC in two cases. In the first case, the steel sample is polarized at an applied potential of −570 mV/SCE, with a tensile stress of 110-132% of the yield stress (σ yield = 290 MPa). In the second case, the same sample is polarized at an applied potential of −640 mV, with a constant loading of 125% σ yield . Nguyen and colleagues used their PF formulation to simulate these two scenarios, as shown in Fig.  4a 1 , a 2 . In the first case, the pit grew wider due to the higher polarization potential and the applied mechanical loading, respectively. In the second case, the corrosion process was dominated by SCC, resulting in a narrower, deeper pit. The authors further observed good agreement between their simulation findings and the experimental results. Under conditions of mechanical loading, the morphology of a pit also depends on the nature of the electrolyte, as demonstrated by these two cases. Therefore, inclusion of the electrostatic free energy density might be important in some studies of SCC, as identified by Lin et al. 100 .
In a case study, Mai and Soghrati 98 investigated the effect of material properties on SCC, as shown in Fig. 4b 1 , b 2 . As they increased the strain hardening (E p ) value of the material from 2 to 10 GPa, they observed a reduction in the stress intensity factor (K) that resulted in a lower crack depth at a higher E p . In a later study, Cui et al. 101 investigated the interactions of multiple pits undergoing SCC, as shown in Fig. 4c. In the illustration of three pits (defects), all but the middle pit is in contact with the electrolyte; material failure occurs as the pit cracks grow towards the defect in the middle. These reported cases demonstrate the ability of the PF method to predict SCC in metallic materials. Although these models make assumptions that might not hold in all practical cases, the proposed methodology can be extended to investigations of SCC in industrial settings.

Microstructure-dependent corrosion
Pit nucleation greatly depends on microstructural features at the exposed metal surface, suggesting that material features of the bulk metal also play an important role in controlling the corrosion rate and morphology. Consequently, most corrosion studies have focused on the effects of these features on corrosion acceleration Fig. 4 Stress-corrosion cracking. a 1 and a 2 Evolution of a pit over time; SCC is dominant in the second case 99 . b 1 Evolution of crack growth and corresponding stress distribution at three different strain hardening (E p ) values. b 2 Difference in crack depths in three different cases 98 . c Study of the effect of SCC on the growth of multiple pits 101 . a 1 and a 2 are reproduced with permission from ref. 99 , published by Elsevier in 2018. b 1 and b 2 are reproduced with permission from ref. 98 , published by Elsevier in 2017. c is reproduced with permission from ref. 101 , published by Elsevier in 2021.
T.Q. Ansari et al. or resistance. Experimental studies have indicated that microstructural features greatly affect the corroding surface morphology. In an early effort, Suter and Alkire investigated the effect of inclusions on pit morphology 103 and reported that deep inclusions resulted in an irregular corroding surface with stable pitting, whereas shallow inclusions were responsible for unstable pit growth. Several experimental studies have determined that the corrosion rate depends on the crystal orientation. Single crystal studies have revealed that the corrosion rate is maximized along the [100] orientation and minimized along the [111] orientation in aluminum 27,104 . Similarly, corrosion was shown to occur preferentially along the sensitized grain boundaries of heat-treated stainless steel and aluminum alloys 105,106 . Here, we briefly discuss progress in the use of the PF method to study these effects.
IGC in sensitized metals: IGC has received considerable research interest over the last few decades both experimentally and theoretically. In a series of publications, Frankel and team reported the IGC kinetics of aluminum alloys and laid the foundation for the future design of alloys with better IGC resistance [107][108][109][110] . Several relevant modeling studies based on the statistical 110 , FEM 111 , CA 33,36 and PD methods 41 have been published. Recently, Ansari et al . 112 proposed a quantitative MPF model for determining the IGC kinetics of sensitized materials. This work showed good agreement with the experimental works of Lim et al. 105 and Bonzom and Oltra 106 . In their work 112 , they considered grains and SGB as materials with different material properties and therefore introduced an MPF formulation that considered grains, SGB and the electrolyte as three distinct phases. In another study, Nguyen et al. 113 proposed a PF model for inter-granular and trans-granular SCC that only used the effect of SGB (i.e., did not consider it as a different phase) and ignored the electrolyte influence.
Ansari et al. 112 compared the predictions made using their MPF model with experimental results available in the literature, as shown in Fig. 5a. In the images, the corrosion depth predicted by the model seems to be in good agreement with the experimental results, although the widths of the corroded SGBs appear larger than those in the experimental images. The authors assumed that the SGB was a few orders of magnitude larger than the real SGB to reduce the computation cost. Ansari et al. 112 also reported a 3-D case study of IGC kinetics, as shown in Fig. 5b; this was also limited to a smaller-sized geometry and thicker SGB. The computation cost associated with such studies can be overcome by using an adaptive meshing technique 114,115 . Most recently, Gao et al. reported a robust strategy that can be used to solve PF models and thus reduce the computation cost 116 .
Crystal orientation-dependent corrosion: Early experimental studies revealed that different exposed crystal orientations exhibit differences in corrosion potential and pitting behavior 104  Some PF models for pitting corrosion can also be used to study the effect of crystal orientation on corrosion surface evolution 64,67,70,71,80 . Mai et al. 64 and Ansari et al. 67 used a PF approach in which all crystal orientations were described by a single-order parameter. Each crystal orientation was given a different corrosion potential, resulting in different corrosion rates along different orientations. Chadwick et al. 71 and Ansari et al. 80 used an MPF approach wherein each crystal orientation is represented by a single-order parameter. Chadwick et al. further discussed the advantages of using multiple versus single-order parameters for all crystal orientations. In all these models, only three principal crystal orientations, [111], [100] and [110], are assumed to face the electrolyte. PF modeling studies are needed to quantitatively verify the simulation results using experimental findings and identify potential key features of corrosion resistance. Crystal orientation engineering has been used to improve the properties of other materials. For example, Li et al. 121 recently used crystal orientation engineering to increase the breakdown strength of ceramic capacitors, which in turn increased the energy storage density. In this field, quantitative PF modeling research might promote the development of corrosion-resistant metallic materials.
Mai et al. 64 studied the evolution of corrosion morphology as a function of three principal crystal orientations, as shown in Fig. 6a.
In the resulting PF model, a single-order parameter is used to describe the metal phase along with all crystal orientations. The model incorporates the effect of the crystal orientation-dependent corrosion rate by considering different values of the kinetic interface parameters along the crystal orientations. Ansari et al. 80 used an MPF approach wherein each orientation was governed by a single-order parameter, as shown in Fig. 6b. Both strategies can be used to predict non-uniform crystal orientation-dependent pitting corrosion. All of these studies suggest that the PF approach is a useful tool for the quantitative prediction of pitting corrosion while considering the crystal orientation. This strategy can be coupled with experimental setups wherein orientation image microscopy (OIM) is used to identify the orientations in each sample, followed by corrosion experiments involving polarization. Such comparisons can enable more quantitative validation of the ability of the PF method to predict crystal orientation-dependent corrosion.
Corrosion in metal matrix composites: Ceramic particles, such as TiB 2 and TiC, are usually embedded in steel to increase the wear resistance, strength and stiffness. Although these particles have very little effect on pitting corrosion in the absence of mechanical loading, the effect becomes non-negligible when the material undergoes SCC. In publications on the PF approach, Mai et al. 64 and Ansari et al. 67 used case studies to demonstrate that pitting corrosion had very little effect on the corrosion rate in ceramic particle-reinforced steel. However, these authors discussed the possibility of accelerated corrosion at the ceramic particle-metal interface due to increased local stress under the condition of external loading. In a later study on SCC, Mai and Soghrati 98 investigated this effect in the presence of external loading. Figure  7a, b show higher corrosion rates of the metal along the interface with the ceramic particles due to increased local stress at the interface. Recently, Lin et al. 100 investigated the effect of variable yield strengths of ceramic particles and fibers on the corrosion rate and found that metals embedded with softer particles (i.e., lower yield strength) had higher corrosion rates, whereas those embedded with stiffer particles showed more corrosion resistance under constant loading, as shown in Fig. 7c.
Porous structure formation during dealloying Typically, dealloying is an undesired corrosion wherein the less noble component of the alloy is dissolved, leading to crack growth, SCC or corrosion fatigue. However, dealloying can also be beneficial to create porous structures with significantly different properties. Previous studies have applied this technique to binary alloys, such as Ag-Au, Cu-Au, Ni-Pt and Mn-Ni, to produce porous structures with many industrial applications 122,123 . Similarly, liquid metal (LM) dealloying has been used to remove one component from a binary alloy by using a liquid bath of another metal element 124 ; here, the element with a lower melting point may be dissolved from the alloy. This technique can also result in the formation of a nano-porous structure. Pattern formation during dealloying remains poorly understood because most progress in this area has been made by experimental determination of the key parameters that control the process. A theoretical modeling framework based on the PF method might elucidate key factors in the formation of nano-porous material structures and enable better control of the process. Geslin et al. 125   interplay between interfacial spinodal decomposition and the diffusion-coupled growth of enriched solid and liquid phases in the alloy. As the PF simulations were performed at a much greater length scale due to computational limitations, the authors derived the scaling laws governing dealloying kinetics based on their experimental results. Recently, Li et al. 56 developed an MPF model to study the chemical dealloying process and nano-porous structure kinetics of an Au-Ag binary alloy in an acidic electrolyte. In this model, the evolution of three phases, namely, solid Au-Ag, Au-rich and electrolyte phases, is governed by three nonconserved order parameters (ϕ i ) that are expressed by a summation term: where L is the kinetic interface parameter, and F is the Gibbs free   energy of the system. The first term on the right-hand side of Eq. (30) is the gradient energy; the second term establishes the double-well potential between the two binary phases; and the third term expresses the chemical energy free energy variation with respect to the non-conserved parameter, as given by the detailed expression in Eq. (31). c j Ag and c j Au represent the concentrations of silver and gold element in phase j, respectively; h j is the interpolation function; and f chem is the chemical free energy in phase j. The last term in Eq. (30) is the noise term used to represent the nucleation of a new solid phase-the gold-rich phase at the alloy-electrolyte interface. The evolution of Ag and Au are governed by a Cahn-Hilliard equation as expressed in Eq. (14), with slight modifications. After verifying their approach against experimental work 126 , Li et al. 56 investigated the influence of the acid concentration on nano-porous structure formation and the corrosion rate. Li et al. 56 found that a more concentrated acid solution results in a much faster dissolution of Ag (as expected) but generates a different nano-porous structure. As shown in Fig.  8f, the Au-rich ligaments are much smaller and the spacing between the ligaments is smaller in a more concentrated acid solution. The model can also predict a threshold initial concentration of Au in the alloy, above which the porous structure cannot form. This prediction agreed well with experimental observations. This case study paved the way toward future attempts to control and manipulate the size, shape and morphology of nano-pores. Li et al. 56 also studied the effect of inert inclusions (usually ceramic particles embedded in alloys to enhance mechanical properties) on the evolution of the nano-porous structure, as shown in Fig. 8g, and are currently investigating whether the inclusion of an external electric field might significantly affect nano-porous structure formation. The findings from these studies will improve understanding of the chemical dealloying process and enable better control and manipulation of the pore formation process.
Hydrogen-assisted cracking Often as a corrosion product, hydrogen can diffuse through the crystal lattice of a metallic material and induce fractures under mechanical loading that grow, resulting in catastrophic failure. Several modeling techniques have been used to study HAC, including dislocation dynamics method 127 , weakest-link approaches 128 and cohesive-zone models 129,130 , that can incorporate reductions in material strength due to hydrogen. Although these theoretical works can model many aspects of HAC, their ability to address important challenges, such as interacting cracks, mixed modes and randomly distributed initial cracks, is limited. On the other hand, the PF method has been used successfully to model brittle [131][132][133] , ductile 134,135 , cohesive 136,137 , dynamic 138,139 and multi-physics fractures 140,141 .
PF models for fracture problems often depend on energy minimization to determine when and where the crack nucleates and how much and along which direction it will propagate 142 . The PF method introduces scalar diffused interface parameters to represent the crack phase and create a diffused band of damage, rather than a sharp band where a stress singularity can yield unrealistic computational results. This aspect of the PF method is discussed in the section "Stress-corrosion cracking". PF models proposed for fracture are well suited for HAC because they can incorporate the fracture energy or fracture stress/strain as a function of the hydrogen concentration. An early PF model 143 studied the evolution of hydrogen concentration around a crack and predicted HAC in ferroelectric ceramics using Ginzburg-Landau strain field equation. Since then, several attempts have been made to model HAC in metallic materials using the PF method [144][145][146][147][148] .
Panẽda et al. 146 proposed a PF model for HAC wherein the fracture process and crack growth are based on a coupled mechanical-diffusion formulation. This model uses a classical PF formulation for fracture, coupled with Fick's law for mass diffusion and DFT calculations based on the relation between the hydrogen concentration and fracture surface energy. Panẽda and colleagues used FEM to numerically solve different HAC situations, such as a cracked plate, notched cylinder and internal cracking in a specimen. They also discussed and presented case studies to elaborate the applications of their model to modeling the degradation of metals in aqueous environments.
As shown in Fig. 9a, a cylindrical bar with a notch in the center is placed under tensile stress. Panẽda and colleagues selected a finer mesh proximal to the notch and a coarser mesh distal from the notch and simulated the effect of a variable hydrogen concentration (in the ppm range) on this cylindrical bar. (Fig. 9(b), plot of the simulated net section strength of the bar vs. hydrogen concentration). The simulation results were consistent with the experimental results reported by Wang et al. 142 . Their work 146 is also limited to brittle fracture, similar to earlier studies on PF modeling of HAC. Recently, Wu et al. 147 proposed a cohesive zone model (CZM) that incorporates the cohesive fracture and evaluates the effect of hydrogen in terms of reducing the fracture energy or cohesive zone strength. To study HAC, Wu and colleagues extended their earlier CZM-PF model, which was intended to address purely mechanical problems, used FEM to numerically solve their proposed model and presented several representative studies of crack growth due to hydrogen. For this work, they chose a 2-D sheet containing a notch at the top surface and a crack in the middle. In Fig. 9c, d, a comparison between two sheets with different interface thicknesses (0.4 and 0.2 mm) in the absence of hydrogen reveals that the crack morphology is independent of the interface thickness parameter in the PF model. Similarly, Fig. 9e, f shows that even when hydrogen is introduced, the model results remain independent of the interface thickness parameter. Wu et al. introduced two models for failure: model-1, which includes a constant ratio of failure strength and fracture energy; and model-2, in which the Irwin length remains constant 149 . These two models were simulated using the PF-CZM method, and the results are shown in Fig. 9g. The PF-CZM method predicted the displacement to be highest in the absence of hydrogen. Compared with model-1, model-2 yielded a smoother transition between applied force and displacement. The authors further determined that an increase in the hydrogen concentration resulted in a decrease in displacement.
The studies described above demonstrate the great potential of the PF method for solving complex evolving geometries and performing validation studies. Additional quantitative PF studies are needed to accurately predict failures, provide engineering guidelines and recommendations to avoid accidents due to HAC and design better hydrogen-resistant materials.
High-temperature metal oxidation A protective oxide layer may or may not form on the surfaces of metals in gas phase at high temperature. Chen and colleagues 150 developed an early PF model to study oxidation kinetics under a dual oxidant atmosphere at high temperature. Specifically, they used Fick's law to model the formation of metal oxide or metal sulfide based on a diffusion process. Their model was limited to 1-D geometry, and thus excluded several important factors, such as the effect of Pilling-Bedworth ratio (the volume ratio between the metal oxide and the consumed metal), defects, and surface morphology. Therefore, this work yielded only qualitative results. Abubakar et al. extended this model to study hot corrosion kinetics in a V 2 O 5 thermal barrier coating 151 , specifically the mechanical behavior of the coating due to phase transformation. Their results suggested that failure of the top-coating layer could be attributed to grain boundary sliding. The research team led by D.N. Fang has developed PF models to study oxidation of alloys 152 , ZrB 2 /SiC composite ceramics 153 , and oxidation at a crack tip 154 . In  a Oxide scale (h ox ) and b interface roughness (Ra) versus time. c Schematic of interference between neighboring oxide bumps as they grow. d Schematic of metal oxide growth in a gaseous environment. a-d are reproduced with permission from ref. 162 , published by nature springer in 2020. these publications, the effect of PB ratio was considered and used to estimate the stresses in the oxide layer, the substrate or at the crack tip. These studies resulted in good agreement with experimental observations and provided useful insight for failure mechanisms.
Extensive investigations have been conducted on the oxidation kinetics of zirconium. Zirconium alloys are used in nuclear power industry as a cladding material at elevated temperatures. Therefore, the zirconium surface must be covered with a protective layer of zirconia (oxide of zirconium). During oxidation, however, the volumetric eigen strain resulting from the change in specific volume from zirconium to zirconia (reason for PB ratio) can lead to crack nucleation and growth. Experimental techniques, such as x-ray diffraction and Raman spectroscopy, can be used to identify these growth stresses, which can be as high as 5 GPa 155 . The growth stress magnitude depends largely on the surface morphology. The PF method has been used to model oxidation kinetics in zirconium alloys. For example, Ammar et al. 156 proposed a PF model wherein the metal-oxide interface is driven by oxygen diffusion. Cheng et al. revisited Wagner theory with their diffuse-interface (PF) model. They reported that Wagnerparabolic law is violated when either the oxide film is of the order of Debye length or interfacial reaction is rate-limiting 157 . Zaeem et al. 158 incorporated the elastic strain energy in their PF model to describe the metal-oxide phase transformation and address the stress-oxidation interaction. Zhao et al. 159 further improved this model by incorporating the inelastic deformation to address hightemperature oxidation. Sherman and Voorhees 160 investigated the electrostatic effect on initial oxidation by solving Poisson equation in their PF formulation. Sait et al. proposed a PF model with a modified Ginzburg-Landau equation coupled with mechanical analysis to study oxide growth in thermal barrier coating. They used rate-dependent Gurson-type plasticity theory to model the inelastic behavior in porous ceramic coating 161 . Recently, Lin et al. 162 extended their earlier model 163 , which addressed outward inhomogeneous oxide growth and surface roughening of NiO, and developed a more comprehensive PF model to study inward and outward oxide growth and surface roughening/smoothening processes. This model not only considered the PB ratio through a volumetric eigen strain, but also considered the eigen strains caused by diffusing species in the oxide such as the anion O 2− and cation Zr 4+ . The model also considered the effect of creep on stress relaxation. As shown in Fig. 10d, a schematic of the oxidation kinetics studied in their work shows that the inhomogeneous oxide can grow outward or inward on the metal surface. Under higher localized stress, metal surface roughening may occur due to inward oxide growth, as shown in Fig. 10a-c. Specifically, Fig. 10a shows that the oxide thickness (h ox ) decreases as the volumetric eigen strain (PB ratio) increases, as verified by the experimental results. Similarly, an increase in volumetric eigen strain (PB ratio) results in increased surface roughness (Ra), as shown in Fig. 10b. The model can also predict the locations of maximum stresses within the oxide, which can often result in oxide cracking. It should be noted that most of these PF methodologies are not limited to zirconium and can be used to study oxidation kinetics in other metals in gaseous environment.

CONCLUDING REMARKS AND FUTURE PERSPECTIVES
This review focuses on progress made in PF modeling of metallic materials evolution under various environmental conditions. In general, PF and MPF models have been developed, which can account for new phase nucleation (such as 125 ), chemical reactions on metal surface and/or in the environment, bulk and surface diffusion, ion migration, electrostatic, and mechanical loading.
The results of most PF models have been verified against earlier numerical modeling results or experimental results obtained using simpler geometries. To date, the validation of PF models for complex geometries has been limited, although some studies have validated the models against established experimental results for complex evolving geometries 112 . However, current PF models for studying corrosion in metallic materials ignore the processes that lead to the breakdown (or pit nucleation) of the passive film, which has long been associated with attack by chemical species (e.g., chloride ions), microstructural features at particular surface sites and/or developed stresses. In principle, these phenomena can be incorporated in PF and MPF models if the real causes of passive film damage are known.
Most of PF models reported in the literature employed both Ginzburg-Landau and Cahn-Hilliard equations, while some used only one of these equations. The reviewed PF models are similar in terms of core methodology (PDEs that govern the morphology and/or composition evolution). However, the major difference lies in the description of bulk chemical free energy, and whether or not to include long-range interaction energy terms such as strain energy and/or electrostatic energy. The definition and description of these free energy terms depend on the degradation mechanism for each case. The phase field equations (Eqs. (14) and (15)) are solved in a multiphysics environment where they are coupled with additional equations that describe (or govern) different phenomena such as mechanical interaction (mechanical equilibrium equations), ion transport (Nernst-Plank equations), electric interaction, and ionic distribution (Poisson or Laplace equation).
Although PF models have been used successfully to predict the evolution of complex degradation processes, research has not yet advanced to the point of guiding the design of degradation-resistant materials. For example, the PF method has not been applied to studies of corrosion fatigue, an aggressive form of material degradation. A quantitative PF model that can predict the service life of a material under cyclic loading and in a corrosive environment may be very useful for solving some industrial problems. Recently, Frankel and colleagues revisited and summarized the corrosion theory to improve understanding of localized corrosion processes [164][165][166][167][168] . This series of publications might be a good reference for future PF studies. Additionally, the computational cost due to the grid size dependence of the interface thickness can be a significant barrier to the application of PF models to larger 3-D geometries. Recently, physics-informed neural networking 169 , a deep learning-based technique for solving PDEs, was introduced. This technique has shown great potential for reducing the computational costs incurred when numerically solving PF models 170 . Possibly, machine/deep learning could be used to train corrosion prediction models based on large datasets generated by high-throughput PF modeling simulations.
The recent development of high-entropy alloys (HEAs) with superior material properties has led to a few modeling studies of the corrosion process in these materials. The successful use of the PF method to model phase transformations in binary and tertiary alloys suggests that a thermodynamically consistent PF model could be used to quantitatively predict localized corrosion kinetics in HEAs and enhance understanding of the process. Similarly, considerable attention has been directed towards additively manufactured (AM) metals over the last decade, as these materials have a diverse range of practical applications. Recent reviews have highlighted the differences in the material properties of these AM alloys from those of conventional alloys, which are due to the internal stresses generated by rapid heating and cooling processes 171,172 . These reviews have also demonstrated the need for a theoretical modeling framework, such as the PF method, that can elucidate the corrosion kinetics in AM alloys.