A computational model of anaerobic reactions in metallic multilayered systems with an equimolar composition of zirconium and aluminum is developed. The reduced reaction formalism of M. Salloum and O. M. Knio, Combust. Flame **157**(2): 288–295 (2010) is adopted. Attention is focused on quantifying intermixing rates based on experimental measurements of uniform ignition as well as measurements of self-propagating front velocities. Estimates of atomic diffusivity are first obtained based on a regression analysis. A more elaborate Bayesian inference formalism is then applied in order to assess the impact of uncertainties in the measurements, potential discrepancies between predictions and observations, as well as the sensitivity of predictions to inferred parameters. Intermixing rates are correlated in terms of a composite Arrhenius law, which exhibits a discontinuity around the Al melting temperature. Analysis of the predictions indicates that Arrhenius parameters inferred for the low-temperature branch lie within a tight range, whereas the parameters of the high-temperature branch are characterized by higher uncertainty. The latter is affected by scatter in the experimental measurements, and the limited range of bilayers where observations are available. For both branches, the predictions exhibit higher sensitivity to the activation energy than the pre-exponent, whose posteriors are highly correlated.

## I. INTRODUCTION

Reactive nanolaminate foils, tens of microns thick, are metallic multilayered systems made up of hundreds of alternating layers of elements that mix exothermically (Fig. 1). Advanced fabrication techniques for such systems include physical vapor deposition whereby the individual layers are sputter deposited on a substrate.^{2–9} Typically, individual layers have thicknesses in the order of tens of nanometers.

Sputtered multilayered nanolaminates offer an interesting setting for studying reaction properties owing to their uniform layering, small atomic diffusion distance, rapid heating, and large concentration gradients. Numerous experimental investigations have focused on the self-propagating nature of intermetallic formation reactions and associated phase transformations.^{2,10–12} The large amount of localized heat due to exothermic mixing in the multilayers has also motivated a number of practical applications, such as welding, brazing, and soldering.^{13–20}

Recently, Joress *et al.*^{21} conducted experiments on Zr-Al nanolaminate foils with varying compositions and microstructures and focused on investigating the oxidation properties of the nanolaminate following the completion of the anaerobic reaction. Their study revealed that for a nanolaminate foil comprised of zirconium and aluminum in a molar ratio of 1:1, reacted in air, the oxidation reaction can be sustained for over a second. A simplified model was subsequently developed to characterize the oxidation behavior in the equimolar Zr-Al system.^{22} However, despite the availability of a body of experimental evidence concerning anaerobic formation reactions in the Zr-Al system (see, e.g., Refs. 23 and references therein), a computational model that can describe the evolution of anaerobic formation reactions in the Zr-Al system is not available. The present work focuses on the development of such a model.

While numerous models for the simulation of transient reactions in reactive multilayered systems have been developed and refined over the past three decades,^{1,3,24–31} the adaptation of existing models to characterize anaerobic reactions in Zr-Al system is limited by a lack of suitable characterization of atomic intermixing rates, and the dependence of these rates on the local temperature. To overcome this hurdle, we rely on recent experimental observations for two reaction regimes. The first concerns essentially homogeneous reactions initiated by means of a current pulse, and the second concerns measurements of the velocity of self-propagating reaction fronts. Following the methodology outlined in Ref. 32, we exploit experimental temperature measurements during uniform ignition of the foil to estimate the atomic diffusivity in a low temperature range. The uniform ignition ranges approximately from the ignition temperature (below which the temperature of the sample is primarily affected by the external stimulus) to approximately the melting point of Al (above which the temperature rise is too fast for accurate information to be drawn). To calibrate intermixing rates at higher temperatures, extending beyond the melting of Al, we rely on experimental measurements of reaction velocity using foils with varying bilayer thicknesses.^{33} The comparison of model predictions and experimental data enables us to determine diffusivity estimates in a suitably wide temperature range.

This paper is organized as follows. In Sec. II, we provide a brief description of the experimental measurements used to support our inference of key parameters needed in the construction of the model for anaerobic Zr-Al reactions in nanolaminates. The formulation of the computational model is then outlined in Sec. III. Section IV then discusses two approaches for inferring the intermixing rates based on the experimental measurements. The first is based on a simple regression methodology that aims at minimizing the square error between model predictions and experimental observations, whereas the second is based on a more elaborate Bayesian inference formalism that enables us to take into account measurement uncertainties and to quantify the discrepancy between the measurements and the corresponding model predictions. Computed results are discussed and analyzed in Sec. V, and major conclusions are highlighted in Sec. VI.

## II. DATA SOURCES

As mentioned in the introduction, our experimental data sources consist of (i) hot plate ignition experiments, (ii) homogeneous ignition experiments, and (iii) measurements of the velocity of self-propagating reaction fronts.

The hot plate ignition experiments consisted in dropping small fragments of reactive Zr-Al foils on a hot plate, and measuring the plate temperature using a thermocouple pressed against its surface. Precisely, 15 different foil fragments were used to account for the variability in the respective ignition points depending upon their size and orientation on the plate surface. The overall ignition point corresponds to the minimum temperature at which all the fragments were ignited. The ignition point is achieved within an accuracy of ±5° by systematically increasing or decreasing the plate surface temperature using an iterative process. Tests were conducted for bilayer thicknesses ranging from about 50 nm to 90 nm and the data recorded for the case of a 20 *μ*m thick nanolaminate foil are illustrated in Fig. 2. Following the analysis developed by Fritz *et al.*,^{34} the ignition temperatures are exploited to estimate the activation energy of atomic diffusion in the low temperature regime. For the equimolar Zr-Al system, it is estimated to be approximately 53 kJ/mol, accurate within ±1.3 kJ/mol. This initial estimate of the activation energy is used in a regression analysis aimed at determining the Arrhenius correlation for atomic diffusivity at low temperatures.

To infer the pre-exponent, we rely on ignition experiments in which the Zr-Al multilayered foils are ignited using a current pulse of finite duration. The process involves holding a Zr-Al foil between two copper electrodes connected to a current-pulse generator. The experimental set-up is the same as that described in Ref. 33. As a result of uniform ohmic heating of the foil, the reaction evolves in an essentially homogeneous fashion. The pulse duration is optimized to the minimum value that enables ignition in the reactive multilayers. A one-color optical pyrometer, focused at the center of the foil, is used to measure the foil surface temperature. The pyrometer calibration error is ±100 K. The experimental observables thus consist in the recorded time evolution of the electrical power input into the foil, and the foil surface temperature. Variations in foil surface emissivity are not expected to have a significant impact on the measurements, due to the low temperatures during ignition and our focus on a temperature range bounded above by the melting temperature of Al. As further discussed in Sec. V below, homogeneous ignition profiles obtained for Zr-Al multilayered foils with *λ* = 68 nm and 72 nm are made available.

Our experimental database also includes reaction velocity measurements for Zr-Al multilayered foils with four different bilayer thicknesses, namely, *λ* = 64 nm, 79 nm, 81 nm, and 87 nm. To achieve this, a given foil is secured firmly between glass slides such that one end of the foil is protruding. The glass slides prevent further oxidation of the foil in the presence of ambient oxygen. An assembly of five optical fibers with known spacing is held above the foil and a self-propagating reaction is initiated by applying an electric spark at the protruding end. As the reaction front propagates along the length of the foil, velocity is estimated using the time taken by the front to pass each fiber in the arrangement. As discussed in Ref. 32, the velocity-bilayer data pairs are used to infer atomic diffusion parameters in the high temperature range.

## III. MODEL FORMULATION

In this section, we provide brief outlines of the models used for the simulation of uniform ignition, and transient reaction fronts in equimolar Zr-Al nanolaminates. As illustrated in Fig. 1, the multilayered system is assumed to have a uniform microstructure, with an alternating arrangement of flat layers of Zr and Al. For a 1:1 molar ratio, the thickness ratio of individual Al and Zr layer is given by

where *ρ* denotes the density, *M* is the atomic mass, and the superscripts refer to individual elements. For compositionally pure layers, *γ* = 1.4. The elemental layers are assumed to be separated by a premixed region of thickness, 2*w* ( ≈ 1.6 nm), that forms during the deposition process.

The evolution of the reactions in the Zr-Al system is described in terms of continuum models that are based on the assumptions that the atomic intermixing can be described in terms of Arrhenius rate expressions, and that volumetric changes during the reaction can be ignored. For uniform ignition experiments, the reaction is assumed to proceed uniformly throughout the foil. This enables us to adopt a simplified formulation consisting of ordinary differential equations that govern the evolution of the enthalpy and the degree of mixing. For self-propagating fronts, a space-time formulation is adopted that also accounts for thermal transport in the layered medium. The formulations of the ignition model and transient front propagation model are discussed in Secs. III A and III B, respectively.

### A. Homogeneous reaction model

The development of the homogeneous reaction model is specifically tailored to uniform ignition experiments in which spatial variations within the reacting multilayers are negligible. The approach based on a lumped parameter formulation of the energy equation is used, namely:^{33}

where $H\u0307$ denotes rate of change of volume-averaged foil enthalpy, whereas $P\u0307,\u2009Q\u0307$, and $L\u0307$, respectively, denote the power input due to Joule heating by the current pulse, the rate of heat release due to chemical reaction, and the rate of heat loss by means of natural convection to the surroundings, all on a volumetric basis. Note that heat loss by radiation has been neglected in the ignition model^{34} in light of the limited range of the temperature measurements, and that heat losses to the surroundings are considered to be dominated by natural convection.^{32} Conductive heat loss to the copper blocks (see experimental set-up in Ref. 33) is not expected to affect the measurements significantly as the pyrometer is focused at the center of the foil and since the length of the foil remains significantly larger than the thermal penetration depth during the short duration of the experiment.

The power input due to Joule heating per unit volume is given by

where $V\u2032$ and *I* are the applied voltage and current across the foil, respectively, and *V _{foil}* denotes the foil volume.

The rate of heat loss to the surroundings per unit of foil volume due to natural convection is expressed as

where *d* is foil thickness, *T*_{0} is the ambient temperature

is the associated heat transfer coefficient,^{35,36} and *L _{c}* is a characteristic length defined as the ratio of the foil surface area and the perimeter of the cross-section.

The heat release rate is related to the rate of atomic intermixing. The latter is described in terms of a conserved scalar, *C*,^{1,27,28} assumed to be governed by a Fickian diffusion process, namely, according to

where *D*(*T*) is the atomic diffusivity. The scalar *C*, also referred to as concentration, quantifies the degree of atomic mixing; it assumes a value of 1 for pure Al, −1 for pure Zr, and 0 for the fully mixed intermetallic. The atomic diffusivity is assumed to follow an Arrhenius dependence on temperature

where *D*_{0} is the pre-exponent and *E _{a}* is the activation energy.

Following the reduced formalism developed in Ref. 1, the PDE in Eq. (5) is transformed into an ODE by introducing the stretched time variable

This enables us to transform Eq. (5) into a canonical form

where $\xi \u2261y\delta $ is a normalized spatial variable. The solution of Eq. (8) is computed in a pre-processing step and made available to the computations.^{1,27,28}

The volume-averaged heat release rate, $Q\u0307$, is related to the evolution of intermixing through^{1,27,28,37}

where $\rho Cp\xaf$ is the mean thermal capacity of the system (J/m^{3}/K) and $\rho \xaf\u2261\u2009(\rho Al+\gamma \rho Zr)/(1+\gamma )$ is the average density. The temperature-dependent molar heat capacities of Zr and Al, depicted in Fig. 3, are used to estimate $Cp\xaf$. Note that the mean-square concentration, $C2\xaf(\tau )$, can also be derived from the canonical solution, as outlined in Ref. 27. Also note that $\Delta Tf\u2261\Delta Hrxn/\rho Cp\xaf$ estimates the temperature rise under adiabatic conditions and is estimated from the heat of mixing, $\Delta Hrxn\u224890\u2009kJ$ per mole-atom.^{38} For the Zr-Al system, $\Delta Tf$ is estimated to be 1758 K.

Substituting individual expressions for $P\u0307,\u2009Q\u0307$, and $L\u0307$ in Eq. (1) results in

Note that $H[tp\u2212t]$ is the Heaviside unit step function and *t _{p}* is the duration of current pulse.

The instantaneous temperature is recovered from the instantaneous enthalpy estimates obtained from the numerical integration of Eq. (10) using the following relationship that accounts for the heats of fusion, adapted from:^{25}

For the Zr-Al system, the individual enthalpies, $H0\u2192H6$, are defined by

where *β* = $C\xaf/(1+\gamma )$ is the fraction of unmixed, pure Al; $HfAl$ = $\rho AlMAlHfAl$, $HfZr$ = $\rho ZrMZrHfZr$, and $HfZrAl$ = $\rho \xafMZrAlHfZrAl$. $HfAl,HfZr$, and $HfZrAl$ denote individual heats of fusion in J/m^{3} for Al, Zr, and ZrAl, respectively. We use $HfAl$ = 10.7 kJ/mol,^{39}^{,}$HfZr$ = 19.25 kJ/mole,^{40} and $HfZrAl$ = 10.75 kJ/(mol-atom),^{41} together with the following values for the melting temperatures of Al, Zr, and ZrAl, $TmAl=933\u2009K$, $TmZr=2125\u2009K$, and $TmZrAl=1758\u2009K$. Note that in the range of temperature measurements in ignition experiments, the temperature falls below the melting temperature of Al, so that the enthalpy falls below *H*_{1}. Consequently, the simulations in this regime are not affected by phase change effects.

### B. Distributed reaction model

The modeling of self-propagating reaction fronts relies on a similar methodology, namely, based on coupling the energy conservation equation with the evolution equation for the conserved scalar. In the present formulation, however, thermal transport must be accounted for, and so the conservation of energy over a region (or computational cell), Ω, is expressed in integral form as^{1,27,28,37}

where ** q** is the heat flux, and the region-averaged chemical heat release term, $\u2202Q\u2202t$, is estimated using Eq. (9). The conduction heat flux is assumed to be given by Fourier's law:

where $k\xaf$ is mean thermal conductivity, obtained using the temperature-dependent thermal conductivities for Zr:^{42}

and Al:^{43}

Specifically, it is estimated using

### C. Numerical simulation

For the case of a homogeneous reaction, the enthalpy and concentration fields degenerate into scalars, and the state is consequently described by the two-dimensional state vector, (*H*, *τ*). The evolution of this state vector is obtained by numerically integrating Eqs. (7) and (1) using a finite difference methodology. The integration scheme developed in Ref. 27 is used for this purpose. A small integration time step, Δ*t* = 100 ns is used to ensure a stable and accurate solution.

Simulation of the non-homogeneous reaction model is based on a finite difference methodology. The computational domain consists of a two-dimensional rectangular region, described using a uniform grid with mesh size Δ*x* = Δ*y* = 1 *μ*m. Adiabatic boundary conditions are imposed, and the coupled system consisting of Eqs. (13) and (7) is integrated over each of the computational cells. A second-order central-difference is used to approximate the conduction heat flux. The numerical scheme developed in Ref. 28 is used in the simulation, with an integration time step, Δ*t* = 50 ns.

## IV. INFERENCE OF DIFFUSIVITY

As discussed earlier, our model formulation for Zr-Al nanolaminates lacks key parameters describing the dependence of intermixing rates on temperature. To address this hurdle, we follow the methodology recently developed in Ref. 32 for the Ni-Al system. Briefly, the analysis in Ref. 32 revealed that a single Arrhenius correlation for the atomic diffusivity cannot provide suitable estimates of intermixing rates in a sufficiently wide temperature range. This problem was overcome by estimating Arrhenius branches for temperatures falling below and above the melting temperature of Al. Specifically, parameters for the first branch were calibrated using homogeneous ignition measurements, whereas parameters for the second branch were estimated based on reaction front velocity measurements.

We shall adopt the methodology developed in Ref. 32 to the Zr-Al system. In the implementation, we explore two approaches for extracting the desired coefficients from measured data. The first is based on a regression analysis, whereas the second relies on a more elaborate Bayesian inference formalism. Both approaches are outlined below.

### A. Regression analysis

The regression approach is based on an optimization procedure that aims at minimizing the root mean square (RMS) error between the numerical predictions and the data. For the low temperature branch, the procedure essentially consists in fixing the value of the activation energy based on the hot plate ignition experiments, and calibrating the pre-exponent, *D*_{0}, using the measured temperature profiles obtained from the uniform ignition experiments.

For the high temperature branch, our implementation aims at minimizing the RMS error between reaction velocity measurements for different bilayer thicknesses and the corresponding model predictions. The procedure iteratively updates both the pre-exponent, *D*_{0}, and the activation energy, *E _{a}*, and yields an optimal pair such that the RMS error falls below a set tolerance limit. In this case, the approach is the same as that outlined in Ref. 32.

### B. Bayesian inference

In contrast to the regression approach, which generally identifies optimal parameter values, the Bayesian inference approach exploits the measured data in determining the full posterior distribution of the parameters being calibrated. This offers the advantage of readily yielding quantitative estimates of the information gained from the data, the impact of measurement uncertainties, and consequently, the level of confidence in calibrated parameters.

In the present implementation, the inference specifically aims at determining the posterior distributions of the diffusion parameters, *D*_{0} and *E _{a}*. Assuming reasonable prior distributions for these parameters, the posterior is determined through the application of Bayes' rule:

^{44,45}

where $X\u2261(D0,Ea,\sigma 2)$ is the extended parameter vector, *σ*^{2} is the hyper-parameter (defined below), $P(X|Ti)$ is the *posterior* of $X$ given the observations *T _{i}*, and $L(Ti|X)$ is the likelihood of observing the data given a parameter $X$.

For the uniform ignition experiments, the data *T _{i}* consist of the temperature measurements at times

*t*from the start of the current pulse. In this case, we assume a likelihood of the form

_{i}where *N* is the number of data pairs. The quantity $(dTdt)$ represents mean slope of the temperature versus time profile. Note that the discrepancy between the model and experimental data, $(ti\u2212Mi)$, is assumed to be normally distributed with zero mean and variance, $\sigma T2$. Here, *M _{i}* is the time at which the simulated temperature coincides with the experimentally recorded temperature,

*T*. Note that the variance is not known

_{i}*a priori*and is consequently treated as hyper-parameter to be determined.

When calibrating the diffusion parameters in the high-temperature Arrhenius branch, the likelihood function is assumed to be given by

where *v _{i}* denotes the measured front velocity,

*M*denotes the computed front velocity, $\sigma v2$ is the hyper-parameter, and

_{i}*N*is the total number of velocity measurements.

The prior $P(X)$ is assumed to be given by the product of individual probabilities, $P(D0),\u2009P(Ea)$, and $P(\sigma 2)$. This is a reasonable assumption since no prior information is given concerning either the physical parameters or the hyper-parameter. The prior distributions of *D*_{0} and *E _{a}* are assumed to be uniform over a finite width interval. We use our knowledge from the regression approach to select appropriate bounds. For the hyper-parameter, an uninformative distribution is assumed based on Jeffrey's prior:

^{46}

To construct the posterior distribution, we rely on a Markov Chain Monte Carlo (MCMC) method and conduct the simulations using the adaptive Metropolis algorithm.^{47} This requires that a sufficiently large number of samples be generated to suitably characterize the likelihood, and consequently the posterior. To mitigate the costs of repeated model evaluations (at different values of the physical parameters),^{48,49} a surrogate is first constructed, which is then sampled in lieu of the actual model in estimating the likelihood. As discussed further below, we rely on a polynomial chaos (PC) methodology^{50–56} for the purpose of constructing suitable surrogates for the experimental observables.

### C. PC surrogate

The PC surrogate methodology is based on parameterizing the *N* uncertain inputs using canonical random variables, *ξ _{i}*, $i=1,\u2026,N$, and expressing the quantities of interest in terms of a truncated series expansion. For a generic quantity of interest (QoI),

*Q*, the surrogate is expanded according to

^{50–56}

where $\xi $ is the vector of uncertain parameters, *P* + 1 is the total number of terms retained in the expansion, and the $\Psi k$'s form an orthogonal polynomial basis in the space of functionals that are square integrable with respect to the probability measure characterizing the canonical random variables. Note that in the present implementation, the germ, $\xi =(\xi 1,\xi 2)$, where *ξ*_{1} and *ξ*_{2} are independent random variables uniformly distributed in the interval [−1, 1], and, respectively, used to parameterize the pre-exponent, *D*_{0}, and activation energy, *E _{a}*. Accordingly, the $\Psi k$'s are multidimensional Legendre polynomials acting on the germ, $\xi $.

^{51,56}As is customary in polynomial approximations, the PC basis is truncated by retaining polynomials with total order less than or equal to

*p*, leading to a basis of size $P+1=(N+p)!N!\u2009p!$.

^{56}

To determine the coefficients of the expansion, a non-intrusive spectral projection (NISP) approach is adopted.^{56} NISP essentially consists in exploiting the orthogonality of the basis functions by expressing the modes *Q _{k}* as integrals, and approximating the integrals using appropriate quadratures. In the implementations below, we rely on a fully tensored quadrature formula

^{56}with sufficiently high order to determine the coefficients accurately.

Finally, note that the availability of the PC representation enables us to readily assess the contribution of the individual parameters to the overall variability of the corresponding quantity. Suitable metrics are provided by the so-called Sobol sensitivity (global) indices.^{56–61} In the analysis below, we rely on the first-order indices and total sensitivity indices associated with the first and second components of the germ. The first-order sensitivity indices, $S1$ and $S2$, essentially quantify the direct contribution of *ξ*_{1} (*D*_{0}) and *ξ*_{2} (*E _{a}*) to the variance of the QoI,

*Q*, whereas the total sensitivity indices, $T1$ and $T2$, combine the direct contribution with that arising from mixed terms.

Note that for the present two-dimensional germ, simplified expressions for the Sobol indices can be obtained. Specifically, defining $K1$ to be the subset of the index set $L\u2261{1,\u2026,P}$, consisting of indices for which the corresponding polynomial is order 0 in *ξ*_{2}, and $K2$ to be the subset of $L$, consisting of indices for which the corresponding polynomial is order 0 in *ξ*_{1}. Then, the $Si$'s and $Ti$'s can be readily obtained from

and

respectively. As further discussed below, the Sobol indices can provide valuable insight into the impact of individual components to the variability in the selected QoI's, even prior to introducing (or collecting) the data or performing the inverse (calibration) analysis.

## V. RESULTS

We now focus on inferring the parameters of the atomic diffusivity correlations based on the temperatures measured during uniform ignition tests, and the reaction front velocity measured during self-propagation. We start by presenting results obtained using the regression analysis and then discuss results obtained using the Bayesian inference. We conclude with a brief discussion of some features of the propagation velocity versus bilayer curve, generated using calibrated diffusivity parameters.

### A. Regression analysis

We focus first on the low-temperature branch, and rely on the measured temperature versus time profiles to calibrate the pre-exponent in the atomic diffusivity law. As discussed earlier, the activation energy is held fixed, *E _{a}* = 53 kJ/mol, which coincides with the estimate obtained from the hot plate ignition experiments. The analysis is conducted for Zr-Al nanolaminates with

*λ*= 68 nm and

*λ*= 72 nm, for which experimental temperature profiles are available.

Figure 4 shows computed temperature profiles for different values of the pre-exponent. Included are plots obtained for both values of *λ*. The experimental temperature profiles are also shown for comparison. In each plot, the computed temperature profiles correspond to a value of *D*_{0} that minimizes the RMS difference between the model predictions and the experimental results (dashed lines), and to values above and below this optimum. The results indicate that for *D*_{0} ≈ 1 × 10^{–10} m^{2}/s, a very good agreement between model predictions and experimental results can be observed, and that there is close agreement between the optimal values of *D*_{0} for both bilayers. This is further examined below using the Bayesian inference computations. Also note that, owing to the limitations of the one-color pyrometer used to measure the temperature, no experimental data are available below a temperature of 473 K. However, this measurement limitation does not impact the analysis because in the range $T\u2009\u2a85\u2009500\u2009K$, the temperature rise is primarily governed by the applied current pulse. This is consistent with the observation that at low temperatures, the temperature rise is approximately linear.

To calibrate the intermixing parameters at higher temperatures, we exploit reaction velocity measurements for *λ* = 64 nm, 79 nm, 81 nm, and 87 nm. To this end, we adopt the two-parameter optimization algorithm developed by Alawieh *et al.*^{32} The optimization algorithm determines the optimal pair of *D*_{0} and *E _{a}* for $T\u2265TmAl$; intermixing rates at lower temperatures are simulated based on the optimized values provided above. The two-parameter optimization yields the following estimates for the high temperature branch,

*D*

_{0}= 3.13 × 10

^{–9}m

^{2}/s and

*E*= 55 kJ/mol.

_{a}To examine our predictions, we plot in Fig. 5, the computed and measured velocities for all four values *λ*. Shown are predictions obtained using the optimal *D*_{0} and *E _{a}* pair, as well as values lying on both sides of the optima. The results indicate that there is good agreement between computed and measured velocities at the optimal point. It is also seen that the velocity predictions increase sharply as

*D*

_{0}increases and decrease rapidly as

*D*

_{0}is decreased. The reverse trends can be observed as

*E*is varied from the optimum.

_{a}Combining the results of two calibration analyses yields a composite curve for the atomic diffusivity. Note that the values of the activation energy for both branches are comparable whereas the pre-exponent in the high-temperature branch is significantly larger. Thus, significant enhancement in intermixing rates around the melting of Al is primarily attributed to the pre-exponent of diffusivity. Though one would expect smaller activation energy for diffusion in the liquid state compared to the solid state, because *D*_{0} and *E _{a}* are inferred simultaneously from the velocity data, the fact that a significant variation in the activation energy across the melting of Al is not observed is simply an outcome of the analysis. A similar result was observed in the analysis of intermixing rates in Ni-Al multilayers,

^{32}which showed small differences between the activation energies of the low and high-temperature branches, but slightly more pronounced than presently observed for the Zr-Al system.

### B. Bayesian inference

We start by constructing PC surrogates that can be used to support and accelerate the Bayesian inference of diffusivity parameters. Specifically, we are interested in constructing functional representations of the dependence on the pre-exponent and activation energy of (a) the temperature–time profile in homogeneous ignition experiments, and (b) the velocity–bilayer curve for self-propagating fronts. In both cases, we assume uniform priors for *D*_{0} and *E _{a}*, and use previous knowledge from the regression analysis to define suitable ranges. For the low-temperature branch, we use $D0\u2208[0.62\xd710\u221210,1.20\xd710\u221210]\u2009m2/s$, and $Ea\u2208[47.83,52.87]\u2009kJ/mol$. For the high-temperature branch, the uniform priors are defined according to $D0\u2208[1.88\xd710\u22129,4.38\xd710\u22129]\u2009m2/s$ and $Ea\u2208[49.5,60.5]\u2009kJ/mol$. These distributions are parameterized in terms of independent, canonical random variables

*ξ*

_{1}and

*ξ*

_{2}that are uniformly distributed over [−1, 1]. We use a sixth-order fully tensored Gauss-Legendre quadrature rule to evaluate the PC coefficients. Thus, the 1-D rule has 7 quadrature points, leading to a total of 49 realizations for the two dimensions.

Figure 6(a) depicts the time-temperature response surface as function of the uncertain inputs, for the case of homogeneous ignition of a Zr-Al foil with *λ* = 68 nm. Specifically, we construct surrogate models of the quantity $t(T;\xi 1,\xi 2)$, which is defined as the time needed for the foil to reach a temperature *T*, given realizations of *D*_{0} and *E _{a}* corresponding to

*ξ*

_{1}and

*ξ*

_{2}. We illustrate in Figure 6(a), a two-dimensional surface corresponding to a selected value of temperature,

*T*= 600 K. Also shown are the predictions generated at the nodes of the Gauss-Legendre quadrature used to determine the PC coefficients. Surrogate models are also obtained for the reaction front velocity,

*V*, as function of the bilayer thickness, given realizations of the Arrhenius parameters of the high temperature branch. In all cases, the parameters of the low-temperature branch are held fixed; these correspond to the maximum a posteriori probability (MAP) estimates obtained from the Bayesian inference analysis using the ignition data (see discussion below). Thus, in this case, we obtain functional representations of the form $V(\lambda ;\xi 1,\xi 2)$, with

*ξ*

_{1}and

*ξ*

_{2}parameterizing the Arrhenius parameters of the high temperature branch. Again, a sixth-order Gauss-Legendre quadrature rule is used to build the velocity surrogate. Response surface based on the resulting PC expansion is shown in Fig. 6(b), for the 64 nm bilayer. Note that as expected, for a fixed value of the bilayer, the reaction velocity increases as the pre-exponent increases and the activation energy decreases. Also note that, consistent with earlier studies on reactive multilayered systems,

^{24–26,62}the reaction velocity generally increases as the bilayer decreases.

It is observed that the distribution of the individual model predictions is in close agreement with the response surfaces, indicating that the surrogate is a faithful representation of the model. To quantify this, we evaluate the relative *l*^{2} norm of the error, *E*, between the model realization and the surrogate at the quadrature nodes

where $R(\xi j)$ are model realizations at the quadrature nodes, and $\xi j$ and *w _{j}* are corresponding weights. For all temperatures, the results yield a value of

*E*, $O(10\u22124)$, which supports the observations made earlier. The surrogate model for reaction velocity is also observed to be in very good agreement with the individual realizations at the quadrature nodes, and the estimate of the relative error,

*E*, is $O(10\u22123)$, for all values of

*λ*.

Figure 7 shows probability density functions (PDFs) of the reaction time for different values of *T*. The PDFs are generated by a large number of Monte Carlo samples of the surrogates, and applying a kernel density estimation to the resulting ensemble. The curves indicate that as the temperature increases, the location of the peak shifts to the right, while its spread increases. This first trend is expected because more time from ignition is needed to reach higher temperature values, while the widening of the PDF reflects an increased variability at higher temperatures.

Figure 8 illustrates the first-order and the total sensitivity indices of the temperature and the reaction velocity to the uncertain inputs, respectively, in uniform ignition and self-propagation experiments. The results indicate that in both cases, the quantities of interest exhibit higher sensitivity to the activation energy than the pre-exponent. For the case of reaction times, one observes a noticeable difference between first-order and total sensitivity indices, highlighting the importance of mixed terms in the corresponding PC expansion. In addition, unlike the first-order index, the total sensitivity index exhibits a significant variation with time, reflecting an increase in the total sensitivity of the pre-exponent and a decrease in the total sensitivity of the activation energy. Thus, the impact of the pre-exponent on the temperature variability is higher for higher temperatures. On the other hand, for the self-propagation velocity, the first-order and total sensitivity indices are comparable, and reflect essentially the same values for all bilayers. The latter observation is consistent with the analysis in Ref. 31, which revealed that the mean concentration versus temperature curves for different bilayers essentially collapse on each other. Thus, for self-propagation, varying the bilayer essentially leads to a scaling of the solution. Consequently, in light of the analysis in Ref. 31, the present observation that the sensitivities to *D*_{0} and *E _{a}* in the self-propagation regime match for different bilayers is in fact expected.

Figures 9(a)–9(c), respectively, show sample chains for *D*_{0}, *E _{a}*, and the hyperparameter, $\sigma T2$. The chains extend for over 10

^{5}steps of the adaptive MCMC algorithm. In all cases, good mixing properties can be observed following the initial burn-in period, indicating that the statistics of the posterior are well captured; this is also confirmed through the decay of the autocorrelation of the signals (not shown). Figure 9(d) shows the joint posterior distribution of the pre-exponent and activation energy. The results reveal the existence of a distinct peak, surrounded by a thin, region of high probability. Thus, the measured temperature profile appears to be quite informative, leading to a highly localized posterior.

Consistent with the observations for the joint posterior, the marginal distributions of *D*_{0} and *E _{a}*, plotted in Fig. 10, reflect the occurrence of sharp peaks. Estimates corresponding to the MAP,

*D*

_{0}= 7.436 × 10

^{−11}m

^{2}/s and

*E*= 50.59 kJ/mol, are in good agreement with the results of the regression analysis. (Note that the values obtained from the regression analysis fall within the high probability region and are close to but do not coincide with the MAP estimates.) To verify the predictions, simulations are conducted using the MAP values of the parameters for a multilayered foil with

_{a}*λ*= 68 nm, and the resulting temperature profiles are contrasted with the corresponding measurements. As can be observed in Fig. 11, there is an excellent agreement between the simulations and experimental data. Also, note that the amplitude of the fluctuations seen in the experimental profiles compares favorably with the most likely value of

*σ*. This provides further confidence in the validity of the predictions.

_{T}Results of the Bayesian inference of the Arrhenius parameters of the high-temperature branch are reported in Figs. 12 and 13. Fig. 12 shows the marginal distributions of *D*_{0} and *E _{a}*, as well as their joint posterior. Note that the extension of the tails of the posterior distributions, beyond the prior range in both cases is an artifact of the kernel density estimation. In contrast to the results discussed earlier, in the present case, the data do not lead to a tight localization of the posterior distribution, as observed for the low-temperature branch. Specifically, the marginal distributions appear to be bimodal, with broad distributions around the local maxima. Consistent with these observations, the joint posterior distribution reflects multiple regions of high probability, indicating

*non-identifiability*of an optimal pair of diffusion parameters in this regime, due to the limitation in the amount of experimental data and the variability in the available measurements. Nonetheless, MAP values can still be clearly determined, $D0\u22483.2\xd710\u22129\u2009m2/s$ and

*E*= 56 kJ/mol. These values are remarkably close to those determined in the regression analysis. While the agreement is reassuring, it underscores potential limitations of simplified optimization approaches, which generally do not readily yield information regarding the uncertainties in the corresponding predictions.

_{a}In light of the variability of the velocity measurements, we rely on a posterior predictive check in order to verify the computations and to quantify the uncertainty in the predictive model. To this end, the MCMC samples are used to simultaneously evaluate the mean and standard deviation of the velocity predictions for each of the bilayers at which velocity measurements are collected. Results are shown in Fig. 13; plotted for comparison are the individual reaction velocity measurements. Clearly, in the present case, limitations in the data do not lead to sharp predictions. Specifically, significant uncertainty in the predictions remains, with deviations that are consistent with the variability in the underlying data.

### C. Propagation velocity

In this section, we briefly analyze details of the velocity–bilayer curve for a wider range of bilayer thicknesses than available experimentally, and contrast the results with recent experiences gained with the Ni-Al system.^{32} In particular, the latter indicated that the velocity bilayer curve for Ni-Al multilayer exhibits an anomalous velocity plateau, similar to what is observed experimentally with the Zr-Al system.^{23} The study in Ref. 32 also suggested a thermal mechanism for the occurrence of this anomaly and observed that for the Ni/Al system the plateau occurs at bilayers such that $w/\delta \u22485%$, where *w* is the premix width. Thus, we also examine whether similar trends occur using the presently developed model for equimolar Zr-Al nanolaminates.

To address these questions, computations were performed for different values of *w* and for bilayer thicknesses ranging from 10 nm to 100 nm, using MAP of the Arrhenius parameters in both temperature branches. Note that *w* can be greatly affected by deposition conditions and can be manipulated by low temperature annealing. Consequently, the range of measured premixed widths can be quite wide.^{23} We have not attempted a systematic analysis of the impact of *w* over the range of possible values, but have rather restricted our attention to two representative values, *w* = 0.5 nm and *w* = 0.8 nm. These values lie within a typical range of experimentally measured premixed widths and are sufficiently separated so as to lead to significant variations within the range of bilayers considered.

The velocity predictions are plotted against the bilayer in Fig. 14, and curves are generated for both values of *w*. Consistent with well known trends for nanostructured multilayered systems, the results indicate that following a peak at small *δ*, the velocity curve generally decays as *δ* increases. However, in the present case, and similar to experiences in Ref. 32, a velocity plateau is also observed. For *w* = 0.5 nm, the plateau can be observed at about *δ* = 10 nm, whereas for *w* = 0.8 nm it occurs near *δ* = 16 nm. Thus, the velocity anomaly in the present velocity-bilayer curves has features that are quite similar to those observed in Ref. 32 for the Ni-Al system. Specifically, the location of the anomaly shifts to higher bilayers as the premix width increases and occurs at $w/\delta \u22485%$. While the present observations support the analysis in Refs. 23 and 32, additional experimental observations are clearly needed in order to test the validity of the thermal mechanism hypothesis and characterize how the occurrence and properties of the velocity anomaly depend on the composition of the multilayered system.

## VI. CONCLUSIONS

A reaction model is developed for the equimolar Zr-Al multilayered foils, also referred to as nanolaminates. Experimental data pertaining to temperature measurements during homogeneous reactions, and reaction velocity measurements during self-propagation, are exploited in order to characterize atomic intermixing rates. The intermixing rates are described in terms of Arrhenius diffusivity correlations, which are first estimated using a regression analysis. A more elaborate Bayesian inference methodology is needed to examine the regression predictions and to assess the impact of uncertainty in measurements as well as potential discrepancies between the model predictions and experimental data. The reaction model incorporating optimal values of the atomic mixing parameters is finally used to briefly examine the behavior of the reaction front velocity for a wide range of bilayer thicknesses. Based on the results and analysis presented in this work, the following conclusions can be drawn:

Diffusivity estimates from a regression analysis exhibit two distinct branches, marked by a sharp jump in atomic intermixing rates around the melting point of aluminum. The jump is attributed to an order of magnitude increase in the diffusivity pre-exponent.

Bayesian inference analysis of the Arrhenius parameters in the low temperature branch yields MAP estimates that are consistent with the regression analysis. The posterior exhibits a tight distribution around the MAP values.

The MAP estimates of the Arrhenius parameters of the high-temperature branch are in close agreement with optimal values determined using regression. However, in this case, the Bayesian inference indicates large uncertainty levels arising due to the scatter in the velocity data and the availability of data in a narrow range of bilayers.

The reaction velocity trends for the equimolar Zr-Al system are observed to be similar to the Ni-Al system. In particular, following a peak at small bilayers, the velocity generally decreases as the bilayer thickness increases. A velocity anomaly is also observed, in the form of a distinct velocity plateau whose location depends on the premix width.

Work is currently underway to further examine the validity of the present correlations for the atomic mixing rates, refine the predictions as additional data is collected, and demonstrate the performance of the reaction model in transient multidimensional settings.

## ACKNOWLEDGMENTS

This work was supported by the Defense Threat Reduction Agency, Basic Research Award No. HDTRA1-11-1-0063, to Johns Hopkins University, and by the US Department of Energy (DOE), Office of Science, Office of Advanced Scientific Computing Research, under Award No. DE-SC0008789. The authors also acknowledge the contributions of Dr. Gregory M. Fritz towards designing the hot plate ignition tests for nanolaminate foils as well as designing and constructing the set-up for reaction velocity measurement.

## References

_{9}Ni

_{2}phase and analysis of its formation