PNNL-20738

Prepared for the U.S. Department of Energy

under Contract DE-AC05-76RL01830

Challenges in Modeling the

Degradations of Ceramic Waste
Forms

R Devanathan
F Gao
X Sun

September 2011

DISCLAIMER
This report was prepared as an account of work sponsored by an agency of the
United States Government. Neither the United States Government nor any agency
thereof, nor Battelle Memorial Institute, nor any of their employees, makes any
warranty, express or implied, or assumes any legal liability or responsibility
for the accuracy, completeness, or usefulness of any information, apparatus,
product, or process disclosed, or represents that its use would not infringe
privately owned rights. Reference herein to any specific commercial product,
process, or service by trade name, trademark, manufacturer, or otherwise does not
necessarily constitute or imply its endorsement, recommendation, or favoring by
the United States Government or any agency thereof, or Battelle Memorial
Institute. The views and opinions of authors expressed herein do not necessarily
state or reflect those of the United States Government or any agency thereof.

PACIFIC NORTHWEST NATIONAL LABORATORY

operated by
BATTELLE
for the
UNITED STATES DEPARTMENT OF ENERGY
under Contract DE-AC05-76RL01830

Printed in the United States of America

Available to DOE and DOE contractors from the
Office of Scientific and Technical Information,
P.O. Box 62, Oak Ridge, TN 37831-0062;
ph: (865) 576-8401
fax: (865) 576-5728
email: [email protected]

Available to the public from the National Technical Information Service,

U.S. Department of Commerce, 5285 Port Royal Rd., Springfield, VA 22161
ph: (800) 553-6847
fax: (703) 605-6900
email: [email protected]
online ordering: http://www.ntis.gov/ordering.htm

This document was printed on recycled paper.

(9/2003)

Waste IPSC Project Report

Challenges in Modeling the Degradation of Ceramic Waste Forms
August 2011
Ram Devanathan1, Fei Gao1, Xin Sun1
1

Pacific Northwest National Laboratory, Richland, WA 99352

Abstract
We identify the state of the art, gaps in current understanding, and key research needs in the area
of modeling the long-term degradation of ceramic waste forms for nuclear waste disposition. The
directed purpose of this report is to define a roadmap for Waste IPSC needs to extend capabilities
of waste degradation to ceramic waste forms, which overlaps with the needs of the subcontinuum
scale of FMM (fundamental methods and method) interests. The key knowledge gaps are in the
areas of i) methodology for developing reliable interatomic potentials to model the complex
atomic-level interactions in waste forms; ii) characterization of water interactions at ceramic
surfaces and interfaces; and iii) extension of atomic-level insights to the long time and distance
scales relevant to the problem of actinide and fission product immobilization.

1. Introduction
The safe long-term immobilization of high-level nuclear waste, actinides and fission products in
spent nuclear fuel, and plutonium from nuclear weapons, is a global Grand Challenge problem.
With growing interest in the expansion of nuclear energy to reduce greenhouse gas emissions
arising from energy use, the issue of nuclear waste disposal has gained renewed importance.
Toxic radioisotopes with long half-lives (~24200 years for 239Pu) have to be isolated from the
biosphere with confidence that they will not be released into the groundwater within a time scale
of the order of tens of thousands of years or longer. Short-lived fission products, such as 90Sr
and 137Cs, produce heat through radioactive decay, which can accelerate chemical alteration of
the waste form during the first 100 years of storage. In addition, highly mobile isotopes, such as
129
I and 99Tc, pose special challenges to the nuclear waste immobilization.
The international consensus on waste disposal [1] favors emplacement of waste in a
stable geological formation with multiple engineered barriers to radionuclide release, including
the waste form, waste container, drip shields, and compacted clay, augmenting natural barriers
such as unsaturated rock. The solid waste form represents the first and most essential safeguard
in any long-term waste disposal strategy. Borosilicate glass is currently the preferred waste form
for high-level waste [2]. Glass has the advantages that it is easy to process; vitrification
technology is well established; it can accommodate a variety of waste streams because of its
chemical flexibility; and it is unlikely to undergo extensive structural changes during storage
given its amorphous structure. However, glass has the disadvantages that the waste loading has
to be small; devitrification is possible; glass tends to crack easily; and glass is not chemically
durable in that it tends to degrade and release radionuclides over the long term [3]. As an
alternative to borosilicate glass, radiation-tolerant ceramics tailored for specific waste streams
and high waste loading have been the subjects of active research for nearly six decades [2-11].
Radionuclides can be distributed in a ceramic waste form in a number of different ways.
They can be accommodated on a specific lattice site, such as the Gd site in Gd2Zr2O7, in a single
phase ceramic. The stability of this arrangement depends on the ionic radius, type of bonding,
and charge balance considerations that may result in creation of lattice defects for the case of
aliovalent substitution [8]. Grains of radioactive phases can be also embedded in a ceramic
matrix [8]. Polyphase ceramics [6], such as Synroc and Supercalcine, are equivalent to a
synthetic rock made of multiple ceramic phases. For instance, Synroc-A consists of hollandite
(BaAl2Ti6O16), perovskite (CaTiO3), zirconolite (CaZrTi2O7), Ba-feldspar (BaAl2Si2O8), kalsilite
(KAlSiO4), and leucite (KAlSi2O6). Synroc-C consists of hollandite to immobilize Cs,
zirconolite to immobilize actinides, perovskite to incorporate Sr, and rutile. It can accommodate
up to 20% loading of high-level waste [11]. Synroc shows good radiation tolerance, and its
resistance to leaching has been shown to be far superior to that of glass in laboratory experiments
[6]. Other ceramic phases proposed for actinide immobilization include fluorite-structured
3

ceramics, pyrochlores, muratite, garnet-structured ferrites, zircon, apatite, sphene, pollucite,

monazite and xenotime [3, 11].
The desirable characteristics proposed for waste forms [11] include high aqueous
durability, high chemical flexibility to accommodate a variety of substitutional elements, high
waste loading to minimize the size and cost of the repository, excellent radiation tolerance over
millennial time scales, low volume swelling to avoid stress accumulation, and the availability of
natural analogs. The last criterion is important to evaluate the long-term performance of a waste
form and to validate computer models. Some durable minerals, such as zircon, occur naturally
with ~100 ppm to a few % of actinides and are known to be the oldest minerals on earth (~4400
Myr old) [12]. However, natural analogs are available only for some ceramics of interest and the
repository conditions may be different from that experienced by the mineral. The waste loading
is also likely to be different. Another major research challenge is the difficulty in extrapolating
from accelerated laboratory experiments, conducted over hours to months at orders of magnitude
higher radiation dose rates and extreme chemical environments relative to anticipated repository
conditions, and achieve a reliable understanding of the changes in physical, chemical, and
mechanical properties of the waste form in a repository over 10000 years or longer. Moreover,
the repository conditions over the life of the waste form are not well characterized and may
include extreme conditions that cannot be reproduced in a laboratory.
The limitations posed by laboratory experiments and the study of natural analogs can be
overcome by performing computer simulations to provide complementary information, model
extreme conditions and what-if scenarios, and perform sensitivity analysis. Experiments in turn
can serve to refine the model and validate the model predictions. Ongoing rapid increases in
computing power, memory, bandwidth, and storage have made such an integrated approach to
waste management feasible. However, existing modeling scheme are limited in their ability to
seamlessly describe phenomena across length scales from nanometer to meter and timescales
from nanosecond to millennium, which has limited the application of computer models to the
evaluation of waste form performance. The aim of this report is to survey the state of current
waste form modeling capabilities, identify knowledge gaps, and recommend future research
directions. The scope has been restricted to waste forms as opposed to the whole repository to
focus the effort to the principal barrier and to avoid uncertainties arising from the fact that an
ultimate repository location has not been selected in most countries.
This report reflects needs at the subcontinuum scale that overlap FMM interests, but its directed
purpose is to define a roadmap for Waste IPSC needs to extend capabilities of waste degradation
to ceramics.

2. Modeling Schemes
4

Modeling of nuclear waste forms is essential to obtain fundamental understanding of

atomic level processes and their influence on long term waste form performance under a variety
of conditions. Given the many orders of magnitude in length and time involved in this endeavor,
it is not feasible to use a single model to span this range. Hierarchical multiscale modeling
schemes that involve several models applicable to different scales of the problem are well-suited
to the task of studying waste forms as shown in Fig. 1. The key challenges in this approach are
in the process of information passing or coupling between models (represented by arrows in Fig.
1), verifying and validating the at-scale models and the coupling, estimating the error
propagation across scales, and finally quantifying the uncertainty of the overall model
predictions.
At the most fundamental level, density functional theory (DFT) [13-19] is widely used to
study the properties of point defects and impurities, such as formation, migration, and binding
energies, in materials. At the heart of this method, is the Kohn-Sham equation [13] that provides
an approximate solution for the electronic structure of materials through the exchangecorrelation functional. Many different functionals are available to the DFT practitioner. It is
important to understand the strengths and limitations of these functionals, and reasons for
discrepancies between results reported in the literature [18, 20]. The strong localization of felectrons in materials that contain actinides is a challenge for DFT, because widely used
functionals cannot treat localized electrons correctly. The use of hybrid functionals improves
accuracy but increases the computational intensity of DFT. The current state of the art for
actinide-bearing materials is the use of an additional mean-field Hubbard parameter, U, to the
Hamiltonian [21-23]. Schreckenbach and Shamov [24] have recently discussed the various
approximations used in computational chemistry of actinides. Due to the limitations of current
functionals, DFT results will quantitatively differ from experimental results [25].
Despite the above mentioned limitations, DFT calculations can serve as the foundation,
i.e. the most accurate and reliable level possible, of the multiscale modeling scheme for
crystalline ceramics. DFT-based molecular dynamics simulation [26] is increasingly being used
to investigate the dynamical behavior of complex chemical systems. It is highly relevant to the
study of nuclear waste forms, especially the interactions of water with model ceramic surfaces.
Unfortunately, the reliability of DFT comes with high computational cost, which limits the
application of this method to systems containing fewer than 1000 atoms for time scales of 1
nanosecond or less. This limitation makes it difficult to model realistic surfaces and interfaces,
which hold the key to ceramic dissolution processes. Gale and Wright [25] have pointed out that
despite the ongoing increases in computational power many problems of interest continue to lie
beyond the scope of DFT. Therefore, the main value of DFT in the study of ceramic waste forms
lies in its ability to provide reliable data to fill gaps where experimental data is not available and
to parameterize the next higher scale (coarser model) in the modeling hierarchy. High
5

throughput DFT calculations [27] can be used to rapidly screen waste form candidates and
optimize waste form composition.
At the next higher level, simulations of systems containing a million to a trillion atoms
are possible with classical molecular dynamics (MD) and Monte Carlo (MC) simulations,
although the current state of the art for ceramics is about 1 billion atoms. In simulations of such
scale the electronic degrees of freedom are ignored and complex interactions between atoms are
simplified in the form of parameterized classical potentials also known as force fields. The
potential plays a crucial role in the reliability of the simulations, and its proper determination,
especially in the case of complex ceramic systems with multiple cation and anion sublattices, is
essential to the success of MD. While important advances have been made towards obtaining
accurate classical potentials for simple metallic systems [28], their construction, especially for
non-metals and complex configurations, still remains a challenging and time consuming process.
Much of the difficulty in generating reliable potentials for ceramics is caused by the lack
of well-defined procedure to obtain the optimal interaction parameters. Traditionally a limited set
of experimental data from only one or, at best, several reference configurations is used in the
fitting process. The resulting potential is not transferable and its accuracy in real dynamical
simulations is questionable. This is especially so if the potential is parameterized with available
equilibrium property data, but is subsequently used to model conditions far from equilibrium.
Due to the time consuming nature of the fitting process, the parameters of interaction, once
determined, remain fixed for the duration of the simulation. This restricts the applicability of
traditional MD simulations to important class of problems that involve significant dynamical
changes in the chemical environment, such as water-mineral surface interactions, species
transport in an aqueous medium, defect clustering and precipitation, and interaction of adsorbates
with nanoparticles.
Recently, van Duin et al [29, 30] have developed reactive force fields that permit bond
breaking to model complex chemical environments without the computational cost of DFT. The
important aspects of this force field are that each atom has a Gaussian charge distribution that
can transfer to other atoms in response to changes in the local environment, the parameters are
obtained from DFT calculations, and there are no predetermined reactive sites [30]. Reactive
force fields are available for barium titanate [30], aluminum oxide [31], and silica [32].
Similarly, charge-optimized many-body potentials have been developed to model Si/SiO2 and
Hf/HfO2 interfaces [33, 34]. However, there remains the need to extend this methodology to the
heavy-element compounds that are being considered as candidate ceramic waste forms.
Several non-reactive potentials have been used to model mineral-water interactions.
CLAYFF is a general force field developed by Cygan et al [35] to describe water, hydroxyl ions,
and surface species. It can be used to model clay, hydroxide and oxyhydroxide phases and their
interfaces with aqueous solutions [35]. Rustad [36] has reviewed molecular models of the
6

oxide-water interface. Recently, a dissociative water model was developed and used to model the
chemisorption of water molecules at a silica surface [37]. Hughes et al [38] have used quantum
mechanical calculations to fit partial charge shell models and rigid ion models to study water
interactions with zeolites. The shell model [39] represents the polarizability of ions by modeling
a given ion as a massive core attached by a spring to a massless shell. The shell model provides
a good representation of the equilibrium phase of ceramics [40-42]. However, shell models are
more computationally intensive than the corresponding rigid ion models.
A key research need to address the limitations of current interatomic potentials is the
development of a new framework that offers a concise, systematic, and efficient procedure for
generating parameters of the interatomic interaction. A promising approach to achieve this goal
is the Force Matching Method (FMM) of Ercolessi and Adams [43]. The FMM uses a large
amount of data obtained from ab initio calculations of various configurations and Car-Parinello
MD simulations [26] to fit the interatomic potential. Atomic forces are determined from first
principles for a large number of configurations including geometries such as clusters, defects,
surfaces, and liquids. The functional form and parameters of the potential are optimized to
match the first principles forces by performing a minimization using more than a hundred
parameters. The method is not limited to presently-used analytical forms, which is an advantage
for the simulation of complex materials for which potentials do not exist at present. FMM has
been successfully applied to the development of potentials for metals [44] and Si [45].
The starting point for this potential development effort is a database of configurations
obtained from experimental data, and DFT calculations of perfect crystals, defect configurations
and the amorphous state. Fig. 2 shows a flow chart of this approach. The fitting of the potential
begins with an initial guess for the form of the interaction. The input configurations will be used
to determine forces and energies using DFT. In addition to these forces, experimentally
determined lattice and elastic constants of various polytypes can be included in the fit. The data
elements of the fit can be assigned various relative weights. Several strategies can be used to
speed up the process at this point. The calculations for different configurations can be delegated
to an independent group of processors giving rise to a very efficient parallel calculation. Based
on ab initio forces, the classical potential can then be refined using a least-squares fit. The
interatomic potential generator can be used as a stand-alone module, and the refined potential
can be used to perform nanosecond to microsecond scale classical MD simulations of radiation
effects in ceramics and ceramic-water interactions.
In addition, the potential-generation algorithm can be dynamically coupled to a classical
MD simulation as shown in Fig. 3. By performing two concurrent simulations a classical MD
simulation of a large ceramic system and a smaller DFT calculation for a small region of interest,
such as a crack tip, one can generate potentials for complex environments on the fly [46]. The
larger MD simulation will define the region of interest and provide its coordinates for the DFTbased simulation, which in turn will refine the parameters of the interaction for the region of
7

interest. By loosely coupling the two simulations, there is the potential to implement the method
efficiently on thousands of processors. Such a scheme would be invaluable to model charge
transfer processes in ceramics, especially at interfaces and near defect clusters or second phases.
Adaptive quantum mechanics/molecular mechanics (QM/MM) schemes [47] in which water
molecules can enter or leave the quantum mechanical region of interest can be quite valuable in
studies of mineral-water interactions. However, major issues remain in the implementation of
embedding methods, such as the need to rigorously define the boundary between the QM and
MM regions and match forces at the interface [48].
At the mesoscale level, the kinetic Monte Carlo (KMC) method [49] and the phase field
(PF) model [50] are very useful to extend the spatial and time regimes of the MD simulations to
deal with defect accumulation and gas bubble nucleation, phase stability and phase
transformation over years. This is essential for upscaling atomic-level calculations to the
macroscopic regime of experimental measurements. The KMC and PF simulations depend on the
inputs such as defect migration energies, entropies, binding energies, and thermodynamics of
different phases, as well as the initial micro-structural features of the material. In KMC, we
consider a set of random and independent events Wi, that occur with average rates ri as a Poisson
process. The sequence of events and the time between events can be constructed from a
probability distribution that is weighted based on the rates of the events. As reactions and
transitions occur during the simulation, the list of possible events and their rates have to be
recalculated to trace the evolution of the system. The scale of the process adapts to the new ri
and the system can be followed over much longer time scales relative to classical MD. The
major stumbling block with KMC is that a list of all possible events and their rates must be
known to carry out the simulation. These inputs can be determined from experimental data
bases, or from data generated by first-principles calculations or classical MD simulations, but it
can become an onerous task for systems with chemically and geometrically complex interfaces.
In contrast to the stochastic KMC method, the phase field model is a deterministic
thermodynamic model that is used to describe microstructural evolution in terms of continuum
phase field. The PF can be a descriptor of microstructure or composition. The field equations
are solved on a two- or three-dimensional grid with grid spacing larger than the atomic scale but
smaller than the grain scale. As in the case of the KMC method, the PF model requires
considerable thermodynamic and kinetic input parameters from experiment or finer length scale
simulations, such as defect formation, migration, binding and interaction energies. Without such
information, reliable quantitative modeling is not possible. The results of mesoscale models can
be used to develop constitutive equation describing material behavior.
3. Modeling of ceramic waste forms
Modeling and simulation studies of waste forms have focused mainly on atomistic details
of radiation damage. Sickafus et al [51] predicted the radiation tolerance of A2B2O7 pyrochlores
8

using molecular statics calculations with fixed charge empirical potentials. Radiation tolerance
was attributed to the ability of the material to tolerate lattice defects produced. Subsequent
molecular statics studies of A4B3O12 compounds [52] by the same group showed that radiation
tolerance can be inferred from the phase diagram and that compounds that accommodate lattice
disorder tend to resist amorphization. MD simulations of A2B2O7 pyrochlore by Chartier et al
[53, 54] have shown that the accumulation of dumbbell interstitials on the B cation sublattice
results in amorphization and that the annealing of the dumbbells is kinetically not feasible at low
temperatures. Recently, Devanathan et al [55] have used molecular statics calculations and MD
simulations to examine the effects of defect accumulation on the structural and mechanical
stability of Gd2Zr2O7 and Gd2Ti2O7 pyrochlores. The results reveal that seemingly minor
differences in radiation damage accumulation and annealing between the two pyrochlores
eventually lead to divergent responses to radiation. The cation Frenkel pair accumulation was
found to drive Gd2Ti2O7 amorphous with a decrease in density of about 13% and decrease in
elastic modulus of 50% while amorphization did not occur in Gd2Zr2O7 due to more effective
dynamic annealing as shown in Fig. 4. Such effective dynamic annealing has also been observed
in MD simulations of energetic recoil damage in yttria-stabilized zirconia (YSZ) [56], which has
the related fluorite structure. Since Gd2Zr2O7 and YSZ are fast ion conductors [57], rapid anion
vacancy migration helps to anneal damage and to prevent radiation damage accumulation from
making the disordered crystal thermodynamically unstable with respect to the amorphous
material. The combination of a mechanism to accommodate radiation damage with minimal
increase in energy and volume and fast anion transport to enable dynamic annealing appears to
be the key to radiation tolerance in compounds based on the fluorite structure [52, 55, 58]. High
values of the displacement threshold energy, i.e. the minimum energy to be imparted to a lattice
atom to produce a stable defect, seem to be good indicators of radiation tolerance of ceramics
[59, 60].
Trachenko et al [61-66] have examined the amorphization resistance of a number of
ceramics using atomistic simulation. They observed that the extent of radiation-induced
structural damage has an effect on activation energy barriers for damage recovery and that the
dynamic recovery processes at the picosecond scale have a strong influence on resistance to
amorphization. Devanathan et al [67-69] performed MD simulations of energetic recoil damage
in zircon and observed direct-impact amorphization. Unlike the case of zirconia and Gd2Zr2O7
pyrochlore in which isolated defects are produced by radiation damage, energetic recoils in
zircon produce an amorphous core surrounded by defects as shown in Fig. 5. At the nanoscale,
radiation damage induces phase separation into zirconia-rich and silica-rich regions [68]. Yu et
al [70] found a 11% volume expansion and 60% decrease in bulk modulus as a result of
amorphization of crystalline zircon. Large volume swelling can induce cracking of wasteforms
and enable the leaching of actinides during long-term storage. Trachenko et al [65-66] have
connected the amorphization resistance to the nature of the chemical bond by stating that a
complex material can be amorphized by irradiation if it can form a covalent network. The
9

authors have indicated that other criteria proposed in the literature such as those based on
coordination number, radius ratio and topology may ultimately be tied to the nature of the
chemical bond.
In addition to DFT calculations and MD simulations, the kMC method has been used to
study the stochastic production of defects in the crystal lattice, recombination of defects, and the
identification of amorphous regions in ZrSiO4 [71]. Within this model, amorphous regions were
identified as those having a critical density of Zr vacancies. The simulated interstitial content and
amorphous fraction as functions of dose were found to be consistent with experimental data at
300K for 238Pu-doped zircon. While this simple model provides reasonable agreement with
experimental data, it does not yet fully describe all the atomic- level processes associated with
cascade overlap and relaxation. Further improvements in the model will require additional
experimental data and realistic atomic-scale information from DFT calculations of defect
reaction kinetics and defect relaxation. More recently, a hybrid kinetic lattice Monte Carlo
(KLMC) model has been integrated with MD simulation and applied to investigate the recovery
and clustering of defects during annealing of a single 10 keV cascade in cubic silicon carbide
[72]. A method of directly transferring the defects created by MD simulations to the KLMC
model was developed, while the KLMC model parameters were obtained from MD simulations
and ab initio calculations of defect migration, recombination, and annihilation. The transfer of
defect configurations between MD and Monte Carlo is important for multiple cascade defect
accumulation simulations since the atomic relaxation during defect accumulation has to be
properly considered. This overcomes an important weakness of previous KMC simulations of
defect accumulation.
In addition to studies of defect accumulation and recoil damage in waste forms, there is
growing interest in studying the response of ceramics to swift heavy ion (SHI) damage.
Simulations of the electronic stopping process, for instance using a thermal spike model, can
provide insights into the competition between amorphization and recrystallization in irradiated
ceramics. Pakarinen et al [73] have studied SHI damage in quartz and amorphous silica. Their
results shed light on the density variations in the radiation damaged region. Areas of lower
density may provide pathways for actinide and other species diffusion. Crocombette [74] has
compared defect production by displacement cascades and thermal spikes and concluded that
while the number of defects produced in simulations of cascades and spikes can be different, the
overall characteristics of the damage are quite similar. Recently, Zhang et al [75] showed that
thermal spike simulations can reproduce the microstructural features of ion tracks and the
amorphization tendency observed in experimental irradiation of pyrochlore compounds. These
findings suggest that thermal spike simulations can be used to quickly screen candidate waste
forms based on their tendency to undergo direct-impact amorphization.
While most simulations of radiation damage to date have focused on single crystals, there
is growing interest in simulating defects and diffusion in polycrystalline materials. The kMC
method has been used to explore the grain boundary diffusion in ceramics [76], which controls
10

many solid-state processes such as sintering, creep, and the growth of oxide films on metals or
ceramics. The basic problem with simulations of grain boundaries in ceramics is that the
activation energies for ion motion are so high (of the order of 1 eV) that it is not feasible to
extract diffusion coefficients from MD simulations except at unrealistic temperatures. The
problem with using static simulations, on the other hand, is that boundaries often have a large
number of different possible hops, all of low symmetry, so that it is difficult to find the saddle
points and obtain the individual hopping rates. While the calculated activation energies for grainboundary diffusion are in reasonable agreement with experiment, problems in comparing the
simulations with experiments performed on polycrystals are still issues to be addressed with
further advances in simulation methods.
In contrast to the extensive atomic-level simulations of radiation damage in waste forms,
there have been limited modeling studies of water interactions with waste form surfaces. Bates
et al [77] have stated that existing models of actinide mobility may underestimate the possibility
of radionuclide release. Nangia and Garrison [78] have reviewed the state of theory in the study
of dissolution at mineral-water interfaces. Most of the studies performed to date have used atscale models and there is a need to span the gap from ab initio studies of small clusters to
continuum models of dissolution. Most of the studies have focused on silicate minerals. Ab
initio molecular orbital calculations by Xiao and Lasaga [79] using disilicic acid, (HO)3Si-OSi(OH))3, to represent the quartz surface have shed light on the rate-determining step in the
dissolution of SiO2 by hydrolysis. This is at odds with the subsequent findings of Criscenti et al
[80] based on a similar study using larger clusters and explicit water molecules, which
demonstrates the need for simulating realistic waste form surfaces. However, modeling large
surface structures using DFT is computationally expensive. At the next level, MD simulations
have been used to study the reactivity of quartz surfaces using models of varying complexity,
such as shell models and dissociative water models [81, 82]. Given the limited ability of MD
simulations to go beyond a few hundred nanoseconds, kinetic models have been developed to
describe dissolution at feldspar surfaces [83]. Recently a time-independent reactive Monte
Carlo-configurational bias Monte Carlo hybrid algorithm has been developed to study
dissolution at silicate-water interfaces [84]. As a part of the NEAMS Waste IPSC, Argonne
National Laboratory has been used ab initio (DFT) and molecular dynamics methods to
calculate salvation and pH values in glass waste forms and to study the dissolution kinetics of
orthoclase feldspar surfaces. These investigations have direct relevance to the ceramics waste
forms at surfaces.
Stochastic methods can be used to improve our understanding of the dissolution process
by testing basic assumptions of the analytical models [85]. In glasses, the simulated results
predict the experimentally observed silica condensation that results in a protective effect of the
alteration layer and an initial silica dissolution rate decreasing as a function of time [85]. The
simulated results also demonstrate that the silica solubility in solution is not that of the highest
energy state (the glass), but depend on the exact parameter values, which is consistent with the
fact that most silica bearing sites at the surface with the solution are in that lowest energy state.
11

This result confirmed that dissolution kinetics leads to different steady state concentrations than
thermodynamic equilibrium, which would suggest that, if two energy sites exist, the observed
solubility in solution corresponds to the site with the highest energy state. Instead, the observed
solubility lies between the solubility corresponding to both sites and its value depends on the
relative amounts of both sites. The silica solubility also depends on the amount of silica particles
in contact with the solution (the surface), implying that replacing part of the silica by another
poorly soluble glass network former lowers the silica solubility. However, the representation of
the surfaces and materials in terms of composition in the Monte Carlo simulations are generally
very simple, and the Monte Carlo predictions are qualitative, but not quantitative. Quantitative
predictions require more atomic level or thermodynamic parameters (e.g. a better description of
the interaction between particles, greater chemical complexity in the model, and more complex
surfaces). This will allow realistic simulations of wasteform dissolution.
Finally, it is worth noting that the phase field model has emerged as a powerful
computational approach at the mesoscale for predicting phase stability and microstructural
evolution under irradiation [86]. Phase field models have been developed to predict gas bubble
microstructure evolution in nuclear fuels [86], void migration and growth kinetics in materials
under irradiation and temperature field [87] shown in Fig. 6, void evolution and swelling in
materials under irradiation [88], and void nucleation and growth in irradiated metals [89]. In
these models, only the diffusion of single vacancies and interstitials is considered. Mobile defect
clusters, such as vacancy clusters or small interstitial loops, are viewed as an assembly of
individual vacancies or interstitials in this approach. Thus, they contribute to the overall defect
concentrations and to the effective mobilities of defects. However, larger and immobile defect
clusters are viewed as sinks or nucleation sites for defect clusters, which may affect the net
generation rate of point defects. A general assumption is that vacancies and interstitials generated
by fission fragments and neutron damage cascades can be absorbed by structural defects, such as
dislocations and grain boundaries, and thus, the net increase of vacancies and interstitials
depends on the generation rates of fission fragments and neutrons. This also depends on the
defect generation rates, and sink strengths, which are determined by the type and density of point
defect sinks. On the other hand, the driving force of vacancy and interstitial diffusion includes
the chemical potential gradients and the concentration gradients, correlated with the
microstructure evolution driven by the minimization of the total free energy of the system.
Moreover, the microstructure evolution path depends on thermal fluctuations, defect generation,
recombination, and defect sink interactions in irradiated materials. The simulated results for
irradiated materials demonstrate the ability of phase-field approach to model the nonequilibrium
process of microstructure evolution.
4. Research needs

12

The development of science-based models to predict the long-term behavior of ceramics

waste forms for nuclear waste disposition represents an important opportunity that requires a
fundamental understanding of atomic-level processes (such as defect generation and migration,
incorporation of fission products, and their interactions with microstructures) and long term
performance under extreme conditions (such as phase stability, phase transformation, radiation
tolerance, and leaching under different temperature and pH conditions). Although significant
progress has been made to model the complex behavior of ceramic nuclear waste forms, ranging
from atomic-level processes to meso-scale mechanisms of radionuclide incorporation, chemical
corrosion, and alteration mechanisms, and the response to different types of radiation fields, the
present understanding of their performance over long time periods is insufficient in detail, even
for the case of atomic-level simulations as outlined above. The ability to develop integrated
models rests not only on our fundamental understanding of the critical processes but also on the
computational ability to link atomic-level understanding to macroscopic phenomena. The
challenge in making predictions about the long-term performance of a waste forms is to identify
the key mechanisms governing fundamental processes in ceramics waste forms. Key research
needs for ceramic nuclear waste forms are recommended below, with the primary objective to
develop fundamental knowledge and models at the atomic, microscopic, and macroscopic levels
in order to enable performance assessments under a variety of conditions.
(i) Advances in the basic science of f electron systems in materials may offer the opportunity to
extend condensed matter physics and reaction chemistry to accurately predict the performance of
ceramic nuclear waste forms. This is a key enabling science for the fundamental understanding
of actinide-bearing materials and solutions.Two major physical effects that challenge the
theoretical modeling of f electron systems are relativistic effects and dynamic correlation.
Computational actinide chemistry has its own challenges, arising from the combined effects of
relativity and electron correlation, along with the size of typical systems, the large number of
electrons, and the close energetic proximity of several electronic states. As described above, a
key research need is to develop a robust theory of the electronic structure or exchange and
correlation functionals in density functional theory for actinide-bearing materials that will
provide an improved understanding of their physical and chemical properties and behavior.
Advances are needed in the application of these new electronic structure methods for f-elementcontaining molecules and solids to calculate the defect properties, incorporation of fission
products and their interaction with microstructures in multi-component systems. Also, these
improved developments will advance the fundamental understanding of related chemical and
physical properties at high temperature.
(ii) Many ceramic forms with different crystal structures, such pyrochlore, zirconia and apatite,
have been proposed as hosts for fission products (e.g. Cs, Sr and I), but they will be subjected to
high ionizing radiation fields for a long period (up to several hundred years), where charge
transfer, electronic excitations and electronic energy loss processes can enhance defect
13

migration, promote defect recovery and influence microstructural evolution in ceramics [90].
Recent studies have used ab initio MD methods (AIMD) to study low energy recoil events in
ceramics [91] and semiconductors [92, 93]. Large-scale ab initio MD methods have also been
used to study ion-solid interactions in SiC [93], and the results demonstrate that significant
charge transfer can occur between atoms, and defects can enhance charge transfer to the
surrounding atoms. Moreover, the charge transfer to and from recoiling atoms can alter the
energy barriers and dynamics for stable defect formation. These simulations have been limited to
systems containing fewer than 1000 atoms and to time scales less than 1 nanosecond.
Alternatively, Sutton et al. [94] have developed a large-scale tight binding molecular dynamics
(TBMD) method capable of dealing with electronic and atomic dynamics for up to 100,000
atoms. This method has been applied to the physics of electron migration in nanowires, currentinduced heating, and electronic excitations and their influence on interatomic forces in radiation
damage in metals [94]. However, the current model considers only the s-band with two sorbitals. While this single s-band model has been used to model radiation damage in metals,
further developments including higher band models are needed to extend this method to
ceramics, which will be crucial for the correct description of directional bonding in materials.
(iii) The use of empirical or semiempirical interatomic potentials makes it possible to simulate
large systems (~ up to a few billion atoms) for long times (microsecond), and thus to tackle such
problems as defect generation, plastic deformation, ion-solid interaction, or atomic diffusion.
However, the reliability of interatomic potentials is a key issue for the success of predictive
models for ceramic waste forms. Due to the extensive use of MD methods in materials science, a
variety of techniques have been utilized over the years to develop reliable atomic potentials. For
ceramic waste forms, long-range Coulombic interaction with partial charge models along with
short range interactions of the Buckingham form are generally employed to simulate radiation
damage and defect accumulation. As discussed above, significant charge transfer can occur
during irradiation, and the defects created can enhance charge transfer to the surrounding atoms.
Furthermore, a major weakness of fixed-charge potentials is that an interstitial dumbbell made
entirely of cations or anions is unlikely to form due to the strong Coulomb repulsion between
like ions, which is in contrast to the findings of ab initio calculations. One possible solution is the
development of charge transfer potentials [95], incorporating variable charge concepts into an
analytical embedded-ion method (EIM) potential, where the equilibrium charges are embedded
in the potential so that the energy minimization step is eliminated. As a result, the method is both
energy conserving and computationally efficient. Using the La-Br system as an example, it has
been demonstrated that this model captures key properties of the LaBr3 compound including the
crystal structure, lattice geometry, and the {1120} cleavage. The model can already provide a
viable interatomic potential for nonstoichiometric ionic systems.
A generalized framework for interatomic potential design has been recently established
[96], as shown in Fig. 7, which allows a researcher to tailor an interatomic potential towards
specific properties. This methodology produces an interatomic potential design map, which
14

contains multiple interatomic potentials and is capable of exploring different nanoscale

phenomena observed in experiments. This methodology is efficient and provides the means to
assess uncertainties in nanostructure properties due to the interatomic potential fitting process.
New advances in the development of interatomic potentials using a large data base of force fields
and atomic configurations, and in variable charge models, are needed to overcome the challenges
of adequately treating defect properties, change transfer and radiation effects. For each potential,
the limits of transferability must be established, and the potential must be validated for far-fromequilibrium conditions, such as unusual coordination, high temperature, and high pressure.
(iv) Waste form evolution over long-time scale (up to millennia) can involve extremely slow
kinetic processes operating over geological time scales, such as leaching, internal radiation
damage leading to amorphization or crystallization. The stability and safety of waste forms will
depend on our knowledge of the details of interfacial chemistry and transport for the specific
waste form and repository. New advances in mesoscale and macroscale modeling of ceramic
waste forms need to focus on the development of methods and models to describe the
fundamentals of surface chemical interactions to allow efficient modeling of multiple phase
interfaces. While significant advances have been made in kMC methods and phase field models,
they have been applied mainly to predict microstructural evolution in relatively simple systems,
such as pure metals and simple alloys. Simulating defect processes and microstructure evolution
of ceramic waste forms using the kMC method and phase field model is challenging in part
because of the inherent chemical and structural complexity.
In addition, current kMC and phase field models of materials are often qualitative, but not
quantitative. The development of quantitative kMC and phase-field models for ceramics is even
more challenging because the systems involve complex reactions consisting of numerous steps,
resulting in the formation of reactive or metastable intermediate species. Moreover, quantitative
kMC and phase-field models require accurate thermodynamic and kinetic properties of the
corresponding ceramics, which may represent long-term grand challenges of atomic-level
simulations. Lacking thermodynamic and kinetic properties, the kMC and phase-field
simulations may use reasonable thermodynamic and kinetic properties which could be obtained
from experiments and atomistic simulations. These mesoscale simulations can be used to build a
thermodynamic database of the effects of thermodynamic and kinetic parameters on the
microstructural evolution of ceramic waste forms. Such a database could play an important role
in uncertainty evaluation of both mesoscale and macro-scale simulations. At present, kMC and
phase field approaches look promising for bridging the atomistic and continuum scales,
particularly because they can be used to reproduce both near- and far-from equilibrium
dissolution and to deal with water-mineral surface interactions, species transport in an aqueous
medium, defect clustering and precipitation, and interaction of adsorbates with nanoparticles.
Finally, there is a need to develop rigorous upscaling techniques to create accurate constitutive
equations for continuum models, which can be used to evaluate the long-term performance of
ceramic waste forms under a variety of environmental conditions. Continuum models are needed
15

to account for microstructural evolution and property changes with radiation damage,
thermochemistry and thermodynamics.
(v) Important issues that should be addressed in the framework of a multiscale and multiphysics
model of ceramic waste forms include cross-verification of data obtained from different
simulation codes, validation of simulation results using experimental data and finer scale models,
and quantitative assessment of uncertainties at each scale and for the handshake between scales.
The development of a standard verification, validation and uncertainty quantification approach is
required for the multiscale scheme, and a similar approach needs to be developed for
determining the reliability of data obtained from experiments. Data used and generated by subcontinuum and mesoscale models with different codes, methods, and algorithms need to be
verified and validated, and properly qualified. In addition to developing methods for upscaling
from finer to coarser scales, quantitative estimation of the propagation of uncertainties through
the hierarchy of scales will be essential.
5. Summary
We have reviewed the state of the art and important research needs in modeling the
degradation of ceramic waste forms. Despite decades of research, modeling studies have
focused on simple crystal structures and simple chemistries, single crystals, and studies of
radiation damage at the nanoscale. There is a need to expand these studies to complex
chemistries and structures, including polycrystalline and polyphase materials. Given the reality
that density functional theory will continue to be limited to short time and distance scales, there
is a need to systematically coarse grain simulations of ceramic waste forms with rigorous
verification, validation and uncertainty quantification. Reactive interatomic potentials have to be
developed to model the complex interactions at ceramic surfaces and interfaces. Simulation of
realistic surfaces, efficient coupling from atomistic simulations to stochastic models, and
development of constitutive equations that can be used in continuum models will be crucial to
advance this field.
ACKNOWLEDGMENTS
This report was supported by the US Department of Energys Nuclear Energy Advance
Modeling and Simulation (NEAMS) Program, within the Waste IPSC, in Pacific Northwest
National Laboratory (PNNL), which is operated by Battelle Memorial Institute for the US
Department of Energy under Contract No. DE-AC05-76RL01830.

16


Figure 1. At-scale models in a multiscale modeling scheme for waste forms.

17

Make initial guess of the function

Processor

Groups
1

Determine classical
forces and energies.

Determine ab initio
forces and energies.

RefinePotential

Good
Fit?

NO

YES

Figure 2. Flow-chart for the generation of transferable potentials

Classical
MD
i l ti

In: Geometry
of the region

Ab-Initio
Calculation

Out: Classical
Interaction Parameters

Figure 3. Schematic diagram of an approach to couple ab initio calculations with classical

molecular dynamics.
18

Figure 4. Atom projections from molecular dynamics simulation of defect accumulation in

pyrochlore [55]. a) Gd2Ti2O7 is amorphous after cation Frenkel pairs were introduced in the
lattice to a level of 0.63 displacements per atom (dpa); b) Gd2Zr2O7 is crystalline after defect
accumulation to 1 dpa. The green, grey and red spheres represent Gd, Ti or Zr and O,
respectively.

19

Figure 5. Perspective projection of defects created by a 30 keV U recoil in ZrSiO4. Zr, Si and O
defects are shown as blue, light brown and red spheres, respectively. The size of the box shown
is about 13 x 5 x 4 nm. The simulation cell is much larger.

20

Figure 6. Phase field modeling of void migration and growth with time in nuclear fuels [86].

21

Figure 7. Schematic of the generalized framework for interatomic potential development

showing the Global and Local Approach stages along with their individual steps [96].

22

References
[1]
[2]
[3]
[4]
[5]
[6]
[7]
[8]

[9]
[10]
[11]
[12]
[13]
[14]
[15]
[16]
[17]
[18]
[19]

Sailor W C, Bodansky D, Braun C, Fetter S and van der Zwaan B 2000 A nuclear solution
to climate change? Science 288(5469) 1177.
Ewing R C, Weber W J and Clinard F W 1995 Radiation effects in nuclear waste forms for
high-level radioactive waste Progress in Nuclear Energy 29(2) 63.
Weber W J, Navrotksy A, Stefanovsky S, Vance E R and Vernaz E 2009 Materials science
of high-level nuclear waste immobilization MRS Bulletin 34(1) 46.
Hatch L P 1953 Ultimate disposal of radioactive wastes Amer. Sci. 41 410.
Roy R 1979 Science underlying radioactive waste management: status and needs Scientific
Basis for Nuclear Waste Management 1 1, McCarthy G J ed., (New York: Plenum
Press).
Ringwood A E, Kesson S E, Ware N G, Hibberson W and Major A 1979 Immobilization of
high level nuclear reactor wastes in SYNROC Nature 278 219.
Ewing R C and Lutze W 1991 High-level nuclear waste immobilization with ceramics
Ceramics International 17 287.
Weber W J, Ewing R C, Catlow C R A, Diaz de la Rubia T, Hobbs L W, Kinoshita C,
Matzke Hj, Motta A T, Nastasi M, Salje E K H, Vance E R and Zinkle S J 1998
Radiation effects in crystalline ceramics for the immobilization of high-level nuclear
waste and plutonium J. Mater. Res. 13 1434.
Ewing R C 1999 Nuclear waste forms for actinides Proc. Nat. Acad. Sci. 96(7) 3432.
Ewing R C, Weber W J and Lian J 2004 Nuclear waste disposalpyrochlore (A2B2O7):
Nuclear waste form for the immobilization of plutonium and minor actinides J. Appl.
Phys. 95 5949.
Lumpkin G R 2006 Ceramic waste forms for actinides Elements 2 365.
Wilde S A, Valley J W, Peck W H and Graham C M 2001 Evidence from detrital zircons
for the existence of continental crust and oceans on the earth 4.4 Gyr ago Nature 409
175.
Kohn W and Sham L J 1964 Self-consistent equations including exchange and correlation
effects Phys. Rev. 140 A1133
Parr R G and Yang W 1989 Density-Functional Theory of Atoms and Molecules (New
York: Oxford University Press)
Dreizler R M and Gross E K U 1990 Density Functional Theory (Berlin: Springer)
KohnW1999 Nobel lecture: electronic structure of matterwave functions and density
functionals Rev. Mod. Phys. 71 1253
Martin R M 2004 Electronic Structure: Basic Theory and Practical Methods (Cambridge:
Cambridge University Press)
Mattsson A E, Schultz P A, Desjarlais M P, Mattsson T R, and Leung K 2005 Designing
meaningful density functional theory calculations in materials sciencea primer
Modelling Simul. Mater. Sci. Eng. 13 R1.
Hafner J, Wolverton C and Ceder G 2006 Toward computational materials design: The
impact of density functional theory on materials research MRS Bulletin 31 659.

23

[20] Devanathan R, Van Brutzel L, Chartier A, Gueneau C, Mattsson A E, Tikare V, Bartel T,

Besmann T, Stan M and Van Uffelen P 2010 Modeling and simulation of nuclear fuel
materials Energy Environ. Sci. 3 1406.
[21] Anisimov V I, Zaanen J and Andersen O K 1991 Band theory and Mott insulators:
Hubbard U instead of Stoner I Phys. Rev. B 44 943.
[22] Dorado B, Freyss M and Martin G 2009 GGA plus U study of the incorporation of iodine
in uranium dioxide European Phys. Journal B 69 203.
[23] Devey A J 2011 First principles calculation of the elastic constants and phonon modes of
UO2 using GGA plus U with orbital occupancy control J. Nucl. Mater 412 301.
[24] Schreckenbach G and Shamov G 2010 Theoretical Actinide Molecular Science Accounts of
Chemical Research 43 19.
[25] Gale J D and Wright K 2010 Lattice dynamics from force fields as a technique for mineral
physics Rev. Mineral. Geochem. 71 391.
[26] Marx D and Hutter J 2009 Ab Initio Molecular Dynamics: Basic Theory and Advanced
Methods (New York: Cambridge University Press)
[27] Jain A, Hautier G, Moore C J, Ong S P, Fischer C C, Mueller T, Persson K A and Ceder G
2011 A high-throughput infrastructure for density functional theory calculations
Comput. Mater. Sci. 50 2295.
[28] Daw MS and Baskes M I 1984 Embedded atom method: Derivation and application to
impurities, surfaces, and other defects in metals Phys. Rev. B 29 6443.
[29] van Duin A C T, Dasgupta S, Lorant F and Goddard W A III 2001 ReaxFF: A Reactive
Force Field for Hydrocarbons J. Phys. Chem. A 105 9396.
[30] Goddard W III, Merinov B, van Duin A, Jacob T, Blanco M, Molinero V, Jang S S and
Jang Y H 2006 Multi-paradigm multi-scale simulations for fuel cell catalysts and
membranes Molec. Simul. 32 251.
[31] Zhang Q, Cagin T, van Duin A C T, Goddard W A, Qi Y and Hector L G 2004 Adhesion
and nonwettingwetting transition in the Al/alphaAl2O3 interface Phys. Rev. B 69
045423.
[32] van Duin A C T, Strachan A, Stewman S, Zhang Q and Goddard W A III 2003 ReaxFF
SiO reactive force field for silicon and silicon oxide systems J. Phys. Chem. A 107
3803.
[33] Yu J, Sinnott S B and Phillpot S R 2007 Charge optimized many-body potential for the
Si/SiO2 system Phys. Rev. B 75 085311.
[34] Shan T-R, Devine B D, Hawkins J M, Asthagiri A, Phillpot S R and Sinnott S B 2010
Second-generation charge-optimized many-body potential for Si/SiO2 and amorphous
silica Phys. Rev. B 82 235302.
[35] Rustad J R 2011Molecular models of surface relaxation, hydroxylation, and surface
charging at oxide-water interfaces Rev. Mineral. Geochem. 42 169.
[36] Mahadevan T S and Garofalini S H 2008 Dissociative chemisorption of water onto silica
surfaces and formation of hydronium ions J. Phys. Chem. C 112 1507.
[37] Cygan R T, Liang J L and Kalinichev A G 2004 Molecular models of hydroxide,
oxyhydroxide, and clay phases and the development of a general force field J. Phys.
Chem. B 108 1255.

24

[38] Hughes Z E, Carrington L A, Raitieri P and Gale J D 2011 A computational investigation

into the suitability of purely siliceous zeolites as reverse osmosis membranes J. Phys.
Chem. C 115 4063.
[39] Dick B G Jr. and Overhauser A W 1958 Theory of the dielectric constants of alkali halide
crystals Phys. Rev. 112 90.
[40] Lewis G V and Catlow C R A 1985 Potential models for ionic oxides J. Phys. C: Solid
State Phys. 18 1149.
[41] Bush T S, Gale J D, Catlow C R A and Battle P D 1994 Self-consistent interatomic
potentials for the simulation of binary and ternary oxides J. Mater. Chem. 4 831.
[42] Minervini L, Grimes R W and Sickafus K E 2000 Disorder in pyrochlore oxides J. Am.
Ceram. Soc. 83 1873.
[43] Ercolessi F and Adams J B 1994 Interatomic potentials from first-principles calculations:
the force-matching method Europhys. Lett. 26 583.
[44] Landa A, Wynblatt P, Siegel D J, Adams J B, Mryasov ON and Liu X.-Y. 2000
Development of glue-type potentials for the Al-Pb system: phase diagram calculation
Acta. Mater. 48 1753.
[45] Lenosky T J, Sadigh B, Alonso E, Bulatov V V, de la Rubia T D, Kim J, Voter A F and
Kress J D 2000 Highly optimized empirical potential model of silicon Modelling Simul.
Mater. Sci. Eng. 8 825.
[46] Lin H and Truhlar D G 2007 QM/MM: what have we learned, where are we, and where do
we go from here? Theor. Chem. Acc. 117 185.
[47] Kerdcharoen T, Liedl K R and Rode B M 1996 A QM/MM simulation method applied to
the solution of Li+ in liquid ammonia Chem. Phys. 211 313.
[48] Huang P and Carter E A 2008 Advances in correlated electronic structure methods for
solids, surfaces, and nanostructures Ann. Rev. Phys. Chem. 59 261.
[49] Fichthorn K A and Weinberg W H 1991 Theoretical foundations of dynamical Monte Carlo
simulations J. Chem. Phys. 95 1090.
[50] Boettinger W J, Warren J A, Beckermann C and Karma A 2002 Phase field simulation of
solidification Ann. Rev. Mater. Res. 32 163.
[51] Sickafus K E, Minervini L, Grimes R W, Valdez J A, Ishimaru M, Li F, McClellan K J and
Hartmann T 2000 Radiation tolerance of complex oxides Science 289 748.
[52] Sickafus K E, Grimes R W, Valdez J A, Cleave A, Tang M, Ishimaru M, Corish S M,
Stanek C R and Uberuaga B P 2007 Radiation-induced amorphization resistance and
radiation tolerance in structurally related oxides Nat. Mater. 6 217.
[53] Chartier A, Catillon G and Crocombette J P 2009 Key role of the cation interstitial
structure in the radiation resistance of pyrochlores Phys. Rev. Lett. 102 155503.
[54] Crocombette J P and Chartier A 2007 Molecular dynamics studies of radiation induced
phase transitions in La2Zr2O7 pyrochlore Nucl. Instrum. Meth. B 255 158.
[55] Devanathan R, Weber W J and Gale J D 2010 Radiation tolerance of ceramicsinsights
from atomistic simulation of damage Energy Environ. Sci. 3 1551.
[56] Devanathan R and Weber W J 2008 Dynamic annealing of defects in irradiated zirconiabased ceramics J. Mater. Res. 23 593.
[57] Tuller H L 1992 Mixed ionic-electronic conduction in a number of fluorite and pyrochlore
compounds Solid State Ionics 52 135.
25

[58] Rushton M J D, Stanek C R, Cleave A R, Uberuaga B P, Sickafus K E and Grimes R W

2007 Simulation of defects and defect processes in fluorite and fluorite-related oxides:
Implications for radiation tolerance Nucl. Instrum . and Meth. B 255 151.
[59] Devanathan R and Weber W J 2005 Insights into the radiation response of pyrochlores
from calculations of threshold displacement events J. Appl. Phys. 98 086110.
[60] Smith R, Bacorisen D, Uberuaga B P, Sickafus K E, Ball J A and Grimes R W 2005
Dynamical simulations of radiation damage in magnesium aluminate spinel, MgAl2O4
J. Phys.: Condens. Matter 17 875.
[61] Trachenko K, Dove M T and Salje E K H 2001 Atomistic modelling of radiation damage in
zircon J. Phys.: Condens. Matter 13 1947.
[62] Trachenko K, Dove M T and Salje E K H 2002 Structural changes in zircon under a-decay
irradiation Phys. Rev. B 65 180102.
[63] Trachenko K, Pruneda M, Artacho E and Dove MT 2004 Radiation damage effects in the
perovskite CaTiO3 and resistance of materials to amorphization Phys. Rev. B 70
134112.
[64] Trachenko K, Dove M T, Geisler T, Todorov I and Smith B 2004 Radiation damage effects
and percolation theory J. Phys.: Condens. Matter 16 S2623.
[65] Trachenko K, Pruneda M, Artacho E and Dove M T 2005 How the nature of the chemical
bond governs resistance to amorphization by radiation damage Phys. Rev. B 71 184104.
[66] Trachenko K, Todorov I T, Dove M T, Artacho E and Smith W 2006 Atomistic simulations
of resistance to amorphization by radiation damage Phys. Rev. B 73 174207.
[67] Devanathan R, Corrales L R, Weber W J, Chartier A and Meis C 2006 Molecular dynamics
simulation of energetic uranium recoil damage in zircon Molec. Simul. 32 1069.
[68] Devanathan R and Weber W J 2007 Radiation effects in a model ceramic for nuclear waste
disposal JOM 59 32.
[69] Devanathan R 2009 Radiation damage evolution in ceramics Nucl. Instrum. and Meth. B
267 3017.
[70] Yu J, Devanathan R and Weber W J 2009 Unified interatomic potential for zircon, zirconia
and silica systems J. Mater. Chem. 19 3923.
[71] Heinisch HL and Weber W J 2005 Computational model of alpha-decay damage
accumulation in zircon Nuclear Instrum. Meth. B 228 293.
[72] Rong Z, Gao F, Weber W J and Hobler G 2007 Monte Carlo simulations of defect recovery
within a 10 keV collision cascade in 3CSiC J. Appl. Phys. 102 103508.
[73] Pakarinen O H, Djurabekova F and Nordlund K 2010 Density evolution in formation of
swift heavy ion tracks in insulators Nuclear Instrum. Meth. B 268 3163.
[74] Crocombette, J P 2009 Can thermal spike calculations reproduce displacement cascades?
Nucl. Instrum. Meth. B 267 3152.
[75] Zhang J, Lang M, Ewing R C, Devanathan R, Weber W J and Toulemonde M 2010
Nanoscale phase transitions under extreme conditions within an ion track J. Mater. Res.
25 1344.
[76] Harding J H and Harris D J 2001 Simulation of grain-boundary diffusion in ceramics by
kinetic Monte Carlo Phys. Rev. B 63 094102.
[77] Bates J K, Bradley J P, Teetsov A, Bradley C R and Buchholtz ten Brink M 1992 Colloid
formation during waste form reaction: Implications for nuclear waste disposal Science
256 649.
26

[78] Nangia S and Garrison B J 2010 Theoretical advances in the dissolution studies of mineral
water interfaces Theor. Chem. Acc. 127 271.
[79] Xiao Y and Lasaga A C 1996 Ab initio quantum mechanical studies of the kinetics and
mechanisms of quartz dissolution: OH- catalysis Geochim. Cosmochim Acta 60 2283.
[80] Criscenti L J, Kubicki J D and Brantley S L 2006 Silicate glass and mineral dissolution:
calculated reaction paths and activation energies for hydrolysis of a Q3 Si by H3O+
using ab initio methods J Phys. Chem. A 110 198.
[81] Du Z and de Leeuw N H 2006 Molecular dynamics simulations of hydration, dissolution
and nucleation processes at the -quartz (0001) surface in liquid water Dalton Trans.
2623.
[82] Mahadevan T S and Garofalini S H 2008 Dissociative chemisorption of water onto silica
surfaces and formation of hydronium ions J. Phys. Chem. C 112 1507.
[83] Lasaga A C and Luttge A 2005 Kinetic Justification of the solubility product: Application
of a general kinetic dissolution model J. Phys. Chem. B 109 1635.
[84] Nangia S and Garrison B J 2009 Advanced Monte Carlo approach to study evolution of
quartz surface during the dissolution process J. Am. Chem. Soc. 131 9538.
[85] Aertsens M 2009 Long-term behaviour of nuclear waste glass: Monte Carlo simulations
Scientific report for WP5 Task 1 of RP.WD.008 SCKCEN-ER-109.
[86] Hu S, Li Y, Sun X, Gao F, Devanathan R, Henager C H and Khaleel M 2010 Application
of the phase field method in predicting gas bubbles microstructure evolution in nuclear
fuels Int. J. Mat. Res. 101 515.
[87] Li Y, Hu S, Sun X, Gao F, Henager C H and Khaleel M 2010 Phase-field modeling of void
migration and growth kinetics in materials under irradiation and temperature field J.
Nucl. Mater. 407 119.
[88] Li Y, Hu S, Sun X, Gao F, Henager C H and Khaleel M 2011 Phase-field modeling of void
evolution and swelling in materials under irradiation, SCIENCE CHINA: Physics,
Mechanics & Astronomy 54 856.
[89] Rokkam S, El-Azab A, Millett P and Wolf D 2009 Phase field modeling of void nucleation
and growth in irradiated metals Modelling Simul. Mater. Sci. Eng. 17 064002.
[90] Devanathan R, Sickafus K E, Weber W J and Nastasi M 1998 Effects of ionizing radiation
in ceramics J. Nucl. Mater. 253 113.
[91] Xiao H Y, Gao F and Weber W J 2010 Threshold displacement energies and defect
formation energies in Y2Ti2O7 J. Phys.: Condens. Matter 22 415801.
[92] Lucas G and Pizzagalli L 2007 Ab initio molecular dynamics calculations of threshold
displacement energies in silicon carbide Phys. Rev. B 72 161202.
[93] Gao F, Xiao H Y, Zu X T, Posselt M and Weber W J 2009Defect-enhanced charge transfer
by ion-solid interactions in sic using large-scale ab initio molecular dynamics
simulations Phys. Rev. Lett. 103 027405.
[94] Sutton A P, Todorov T N, Cawkwell M J and Hoekstra J 2001 A simple model of atomic
interactions in noble metals based explicitly on electronic structure Phil. Mag. 81 1833.
[95] Zhou X W and Doty F P 2008 Embedded-ion method: An analytical energy-conserving
charge-transfer interatomic potential and its application to the La-Br system Phy. Rev.
B 78 224307.

27

[96] Tschopp M A, Solanki K N, Baskes M I, Gao F, Sun X and Horstemeyer M F 2011

Generalized framework for interatomic potential design: Application to Fe-He system
J. Nucl. Mater. doi:10.1016/j.jnucmat.2011.08.003.

28