**1. Introduction**

Corrosion of metallic structural materials is a pervasive phenomenon in industry and technology [1–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–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–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–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–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 *δ* Zr2H3 displays good thermo-mechanical stability, it is also an exceedingly brittle phase [28–31] that can compromise the clad's mechanical integrity [18,31–33].

A key observation of the hydride microstructure is the elongated shape of the precipitates up to a few microns in length [34–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–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–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–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–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.

#### **2. Chemical Reaction Kinetics Model**

#### *2.1. 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:

$$2\text{HO}\_2 + 4\text{e}^- \rightarrow 2\text{O}^{2-} + 2\text{H}\_2\tag{1}$$

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–64]. This, known as the *pickup* fraction, sets the boundary condition for the adsorption of hydrogen at the clad's surface. Adsorbed H2 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:

$$\text{Zr} + \text{x} \text{H} \to \text{ZrH}\_x \tag{2}$$

where *x* is the atomic hydrogen concentration.

By way of illustration, Figure 1 shows representative hydridized microstructures in Zirc-4 and Zirconioum.

**Figure 1.** (**a**) Optical micrographs of hydride morphologies in Zircaloy-4. 'TD' and 'RD' indicate the tangential and radial directions in the clad (from Ref. [67], reproduced with permission of the International Union of Crystallography). (**b**) Electron micrograph detail of a needle-like Zr hydride (reproduced with permission from Ref. [68]).

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–72], here we consider a single phase (monoclinic) ZrO2 with thickness defined by the variable *<sup>s</sup>*(*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.

**Figure 2.** 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.

#### *2.2. 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:

$$\frac{\partial \mathbf{c\_H}}{\partial t} = \nabla \left( D\_\mathbf{H} \nabla \mathbf{c\_H} \right) - \frac{lI\_\mathbf{H} D\_\mathbf{H}}{kT^2} \nabla c\_i \nabla T + \frac{qD\_\mathbf{H}}{kT} \nabla \left( \mathbf{c\_H} \nabla \phi \right) \tag{3}$$

This equation includes the following contributions:


$$
\nabla^2 \phi = -\frac{\rho}{\varepsilon} \tag{4}
$$

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:

$$\begin{aligned} c\_{\rm H}(\mathbf{x},0) &= \, ^0 \mathbf{I} \\ J\_{\rm H}(0,t) &= \, ^0 D\_{\rm H} \frac{\partial c\_{\rm H}(0,t)}{\partial \mathbf{x}} = 2f\_{\rm H} \mathbf{C}\_0 \frac{\partial \mathbf{s}}{\partial t} \end{aligned}$$

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, *s*˙. 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<sup>n</sup>* such that the growth rate can be directly expressed as

$$\frac{\partial \mathbf{s}}{\partial t} = a m t^{n-1} \tag{5}$$

*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, *s*˙ ≈ 0.11*t* −0.66 (microns per day).

#### *2.3. 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 *Ci* in an infinite volume, SCD evolves an integer-valued defect population *Ni* in a finite material volume *V*, thus limiting the number of 'active' ODEs at any given moment. Mathematically, SCD recasts the standard ODE system:

$$\frac{d\mathbf{C}\_i}{dt} = \mathbf{g}\_i - \sum\_i \mathbf{s}\_{ij}\mathbf{C}\_i + \sum\_i \mathbf{s}\_{ji}\mathbf{C}\_j - \sum\_{i,j} k\_{ij}\mathbf{C}\_j\mathbf{C}\_i + \sum\_{j,k} k\_{jk}\mathbf{C}\_i\mathbf{C}\_k \tag{6}$$

into stochastic equations of the form:

$$\frac{dN\_i}{dt} = \vec{g}\_i - \sum\_i \vec{s}\_{ij} N\_i + \sum\_i \vec{s}\_{ji} N\_j - \sum\_{i,j} \tilde{k}\_{ij} N\_j N\_i + \sum\_{j,k} \tilde{k}\_{jk} N\_i N\_k \tag{7}$$

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:

$$\mathfrak{g} \equiv \mathfrak{g}V, \mathfrak{s} \equiv \mathfrak{s}, \,\,\check{k} \equiv kV^{-1}$$

The value of *V* chosen must satisfy

> *V*1 3 > -

where - is the maximum diffusion length *li* of any species *i* in the system, defined as:

$$\ell = \max\_{i} \{ l\_{i} \} \tag{8}$$

$$l\_i = \sqrt{\frac{D\_i}{R\_i}}\tag{9}$$

Here, *Di* and *Ri* 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), *Ri* = *s*˜ + ∑*j* ˜ *kijNj*. 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 *sij* in Equation (6) for dislocation and grain boundary defect sinks is:

$$s\_{i\bar{j}} = s\_d + s\_{\bar{\varsigma}}b = \rho + \frac{6\sqrt{\rho}}{d} \tag{10}$$

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:

$$\frac{d\mathbb{C}\_i}{dt} = \nabla \left( D\_i \cdot \nabla \mathbb{C}\_i \right) + f\left( t; \mathbb{C}\_1, \mathbb{C}\_2, \mathbb{C}\_3 \dots \right) \tag{11}$$

where *Di* is the diffusivity of species *i*, and *f*(*<sup>t</sup>*; *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:

$$D\_i \frac{N\_i^a - N\_i^\beta}{I^2}$$

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:

$$\frac{dN\_i^a}{dt} = D\_i \sum\_{\beta} \frac{N\_i^a - N\_i^{\beta}}{l^2} + \tilde{g}\_i - \sum\_i \tilde{s}\_{i\bar{j}} N\_i^a + \sum\_i \tilde{s}\_{\bar{j}i} N\_{\bar{j}}^a - \sum\_{i,j} \tilde{k}\_{i\bar{j}} N\_{\bar{j}}^a N\_i^a + \sum\_{\bar{j},\bar{k}} \tilde{k}\_{\bar{j}k} N\_i^a N\_k^a \tag{12}$$

Further, the model assumes the following:

	- (a) H diffusion
	- (b) Immobilization of H atoms through formation of Zr2H3 molecules (equivalent to nucleation of hydride platelets).
	- (c) Growth of Zr2H3 clusters.
	- (d) Thermal dissolution of Zr2H3 clusters.

Next, we provide suitable expressions for each of the kinetic processes just listed.

#### 2.3.1. H Atom Diffusion

The hydrogen diffusivity in both the oxide and the metal is assumed to follow an Arrhenius temperature dependence:

$$D\_{\rm H}^{a}(T) = D\_{0}^{a} \exp\left\{ \left( -\frac{e\_{\rm m}^{a}}{kT} \right) \right\} \tag{13}$$

where *D*0 is the exponential pre-factor, *em* 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−<sup>14</sup> <sup>m</sup>·s<sup>−</sup><sup>1</sup> and 0.41 eV for *D*ox0 and *<sup>e</sup>*ox*m* , respectively. Alloy composition, however, has been shown to have a significant impact on diffusion parameters. For example, values of *<sup>e</sup>*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−<sup>12</sup> <sup>m</sup>·s<sup>−</sup><sup>1</sup> [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 (*Di* in Equation (12)) are those by Kearns [83] in the 200–700 ◦C temperature range, with values of *D*m0 = 7.90 × 10−<sup>7</sup> <sup>m</sup>·s<sup>−</sup><sup>1</sup> and *<sup>e</sup>*m*m* = 0.46 eV. Earlier literature on these measurements [25,84,85] reveals pre-factors ranging from 7.00 × 10−<sup>8</sup> to 4.15 × 10−<sup>7</sup> <sup>m</sup>·s<sup>−</sup><sup>1</sup> 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.

#### 2.3.2. Nucleation of Zr2H3 Hydride

As shown in Figure 2, once hydrogen penetrates into the metal clad, the hydration reaction Zr + *x*H → 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:

$$\mathbf{K}\_{\delta} = \frac{\mathbf{[ZrH\_{1.5}]}}{\mathbf{[Zr]}\,\mathrm{[H]}^{1.5}}$$

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:

$$k\_{\delta} = 4\pi \left( r\_{\text{H}} + r\_{\text{Zr}} \right) \left( V^{-\frac{1}{3}} \rho\_{\text{Zr}}^{\frac{2}{3}} \right) D\_{\text{H}} N\_{\text{H}} p(\mathbf{x}) \exp\left( -\frac{\Delta E\_{\delta}}{kT} \right) \tag{14}$$

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:

$$p(\mathbf{x}) = \frac{\mathbf{x} - \mathbf{x}\_{\text{TTS}}}{\mathbf{x}\_{\delta} - \mathbf{x}\_{\text{TTS}}} \tag{15}$$

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 *<sup>x</sup>δ*<sup>≈</sup>60.0% at.

**Figure 4.** Phase diagram of the Zr-H system in the temperature and concentration region relevant to the present study (adapted from several sources [21,89,90]).

Equation (15) ensures that *p*(*<sup>x</sup>*TTS) = 0 and *p*(*<sup>x</sup>δ*) = 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:

$$\mathbf{x} = \frac{N\_{\rm H}}{N\_{\rm H} + \rho\_{\rm Zr}V}$$

#### 2.3.3. Growth of Zr2H3 Hydride

Once hydride nuclei appear in the clad, their growth is treated as a standard coagulation process in 3D with rate constant:

$$k\_{\rm II} = 4\pi V^{-1} \left( r\_{\rm II} + r\_{\delta}(n) \right) D\_{\rm II} N\_{\rm II} N\_{\rm (Zv\_{0.6\delta\alpha} \text{H}\_{\rm w})} \exp\left( -\frac{\Delta E\_{\delta}}{kT} \right) \tag{16}$$

where *rδ* is the interaction radius of the hydride clusters, and *<sup>N</sup>*(Zr0.66*n*H*n*) is the concentration of a hydride cluster containing *n* H atoms (which implies having 0.66*n* 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]:

$$r\_{\delta}(n) = \sqrt{\frac{n\Omega\_{\mathbb{H}}}{\pi d}}$$

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−<sup>3</sup> nm<sup>3</sup> 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.

#### 2.3.4. Dissolution of Zr2H3 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–94]. The dissociation rate is a first-order process that can be expressed as:

$$s\_{\rm tr} = 4\pi r\_{\delta}(n) D\_{\rm H} N\_{(\rm Zr\_{0.66n} \,\text{H}\_{\rm v})} \exp\left\{ \left( -\frac{c\_{\delta}(n)}{kT} \right) \right\} \tag{17}$$

where *eb* is the binding energy between a H monomer and a cluster containing *n* hydrogen atoms. Here, we assume a capillary approximation for *eb* [55,56]:

$$e\_b(n) = e\_s - 0.44\left[n^{\frac{2}{3}} - (n-1)^{\frac{2}{3}}\right]$$

with *es* 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].

#### 2.3.5. 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, *ri*, one simply normalizes it by the interface thickness, *s*.

$$r\_i = \left(\frac{1}{s}\right)\frac{ds}{dt} = \frac{ant^{n-1}}{at^n} = \frac{n}{t}$$

which results in the following expression for *ri*:

$$r\_l = 0.34t^{-1} \tag{18}$$

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.

#### *2.4. 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.


**Table 1.** Zr-H energetics used in the model with the respective source.

**Table 2.** Physical constants for the Zr-H system employed here. In actuality, the interaction radii of Zr and H atoms are extended by a distance equal to the Burgers vector a in *α*-Zr, which is equal to 3.23 Å.


**Table 3.** Numerical parameters used in the model.

