val-2g

Deuterium Transport in Proton-Conducting Ceramics

Case Description

This case reproduces the analysis published in Yang et al. (2026), which proposed a new model that accurately reproduces deuterium transport in a proton-conducting ceramics in both dry and wet environment. This validation effort involves deuterium gas (D) and heavy water (DO) transport in proton-conducting ceramics (PCC), specifically yttrium-doped barium zirconate (BaZrYO, also known as BZY10). The experimental data used for validation comes from thermal desorption spectroscopy (TDS) measurements performed by Hossain et al. (Hossain and Hashizume, 2022). The primary objective is to understand and model the mechanisms of hydrogen isotope transport in PCC materials, which are of interest for tritium extraction systems in fusion energy applications.

In proton-conducting ceramics, hydrogen isotopes are transported as protonic defects (OD) through a hopping mechanism between oxygen sites in the crystal lattice. Unlike metals where hydrogen diffuses as interstitial atoms, hydrogen in PCC is incorporated into the oxygen sublattice. This case models the TDS experiment where a BZY10 sample (0.5 mm thick, 7.7 mm 2.2 mm surface area) is exposed to either D gas (dry condition) or DO vapor (wet condition), followed by thermal desorption analysis.

The TDS simulation consists of three phases:

  1. Dissolution phase (1 hour at 873 K): The sample is exposed to D at 1.33 kPa (dry condition) or DO at 2.8 kPa (wet condition), allowing deuterium to dissolve into the material.

  2. Cooldown phase (1 hour): The sample temperature decreases exponentially from 873 K to approximately 300 K while the gas pressure is reduced to vacuum levels.

  3. Desorption phase: The sample is heated from 300 K to 1400 K at a constant rate of 0.5 K/s, causing trapped deuterium to be released.

The temperature and pressure histories are shown in Figure 1.

Temperature and pressure history during the TDS simulation for the dry (D$_2$) and wet (D$_2$O) condition.

Figure 1: Temperature and pressure history during the TDS simulation for the dry (D) and wet (DO) condition.

Model Description

Diffusion of Mobile Species

In BZY10, three mobile species are considered: protonic defects (OD), oxygen vacancies (V), and electrons (e). The hydrogen isotope transport model for PCCs considers diffusion, trapping within the material, and kinetic chemical reactions on the surface. These transport mechanisms are illustrated in Figure 2.

Schematic of the hydrogen isotope transport in PCCs.

Figure 2: Schematic of the hydrogen isotope transport in PCCs.

The governing equations for OD in PCC materials is described as

where is the concentration of mobile OD, is the diffusivity, is the concentration of trapped OD in the material. The diffusivity follows an Arrhenius relationship:

where is the pre-exponential factor, is the activation energy, is the gas constant, and is temperature. The governing equations for the diffusion of the other mobile species, V and e, are the same as those for OD, except that no trapping effects are considered.

Trapping and Detrapping

The model includes a single trapping site to capture deuterium retention effects observed in the TDS spectra. The trapped concentration evolves according to:

where is the lattice site density, and is the empty trap concentration with being the trap site fraction.

The trapping and release rate coefficients follow Arrhenius relationships:

where and are pre-factors of trapping and release rate coefficients, and are trapping and release energies, and is the Boltzmann constant.

Surface Reactions

In addition to diffusion and trapping within the bulk material, hydrogen isotope transport also involves chemical reactions at the surface. The chemical reactions using the Kröger–Vink notation are:

DO hydration reaction (wet condition):

D incorporation reaction (dry condition):

The net flux for each reaction is:

where and are the pressure and concentration of DO or D, respectively, and are the forward reaction rate and the reverse reaction rate for reaction , respectively. At equilibrium, the surface concentrations depend on the chemical reaction constant, which is a function of temperature:

(1)

and

(2)

where is the chemical reaction constant of reaction , and are the entropy and enthalpy of reaction , respectively.

commentnote:The comprehensive mechanism and analysis is in Yang et al. (2026).

This model description is a simplified version of the Method section in Yang et al. (2026). The comprehensive explanation of mechanism and analysis can be found in Yang et al. (2026).

Case and Model Parameters

Table 1 summarizes the detail of sample and experimental conditions from Hossain et al. (Hossain and Hashizume, 2022), as well as the model parameters from Hossain et al. (Hossain and Hashizume, 2022), Kreuer (Kreuer, 2003), and estimated from validation cases in TMAP8. Table 2 includes the trapping parameters from Karmonik et al. (Karmonik et al., 1995) and estimated based on existing validation cases in TMAP8.

Table 1: Experimental set up from Hossain et al. (Hossain and Hashizume, 2022) and modeling parameters for deuterium transport in a BZY membrane at 873 K from Hossain et al. (Hossain and Hashizume, 2022) and from Kreuer et al. (Kreuer, 1999; Kreuer, 2003).

ParameterDescriptionValueUnitsReference
Initial/dissolution temperature873KHossain and Hashizume (2022)
Cooldown final temperature300KHossain and Hashizume (2022)
Desorption final temperature1400KHossain and Hashizume (2022)
Heating rate0.5K/sHossain and Hashizume (2022)
Dissolution duration1hHossain and Hashizume (2022)
Cooldown duration1hHossain and Hashizume (2022)
Sample thickness0.5mmHossain and Hashizume (2022)
Sample surface area16.94mmHossain and Hashizume (2022)
Sample density5.98g/cmHossain and Hashizume (2022)
Initial electron concentrationg/cmHossain and Hashizume (2022)
D pressure (dry condition)1.33kPaHossain and Hashizume (2022)
DO pressure (wet condition)2.8kPaHossain and Hashizume (2022)
Entropy for DO reaction-88.90J/mol/KKreuer (2003)
Enthalpy for DO reactionJ/molKreuer (2003)
Entropy for D reaction-124.53J/mol/KEstimated from Kreuer (2003)
Enthalpy for D reactionJ/molEstimated from Kreuer (2003)
Forward reaction rate for DOm/atom/sEstimated from TMAP8
Forward reaction rate for Dm/atom/sEstimated from TMAP8
Diffusivity coefficient for ODm/sHossain and Hashizume (2022)
Diffusivity energy for OD0.23eVHossain and Hashizume (2022)
Diffusivity coefficient for Vm/sKreuer (1999)
Diffusivity energy for V89216.77J/molKreuer (1999)
Diffusivity coefficient for em/sKreuer (1999)
Diffusivity energy for e103818.22J/molKreuer (1999)

Table 2: Trapping parameters for deuterium transport in BZY membrane from Karmonik et al. (Karmonik et al., 1995) and estimated from TMAP8.

ParameterDescriptionValueUnitsReference
Trapping rate coefficient1/sKarmonik et al. (1995)
Release rate coefficient1/sKarmonik et al. (1995)
Trapping energy0.38eVEstimated from TMAP8
Release energy1.60eVEstimated from TMAP8
Trapping site atom fraction-Estimated from TMAP8

Results

Results Without Trapping Model

Initial simulations were performed without the trapping model using diffusivity and surface reaction parameters from literature. Figure 3 shows the comparison between TMAP8 predictions and experimental data for the dry condition (D exposure), while Figure 4 shows results for the wet condition (DO exposure).

As shown in Figure 5 and Figure 6, the model does not capture the experimental data from Hossain et al (Hossain and Hashizume, 2022). The fluxes of D and DO under both dry and wet environments exhibit significant discrepancies between the simulations and experimental results, except for the DO flux under wet environment. The temperature at which the deuterium flux peak occurs for D under dry condition is drastically lower than that of the experimental peak, and the predicted peak is a lot wider than experimentally observed. Furthermore, the deuterium flux from DO under dry environment and from D under wet environment remains close to 0 and far from the experimental results. The root mean square percentage error (RMSPE) values exceed 100% for most cases, indicating that this model without traps does not accurately reproduce the experimental results.

Comparison of TMAP8 calculations (without trapping) with experimental data for the dry condition. RMSPE values are shown on the plot.

Figure 3: Comparison of TMAP8 calculations (without trapping) with experimental data for the dry condition. RMSPE values are shown on the plot.

Comparison of TMAP8 calculations (without trapping) with experimental data for the wet condition. RMSPE values are shown on the plot.

Figure 4: Comparison of TMAP8 calculations (without trapping) with experimental data for the wet condition. RMSPE values are shown on the plot.

Results With Trapping Using Initial Parameters

Adding the trapping model with parameters from literature improves the predictions, as shown in Figure 5 and Figure 6. The trapping model captures the delayed release of deuterium and produces peak shapes more consistent with the experimental data. However, the peak positions and magnitudes still show significant discrepancies.

Comparison of TMAP8 calculations (with trapping, initial parameters) with experimental data for the dry condition. RMSPE values are shown on the plot.

Figure 5: Comparison of TMAP8 calculations (with trapping, initial parameters) with experimental data for the dry condition. RMSPE values are shown on the plot.

Comparison of TMAP8 calculations (with trapping, initial parameters) with experimental data for the wet condition. RMSPE values are shown on the plot.

Figure 6: Comparison of TMAP8 calculations (with trapping, initial parameters) with experimental data for the wet condition. RMSPE values are shown on the plot.

Results With Trapping Using Calibrated Parameters

To improve the agreement with experimental data, the model parameters were calibrated using the Parallel Subset Simulation (PSS) approach in MOOSE's stochastic tools module. The PSS method systematically explores the parameter space to minimize the difference between simulation predictions and experimental measurements. A detailed explanation of this method is provided in Yang et al. (2026). In this PSS optimization process, an optimization function using RMSPEs is calculated from the four flux results: DO and D fluxes under both dry and wet environments. We calculate the difference from the four fluxes during desorption simultaneously, since the experiments used the same BZY sample. The multi-objective optimization function (log inverse error) is described as:

where is the RMSPE for or fluxes under dry or wet environment, and is the results from simulations or experiments. A high means a better accuracy of the simulation results. Figure 7 shows the evolution of optimization function from a simplified PSS optimization process using initial parameters close to the calibrated values.

The evolution of the objective function (log inverse error) during a simplified PSS optimization.

Figure 7: The evolution of the objective function (log inverse error) during a simplified PSS optimization.

Figure 8 shows the deviation of all calibrated parameters from the input parameters. All deviations are within three standard deviations of their corresponding normal distribution ranges.

Deviation of calibrated material parameters from initial values, expressed in standard deviations of their respective normal distribution ranges.

Figure 8: Deviation of calibrated material parameters from initial values, expressed in standard deviations of their respective normal distribution ranges.

The calibrated parameters are listed in Table 3, Table 4, and Table 5. These calibrated parameters are taken directly from Yang et al. (2026).

Table 3: Calibrated diffusivity parameters.

ParameterDescriptionValueUnitsReference
OD diffusivity pre-factor1.90 m/sYang et al. (2026)
OD activation energy0.122eVYang et al. (2026)
V diffusivity pre-factor1.24 m/sYang et al. (2026)
V activation energy1.04eVYang et al. (2026)
e diffusivity pre-factor2.06 m/sYang et al. (2026)
e activation energy0.99eVYang et al. (2026)

Table 4: Calibrated trapping parameters.

ParameterDescriptionValueUnitsReference
Trap site fraction1.39 -Yang et al. (2026)
Trapping rate pre-factor3.95 1/sYang et al. (2026)
Trapping energy0.467eVYang et al. (2026)
Release rate pre-factor1.60 1/sYang et al. (2026)
Release energy1.24eVYang et al. (2026)

Table 5: Calibrated thermodynamic parameters.

ParameterDescriptionValueUnitsReference
DO reaction enthalpy-156.4kJ/molYang et al. (2026)
DO reaction entropy-137.4J/mol/KYang et al. (2026)
DO forward rate constant2 m/at/sYang et al. (2026)
D reaction enthalpy-112.2kJ/molYang et al. (2026)
D reaction entropy-37.0J/mol/KYang et al. (2026)
D forward rate constant2 m/at/sYang et al. (2026)

With these calibrated values, TMAP8 achieves good agreement with the experimental TDS data for both dry and wet conditions, as shown in Figure 9 and Figure 10. The root mean square percentage errors (RMSPE) are 16.59%, 20.49%, 37.60%, and 40.49% for the D flux (dry), DO flux (dry), D flux (wet), and DO flux (wet), respectively.

Comparison of TMAP8 calculations (with trapping, calibrated parameters) with experimental data for the dry condition. RMSPE values are shown on the plot.

Figure 9: Comparison of TMAP8 calculations (with trapping, calibrated parameters) with experimental data for the dry condition. RMSPE values are shown on the plot.

Comparison of TMAP8 calculations (with trapping, calibrated parameters) with experimental data for the wet condition. RMSPE values are shown on the plot.

Figure 10: Comparison of TMAP8 calculations (with trapping, calibrated parameters) with experimental data for the wet condition. RMSPE values are shown on the plot.

The calibrated model successfully reproduces the key features of the TDS spectra:

  • The peak temperature positions for both D and DO desorption

  • The relative magnitudes of the desorption peaks

  • The peak shapes and widths

This validation demonstrates TMAP8's capability to model hydrogen isotope transport in proton-conducting ceramics in both dry and wet environments, including the complex interplay between diffusion, trapping, and surface reactions.

Input files

The input files for this validation case are:

More information about these tests can be found in the test specification file for this case, namely (tmap8/test/tests/val-2g/tests).

References

  1. M Khalid Hossain and Kenichi Hashizume. A comparative study on the hydrogen dissolution and release behaviors in the zirconate proton conductors by TDS and TMAP4 analysis. Journal of Alloys and Compounds, 907:164436, 2022.[BibTeX]
  2. Ch Karmonik, R Hempelmann, Th Matzke, and T Springer. Proton diffusion in strontium cerate ceramics studied by quasielastic neutron scattering and impedance spectroscopy. Zeitschrift für Naturforschung A, 50(6):539–548, 1995.[BibTeX]
  3. KD Kreuer. Aspects of the formation and mobility of protonic charge carriers and the stability of perovskite-type oxides. Solid State Ionics, 125(1-4):285–302, 1999.[BibTeX]
  4. Klaus-Dieter Kreuer. Proton-conducting oxides. Annual Review of Materials Research, 33(1):333–359, 2003.[BibTeX]
  5. Lin Yang, Pierre-Clément A Simon, Wei Tang, Meng Li, Zeyu Zhao, Dong Ding, and Thomas Fuerst. Elucidating hydrogen isotope transport mechanisms in proton-conducting ceramics with trapping effects using TMAP8. International Journal of Hydrogen Energy, 210:153551, 2026. doi:10.1016/j.ijhydene.2026.153551.[BibTeX]