Kinetic Model of Incipient Hydride Formation in Zr Clad under Dynamic Oxide Growth Conditions

The formation of elongated zirconium hydride platelets during corrosion of nuclear fuel clad is linked to its premature failure due to embrittlement and delayed hydride cracking. Despite their importance, however, most existing models of hydride nucleation and growth in Zr alloys are phenomenological and lack sufficient physical detail to become predictive under the variety of conditions found in nuclear reactors during operation. Moreover, most models ignore the dynamic nature of clad oxidation, which requires that hydrogen transport and precipitation be considered in a scenario where the oxide layer is continuously growing at the expense of the metal substrate. In this paper, we perform simulations of hydride formation in Zr clads with a moving oxide/metal boundary using a stochastic kinetic diffusion/reaction model parameterized with state-of-the-art defect and solute energetics. Our model uses the solutions of the hydrogen diffusion problem across an increasingly-coarse oxide layer to define boundary conditions for the kinetic simulations of hydrogen penetration, precipitation, and dissolution in the metal clad. Our method captures the spatial dependence of the problem by discretizing all spatial derivatives using a stochastic finite difference scheme. Our results include hydride number densities and size distributions along the radial coordinate of the clad for the first 1.6 h of evolution, providing a quantitative picture of hydride incipient nucleation and growth under clad service conditions.


Introduction
Corrosion of metallic structural materials is a pervasive phenomenon in industry and technology [1][2][3][4]. In nuclear reactors, understanding the kinetics of corrosion of metallic components is grand materials science challenge due to the synergistic combination of high temperature, mechanical stresses, complex coolant and fuel chemistry, and irradiation [5][6][7]. In light-water nuclear reactors (LWR) zirconium alloys are used as cladding material in fuel elements to provide mechanical integrity between the coolant (water) and the fuel while keeping low levels of neutron absorption [8][9][10][11]. In principle, Zr clad is subjected to corrosion from the coolant (water) and fuel sides, both by way of oxygen and hydrogen penetration. The oxidation and hydrogenation of zirconium fuel components in LWR may affect reactor safety and efficiency, which makes corrosion a critical design aspect of Zr materials response in nuclear environments [12][13][14][15][16][17][18].
While the majority of the focus of corrosion studies has centered on oxidation and oxygen transport and chemistry in the clad, in the corresponding temperature range, zirconium is known to absorb hydrogen and form hydrides once the critical concentration is reached in the interior of the clad. The accumulation of hydrides during operation plays an important role in fuel performance and safety during steady-state operation and transients, accident conditions, and temporary and permanent fuel storage [19,20]. Examination of the Zr-H phase diagram below 810 K [21][22][23][24] indicates that the first stable compound that appears after the metal solid solution (α-Zr) is a cubic phase with a nominal stoichiometry of 1:1.5 (atomic) known as δ-hydride. Generally, a range of stoichiometries between 1.52 and 1.66 is accepted experimentally as corresponding to this phase [25,26]. While the presence of other metastable hydrides has been reported depending on temperature, aging time, or alloy composition [27], it is now well accepted that the needle-shaped structures that form in the metal region beneath the oxide layer are δ-hydride precipitates. Precipitation first starts when the hydrogen concentration reaches the terminal solubility limit, which ranges from zero at 523 K to approximately 7% at. at 810 K [21]. Although δ Zr 2 H 3 displays good thermo-mechanical stability, it is also an exceedingly brittle phase [28][29][30][31] that can compromise the clad's mechanical integrity [18,[31][32][33].
A key observation of the hydride microstructure is the elongated shape of the precipitates up to a few microns in length [34][35][36][37], typically aligned along directions consistent with the stress distribution within the clad. Calculations and experiments point to the large misfit strains between the cubic δ-hydride and the host α-Zr as the reason behind such preferential alignment [35,[38][39][40][41][42], which may also impact the mechanical response of the clad. Indeed, the formation of brittle hydride phases is a principal cause of delayed hydride cracking (a subcritical crack growth mechanism facilitated by precipitation of hydride platelets at the crack tips in Zr clad [19,[43][44][45]).
The phenomenology of corrosion is such that oxidation and hydrogenation are typically treated separately, despite some evidence suggesting that there might exist synergisms between oxygen and hydrogen pickup and transport that must be considered jointly in corrosion of Zr [46][47][48][49][50]. This is partially due to the formation of a clearly distinguishable outer oxide scale and and inner region where hydride platelets accumulate. In keeping with this distinction, existing models of hydrogen pickup and precipitation have been developed assuming no cooperative effects from oxygen on hydrogen transport and reaction [35,51]. Models based on hydrogen supersaturation of the α-Zr metal [52][53][54] assume a binary partition of hydrogen in the clad, either as solid solution or as part of precipitates without specification of their size, number, or orientation. Detailed cluster dynamics (CD) modeling offers a more accurate alternative to obtain hydride size distributions and number densities by solving the complete set of differential balance equations with one-dimensional spatial resolution [55,56]. Phase field methods can capture extra detail by furnishing the shape and orientation of hydrides in addition to concentrations and sizes [27,57,58].
An important aspect often overlooked in the models when studying hydrogen transport and hydride formation in the clad is that it occurs in a dynamic setting, with the oxide scale growing in time and hydrogen traversing an increasingly thicker layer before it can reach the interface. This is rationalized in terms of sluggish H diffusion through the oxide in the relevant temperature range (<300 • C), suggesting that this then would be the rate limiting step [34,36]. This is the accepted picture during the pre-transition regime, as, after that, fast H transport then occurs through percolated crack networks formed in the oxide layer [59]. However, there is contradicting evidence in the literature about this [60,61], and it is not clear what effect a dynamic boundary condition might have on hydrogen precipitation in the metal substrate at higher temperatures, and what the evolution of the hydride microstructure will be in those conditions. With the objective of shedding new light on these and other issues by using new computational and experimental understanding, in this paper, we present an comprehensive hydrogen transport and precipitation model in Zr formulated from first principles reaction kinetics and fundamental thermodynamics and mechanics. The model is parameterized using electronic structure calculations and experiments and captures both transport across the oxide layer growth and precipitation in the clad under dynamic hydrogen concentration profiles at the oxide/metal interface. First, we describe the fundamental chemistry and phenomenology of the hydrogen evolution in the clad followed by a mathematical formulation of the model. We then provide numerical results under a number of conditions relevant to LWR operation. We finalize with a discussion of the results and the implications of our modeling approach for zircalloy behavior.

Zr-Clad Hydrogen Chemistry
The formation of hydrides in the clad is predicated on exposure of its outer surface to bi-molecular hydrogen. This can occur as a consequence of exposure to water or steam, from the reduction of water molecules as: or directly from exposure to hydrogen gas. It is well known that only a fraction of the hydrogen produced in this way is absorbed by the clad, ranging between 5 and 20% of the total hydrogen uptake (the total amount of hydrogen obtained stoichiometrically from reaction (1) [35,50,[62][63][64]. This, known as the pickup fraction, sets the boundary condition for the adsorption of hydrogen at the clad's surface. Adsorbed H 2 molecules can split into atomic hydrogen by a number of processes [65,66], although whether this atomic H appears in a neutral or charged state in the metal is still an issue under debate [61,65]. Hydrogen atoms diffuse through the oxide layer and reach the oxide/metal interface, from which they can enter the α-Zr substrate and undergo a number of processes depending on temperature and concentration. Above the terminal solubility limit, hydrogen and zirconium react to form a hydride: where x is the atomic hydrogen concentration. By way of illustration, Figure 1 shows representative hydridized microstructures in Zirc-4 and Zirconioum. ed to s Q in scatsmallwhich uthal uthal dified ð2Þ io b/a attern ld & ult of three onescissa ansfer Q*. Since < 1 for elliptically symmetric scattering, the Q* range is shifted to lower values. An example of the mapping of Q* onto the Q x Q y detector plane based on equation (2) is shown in the inset of Fig. 5. A one-dimensional dAE/d versus Q* plot can then be constructed provided is not Q dependent for a given sample. The values from iso-intensity contour fitting listed in Table 2 support this provision. Optical microscopy has been performed to image the deuteride-phase particles at a deuterium concentration of 394 w.p.p.m., within the resolution of the microscope (estimated to be equal to or slightly below 1 mm). A series of optical micrographs are shown in Fig. 6 for the zero-tensileload material (Zy4-33-394d-00MPa, left-hand micrographs) and the isothermal-tensile-load samples (Zy4-32-394d-170MPa, right-hand micrographs). All micrographs show the ND-TD plane with the RD out of the paper normal to these directions. The sample orientation during SANS analysis was an incident beam direction along the ND with the TD-RD plane coinciding with the measured Q plane. The deuteride particle morphologies are similar in appearance and do not support conclusive statements with respect to changes induced by isothermal load. This is consistent with the small amount of deuterium ($40-100 w.p.p.m.) put into solution at 473 K relative to the total concentration (394 w.p.p.m.). The same explanation will hold for temperature cycling under load, although more deuterium will go into solution at the higher    In view of this picture, and to be consistent with our recent work on oxide layer growth modeling [69], we split our model into two connected elements: (i) a transport part involving H diffusion through an evolving oxide layer, and (ii) a kinetic model of hydride formation and growth in the metal with a dynamic boundary condition set by the first part (i). Figure 2 shows a schematic diagram of the geometry considered for this study and the principal chemical processes taking place in the material. Although it is well known that the Zr oxide layer is not monolithic, containing various Zr-O phases depending on the external conditions and alloy composition [69][70][71][72], here we consider a single phase (monoclinic) ZrO 2 with thickness defined by the variable s(t). Hydrogen's diffusion through this layer is thought to occur mostly along grain boundaries, in microstructures ranging from columnar in out-of-pile [71] tests to roughly equiaxed for in-pile conditions [70]. This phenomenon takes place during the pre-transition regime, before the oxide layer cracks and/or develops porosity due to Pilling-Bedworth stresses developed during the metal-to-oxide transformation [73,74]. Once cracking occurs, new diffusion avenues open up for hydrogen to reach the interface and diffusion is no longer seen as a rate limiting step. Our model applies only up to this transition point but not beyond.
Schematic diagram (not to scale) of the geometry considered for the hydrogen penetration and hydride model developed in this work. x is the depth variable, s is the thickness of the oxide scale, and L is the total thickness of the clad. The chemical processes occurring at each interface are shown for reference.

Diffusion Model of Hydrogen in ZrO 2
The goal of this part of the model is to determine the hydrogen concentration at the metal/oxide interface as a function of time. For this, a generalized drift-diffusion equation is solved: This equation includes the following contributions: • The first term is standard Fickian diffusion in the presence of a concentration gradient.
• The second term is the so-called thermo-migration contribution, which depends on the temperature gradient and where U H is the activation energy for diffusion. The convention is for interstitial solutes to move in the direction opposing the gradient, i.e., a 'negative' drift contribution in the equation. • The third term represents electro-migration, where q is the charge of the diffusing species (+1 for protons), and φ is the electrical potential, which can be determined by solving Poisson's equation: where ρ is the charge density and ε is the dielectric permittivity.
Consistent with our previous work [69] and other studies [61], we assume the existence of a charge gradient across the oxide layer that originates from the onset of an electron density profile [75]. As well, in this work we consider autoclave conditions and thus neglect the thermomigration contribution.
Equation (3) is solved in one dimension (x) using the finite difference model with the following dynamic boundary conditions: where C 0 is the amount of oxygen (per unit volume) absorbed into the clad to form Zr oxide.The first condition trivially states that the hydrogen content in the clad at the beginning of time is equal to zero, while the second one prescribes the flux of hydrogen at the water/oxide interface. This condition is time-varying as indicated by the growth rate of the oxide layer,ṡ. As well, it depends on the H pickup fraction, f H , which albeit may also be time dependent [62,64], we fix at 15% for the remainder of this work. The factor of '2' represents the fact that there are two atoms of hydrogen per oxygen atom available to penetrate the clad. Under homogeneous oxide formation conditions, C 0 ≈ 2ρ Zr , with ρ Zr the Zr atomic density. Expressions for s(t) have been provided in our previous study for a number of nuclear-grade Zr alloys [69]. In general, s(t) = at n such that the growth rate can be directly expressed as a values range between 0.33 (pure Zr) and 0.37 (Zirc-4), while n = 0.34 in both cases. These values give s in microns when t is entered in days. With this,ṡ ≈ 0.11t −0.66 (microns per day).

Stochastic Cluster Dynamics Model with Spatial Resolution
Here we use the stochastic cluster dynamics method (SCD) [76] to perform all simulations. SCD is a stochastic variant of the mean-field rate theory technique, alternative to the standard implementations based on ordinary differential equation (ODE) systems, that eliminates the need to solve exceedingly large sets of ODEs and relies instead on sparse stochastic sampling from the underlying kinetic master equation [76]. Rather than dealing with continuously varying defect concentrations C i in an infinite volume, SCD evolves an integer-valued defect population N i in a finite material volume V, thus limiting the number of 'active' ODEs at any given moment. Mathematically, SCD recasts the standard ODE system: into stochastic equations of the form: The set {g,s,k} represents the reaction rates of 0th (insertion), 1st (thermal dissociation, annihilation at sinks), and 2nd (binary reactions) order kinetic processes taking place inside V, and is obtained directly from the standard coefficients {g, s, k} as: The value of V chosen must satisfy where is the maximum diffusion length l i of any species i in the system, defined as: Here, D i and R i are the diffusivity and the lifetime of a mobile species within V. The above expression is akin to the stability criterion in explicit finite difference models. From Equation (7), R i =s + ∑ jkij N j . The system of Equation (7) is then solved using the kinetic Monte Carlo (residence-time) algorithm by sampling from the set {g,s,k} with the correct probability and executing the selected events. Details of the microstructure such as dislocation densities and grain size are captured within SCD in the mean-field sense, i.e., in the form of effective sink strengths for hydrogen atoms. For example, the value of s ij in Equation (6) for dislocation and grain boundary defect sinks is: where ρ is the dislocation density and d is the grain size. SCD has been applied in a variety of scenarios not involving concentration gradients [76,77].
The SCD code has been developed in house and is available at: http://jmarian.bol.ucla.edu/packages/packages.html. However, Equation (6) must be expanded into a transport equation (i.e., a partial differential equation, or PDE) by adding a Fickian term of the type: where D i is the diffusivity of species i, and f (t; C 1 , C 2 , C 3 . . .) is used for simplicity to represent all of the terms in the r.h.s. of Equation (6). To cast Equation (11) into a stochastic form, the transport term must be converted to a reaction rate in the finite volume V. As several authors have shown, this can be readily done by applying the divergence theorem and approximating the gradient term in terms of the numbers of species in neighboring elements [78,79]. For a one-dimensional geometry such as that schematically shown in Figure 3, the Fickian term simply reduces to: where Greek superindices refer to neighboring elements, i is the cluster species, and l is the element size. When summed over all neighboring elements, this term then represents the rate of migration of species i from volume element α to β, which can be now added to the r.h.s. of Equation (7) and sampled stochastically as any other event using the residence-time algorithm. Figure 3. Schematic diagram of two volume elements of a 1D space discretization used to calculate spatial gradients within the stochastic cluster dynamics method (SCD). The superindex α refers to the physical element, while the subindex i refers to the cluster species. The full stochastic PDE system, written for a generic species i in volume element α, takes then the following form: Further, the model assumes the following: The only mobile species considered are hydrogen atoms. (ii) The source termg i only applies to element 0 (oxide/metal boundary) and is calculated from the hydrogen arrival flux calculated from the model in Section 2.2.
The only processes considered in the metal are: Next, we provide suitable expressions for each of the kinetic processes just listed.

H Atom Diffusion
The hydrogen diffusivity in both the oxide and the metal is assumed to follow an Arrhenius temperature dependence: where D 0 is the exponential pre-factor, e m is the migration energy, k is Boltzmann's constant, and the superscript α can refer to the oxide ('ox') or the metal ('m'). The diffusivity of H in Zircaloy-4 oxides has been recently measured by Tupin et al. [80], which give values of 2.5 × 10 −14 m·s −1 and 0.41 eV for D ox 0 and e ox m , respectively. Alloy composition, however, has been shown to have a significant impact on diffusion parameters. For example, values of e ox m = 0.55 and 1.0 eV have been reported for Zr-2.5%Nb and pure Zr, respectively, with D 0 numbers in as high as 1.1 × 10 −12 m·s −1 [81,82]. Here, we use the parameters for Zirc-4 given by Tupin et al. Similarly, the only mobile species in the metal is monoatomic hydrogen. The most widely used parameters for hydrogen diffusion in metal Zr, and Zircaloy-2 and -4 (D i in Equation (12)) are those by Kearns [83] in the 200-700 • C temperature range, with values of D m 0 = 7.90 × 10 −7 m·s −1 and e m m = 0.46 eV. Earlier literature on these measurements [25,84,85] reveals pre-factors ranging from 7.00 × 10 −8 to 4.15 × 10 −7 m·s −1 and migration energies between 0.3 and 0.5 eV, all in a similar temperature range. More recent experiments and molecular dynamics simulations are also consistent with these values [86,87].
The values chosen here for each case (diffusion in the oxide and in the metal) are given un Table 1.

Nucleation of Zr 2 H 3 Hydride
As shown in Figure 2, once hydrogen penetrates into the metal clad, the hydration reaction Zr + xH → ZrH x starts occurring. Although the formation of the δ-hydride is seen for a range of x values, here we assume a perfect stoichiometry of x=1.5. Consequently, the governing equilibrium constant for the reaction can be expressed as: However, it is more convenient to use an expression that is linear in the hydrogen concentration. From this, one can write the reaction rate as: which is simply a coagulation rate for two species -H and Zr-in the proportions indicated by the exponents of ρ Zr and N H . ∆E δ is the formation energy of a molecule of δ hydride (≈ 0.52 eV at 350 • C according to Blomqvist et al. [88]) and p(x) represents the thermodynamic probability for this reaction to occur, which can be directly extracted from the Zr-H phase diagram using the lever rule: where x TTS is the terminal thermal solubility at the temperature of interest, and x δ is the phase boundary. A phase diagram of the Zr-H system in the temperature and concentration region relevant to the present study is shown in Figure 4. By way of example, at 660 K (horizontal dashed line in the figure) x TTS is approximately 1.6% at. and x δ ≈60.0% at. Equation (15) ensures that p(x TTS ) = 0 and p(x δ ) = 1, i.e. the nucleation probability is zero at the phase boundary between the (α) and (α + δ) regions, and unity at the (α + δ) and (δ) boundary. This simple factor captures the thermodynamic propensity for the hydride reaction to take place, and thus ties the thermodynamics and kinetics of hydration together. r H and r Zr in Equation (14) are the interaction radii of H and Zr, respectively, their values given in Table 2. In the SCD calculations, the atomic fraction x is simply defined at any instant in time as:

Growth of Zr 2 H 3 Hydride
Once hydride nuclei appear in the clad, their growth is treated as a standard coagulation process in 3D with rate constant: where r δ is the interaction radius of the hydride clusters, and N (Zr 0.66n H n ) is the concentration of a hydride cluster containing n H atoms (which implies having 0.66n Zr atoms). It is assumed that hydride clusters are immobile. In accordance with previous works, hydrides grow as circular platelets whose size is directly related to the number of hydrogen monomers contained in it [56]: with Ω H and d the formation volume of hydrogen and the thickness of the platelet, respectively. As given in Table 2, here we use Ω H = 2.8 × 10 −3 nm 3 per atom [91], and d ≈ 0.28 nm [37]. The growth of hydride platelets is known to be highly directional, and influenced by stress and microstructure. Typically hydrides align themselves along the direction of the dominant axial stress components and grow preferentially in-plane on grain boundaries [27,36,37,53]. These details are not captured in our model at present.

Dissolution of Zr 2 H 3 Hydride
The last process considered in our model is the thermal dissolution of the hydrides, as Figure 4 shows, strictly speaking, hydrides are stable up to 550 • C (eutectoid temperature), although there is ample evidence of their decomposition at much lower temperatures, as well as the observation of thermal hysteresis during heating/cooling cycles [37,[92][93][94]. The dissociation rate is a first-order process that can be expressed as: where e b is the binding energy between a H monomer and a cluster containing n hydrogen atoms.
Here, we assume a capillary approximation for e b [55,56]: with e s being the heat of solution of H in the α-Zr matrix. This parameter has been found to be approximately 0.45 eV in electronic structure calculations [95,96], compared to 0.66 eV in experiments [97].

Metal/oxide Interface Motion
Finally, the motion of the interface must also be considered as a viable stochastic event. To turn the interface velocity, Equation (5), into an event rate, r i , one simply normalizes it by the interface thickness, s.
which results in the following expression for r i : This is added to the event catalog and sampled with the corresponding probability as given by Equation (18). As the equation shows, this is a time-dependent rate that reflects the nonlinear growth of the oxide layer with time. In the context of the SCD model, it implies that the one-dimensional mesh shown in Figure 3 must be dynamically updated with time because the physical dimensions of the simulation domain are dynamically changed. To our knowledge, this has not been attempted in any prior models of hydride formation and buildup.

Parameterization, Physical Dimensions, and Boundary Conditions
All the material constants used in the present model are given in Tables 1 and 2. External parameters representing the geometry and the boundary conditions are given in Table 3.

Results
The first two figures show results intended to set the stage for the SCD calculations. We begin with the time evolution of the hydrogen concentration at the oxide/metal interface. This results from solving the diffusion equation in the the oxide layer subjected to a moving boundary as explained in Section 2.2. Figure 5 shows the buildup of hydrogen up to the first 580 h. This represents a dynamic Dirichlet boundary condition for the spatially-resolved SCD calculations of hydride nucleation and growth in the metal substrate (g term in Equation (12)). Second, we track the sampling rate r i defined in Section 2.3.5 to confirm that it matches Equation (18). Figure 6 shows a comparison between both, indeed demonstrating their equivalency and confirming the correctness of its implementation in the code. The effect of this interface motion is that, over the course of the time scale covered in the SCD simulations, the oxide layer effectively sweeps over the first mesh element of the metal depth profile (recall that we assume that such sweep results in dissolution of the hydrides existing within that element at that point, and re-solution of the immobilized hydrogen in the metal). In practice, this allows us to subsequently discard the first spatial element of the 1D mesh. That is the reason why in the figures shown next the spatial range shown spans 800 (as opposed to the original 900) nm.   Next, we study the generation of hydride molecules in the metal layer as a function of time and depth. The results are shown in Figure 7a, which shows a histogram with the concentration of hydride molecules at several instants in time for each of the mesh elements of the metal region. As discussed in Section 2.3.2, the probability that a new hydride molecule will form depends primarily on the relative H concentration at the interface and the heat of formation of δ-hydride. With a probability per unit time k δ (Equation (14)), freely-diffusing H atoms are immobilized to form Zr2 /3 H molecules that act as incipient hydride nuclei. The concentrations of such nuclei are strongly depth-dependent, as shown in the figure, ranging over two orders of magnitude over the entire specimen thickness L of 900 nm. As well, the nucleation rate, i.e., the derivative of the evolution curves shown in Figure 7b (which display the same data as Figure 7a but plotted as a function of time), can be seen to decrease gradually in time across the entire depth profile.  Table 3, each element is 100-nm thick. Subsequent growth of these embryos occurs at a rate given by the combination of the rates of H-atom absorption (Equation (16)) and dissolution (Equation (17)), i.e., (k n − s n ), as shown in Figure 8. Rapid net growth is seen in the initial stages of hydridization close to the oxide/metal interface. However, these rates gradually abate both in time and with increasing depth until almost no net growth is observed, particularly at depths greater than 700 nm after 1.4 h of evolution.
The resulting hydride concentrations across the 900-nm metal layer at the end of the simulated time can be found in Figure 9a. As the graph indicates, the hydride number densities suffer almost a 100-fold decrease through the metal layer studied. In relative terms, these are large concentrations of small clusters, so it is to be expected that further time evolution of the hydride subpopulations will be dominated by growth, perhaps by way of some type of coarsening or ripening mechanism. The associated size distributions of the hydride clusters are shown in Figure 9b, where both the average and maximum cluster sizes are shown. We emphasize that, during the incipient nucleation of the hydrides, they grow as circular discs, and so the sizes simulated (≈ 50 nm or less), correspond to the regime prior to the acicular growth of the hydrides.

Discussion
Several of the most important features of the model presented here are: (i) consideration of a moving interface representing the growth of the oxide scale during operation in corrosive conditions; (ii) using a hydride nucleation criterion that is consistent with the thermodynamics of the Zr-H system; (iii) using a mean-field growth/dissolution model that respects; (iv) a completely physics-based parameterization based on calculated atomistic data. Some of these features were part of a comparable study [56], to which new ones have been added and existing ones augmented. All these features combined are the basis of a model that has been developed as an attempt to break the phenomenological vicious cycle in which models of materials degradation in nuclear environments are often found.
To study the nucleation of the hydride clusters, our method samples discrete kinetic processes defined by the corresponding energetics and thermodynamics. For example, hydride nucleation is simulated by considering the interplay between (i) aggregation, (ii) growth, and (iii) dissolution processes, which together determine the net nucleation and growth rates. Processes (i), (ii), and (iii) are embodied in Equations (14), (16) and (17), respectively. Each one of these processes is treated as a stochastic event sampled with the probabilities given by each respective rate. If the conditions are such that dissolution would dominate over nucleation, the clusters would never form. If growth dominates over nucleation, the clusters would grow bigger, etc. All the energetics are given by the parameters in each of those equations.
As is often the case, the price paid for an increased physical fidelity in the simulations is computational efficiency. For this reason, our simulations can only extend to times of several thousand seconds (<2 h), which is of course only representative of the initial stages of hydridation in Zr clad (and, of course, part of the pre-transition corrosion regime)). In these relatively short time scales (which are still orders of magnitude higher than what direct atomistic methods can cover), one can only claim to faithfully study the incipient nucleation phase of the hydride microstructure. In this sense, our results do not include important features of the Zr hydride particles such as their elongated shape and/or their orientation. Excellent recent examples of experimental characterization displaying all of these structure complexities exist now in the literature [99][100][101]. They can act, however, as a good springboard from which to connect to other methods such as phase field simulations [27,102,103], or orientation-dependent precipitation models [104,105]. Therefore, it is reasonable to assume that the time scale of the next phase of hydride formation/growth kinetics would be one dominated by coarsening/ripening, where population densities suffer a gradual decline at the expense of an increased average precipitate size. Thus, it is important to emphasize this aspect of the work: our results correspond to the incipient hydride nucleation and growth phase, before steady state populations are established. Steady state sizes and concentrations in corroded Zr specimens range from 100 nm to 1 µm [101] and ∼10 24 m −3 . On this aspect, it is also difficult to reconcile calculated H-atom diffusivities in the clad with almost cross-clad uniform hydride distributions observed experimentally [106]. Calculated migration energies suggest a much more sluggish diffusion in the metal, and screening of the clad interior by hydrides formed near the oxide metal interface, as seen in this study, compared to experimental results. While validation on the time and length scales covered in this work is always difficult, it is encouraging to see reasonable qualitative agreement with experimental studies, e.g., hydride precipitation completion fractions in ref. [36] (Figure 3) vs. Figure 7b in this paper. As well, our predictions for the size (long axis) of the precipitates in Figure 9b are in good agreement with in situ SEM observations [107].
As reviewed in the Introduction section, the formation of Zr hydrides in the metal clad is considered to be highly detrimental to reactor performance due to their embrittling effect. However, the high thermal stability of these hydride phases also makes them a matter of concern for reactor safety due to the potential for hydrogen storage and release during loss-of-coolant conditions and core meltdown. Palliative measures such as increasing the enthalpy of formation of δ-ZrH by microstructure tailoring [108], or by hindering H diffusion in Zr oxide by selective alloying in the clad [109,110], have been proposed for future candidate materials in novel nuclear fuel designs.

Conclusions
We end this paper with a list of the most important conclusions: • We have developed a spatially-resolved kinetic model of hydrogen transport/accumulation in Zr-metal clad. The model includes state-of-the-art hydride energetics data from atomistic calculations and is formulated as a stochastic version of the cluster dynamics method. Notably, boundary conditions are dynamically updated in time during the simulations, by accounting for oxide/metal interface motion due to the time-dependent growth of the oxide scale. • In doing so, our model is consistent with the oxidation in the clad, as well as with the equilibrium thermodynamics of the Zr-H system. • As most cluster dynamics models based on mean-field rate theory, our model does not capture the orientation dependence of elongated hydride platelets observed experimentally, and microstructural information such as grain sizes and dislocation densities is included only in an effective way. As such, our results are representative of the 'average' structure along the depth direction. • Our results show that high concentrations of small hydride nuclei form across the entire metal clad. This results in a very fine microstructure that sets the stage for the next kinetic phase, likely to be one of ripening and coarsening. • Gaps in our knowledge identified in this work include, among others: (i) how to model the H dissolved from hydrides swept by the growing oxide layer, (ii) how to reconcile existing H-atom diffusion energies with almost cross-clad uniform hydride distributions, and (iii) the reasons for the acicular (or capsular) growth of the precipitates are still not clear and, while such geometries can be adopted in the models, a physical approach that yields these geometric features is still lacking.