Data compilation

The datasets were gathered from the JAEA-DDB and 16 published resources, covering the period from 1982 to 2024. The dataset comprised 16 input features and 324 experimental instances, including 304 instances obtained from Wu et al. (2024) and 20 experimental instances from three other studies [17, 20, 27]. Notably, the absence of pH values in 514 instances of the JAEA-DDB resulted in a significantly reduction in data size. To address this, regression imputation techniques using ML models were applied to predict the pH values based on a dataset of 324 instances, thereby expanding the dataset to 838 instances.

The dataset included 16 input features, which were categorized into three groups: (i) porewater properties, comprising the ionic strength (I), temperature (T), and pH; (ii) bentonite properties, including the montmorillonite content (m), external surface area (A_ext), dry density (ρ_d), grain density (ρ_s), total porosity (ε_tot, and montmorillonite stacking number (n_c); and (iii) radionuclide properties, encompassing the ion diffusion coefficient in water (D_w), molecular weight (MW), ion molar conductivity (λ), ionic radius (r), ionic charge (z), distribution coefficient (K_d), and rock capacity factor (α).

Data preprocessing

The presence of outliers can reduce the predictive accuracy of ML models. To address this issue, the Mahalanobis distance (MD) method was employed to identify and remove outliers. The cutoff point (di) is given as: $d_i=\sqrt{(x-\mu )\cdot S^{-1}\cdot (x-\mu )}\mathrm{,}$ (1) where x represents the object vector, μ denotes the mean arithmetic vector, and S is the covariance matrix of instances. The cutoff point was set to eight to ensure that the skewness of all input features was less than 10.

Three datasets were used to enhance the prediction of radionuclide diffusion. An overview of the features and instances of each dataset is summarized in Table 2. Dataset I included 15 input features, with pH as the output feature. To ensure the data quality and reduce noise, eight instances were removed using the MD method. This process yielded Dataset I, comprising 316 instances. The statistical details of Dataset I are presented in Table S1 of the supporting information. Datasets II and III comprised 16 input features, including the basic features (15 input features of Dataset I) and pH. The output feature for Datasets II and III was the D_e. Dataset III, comprising 813 instances, was obtained after removing 17 instances. It is noteworthy that these datasets comprised parameters at the micro-mesoscopic level. Specifically, the montmorillonite stacking number and ionic radius were classified as microscopic parameters, whereas the other parameters were considered as mesoscopic.

Imputation methods

Four decision-tree models, namely LGBM, CatB, XGB, and RF, were used as regression imputation methods to predict the pH values of Dataset I. LGBM exhibited superior predictive accuracy compared with the other models. This was consistent with the results of our previous study [21]. Dataset III was established by incorporating additional 514 instances with Dataset II using the LGBM for data imputation. Table S2 of the supporting information summarizes the statistical results of the input and output features for Dataset III.

Regression imputation for predicting pH

Handling missing data is a crucial step affecting the quality and reliability of the data analysis. Various regression imputation techniques have been applied to impute missing data, such as ANNs, multivariate imputation by chained equations, k-nearest neighbors, time-series deep learning models, generative broad Bayesian imputation, principal component analysis imputation, and simple arithmetic averages. These methods have been applied to datasets with missing data percentages ranging from 0 to 80% [24, 26, 32-36]. Generally, three types of missing data mechanisms are recognized: missing completely at random, missing at random, and missing not at random [23]. Each mechanism presents different challenges and implications for imputation, highlighting the importance of identifying the underlying pattern of missingness before selecting an appropriate imputation strategy.

The JAEA-DDB database collected data from the literatures and reports covering 1982 to 2009. The instances have been derived from various diffusion experimental methods and numerous researchers. The absence of pH values in 514 instances within the JAEA-DDB database can be explained that these researches ignored the importance of pH values in their studies. In the JAEA-DDB database, missing data primarily resulted from ignoring or inadequately measuring the parameters that related to the radionuclide diffusion. The missing mechanism in the JAEA-DDB database was assumed to be missing completely at random, corresponding to a non-continuous missingness. Based on the selected 16 input features, more than 60% of the dataset (514 instances) lacked pH values. Decision-tree models were employed to predict the missing pH values to augment the dataset and enhance the robustness of the ML models. Specifically, LGBM, CatB, XGB, and RF were employed to predict the pH values of Dataset I.

The predicted performances are summarized in Table 3. The LGBM exhibited superior robustness compared with the other models. For instance, the $R_{\text{cv}}^2$ values for the test sets were ranked in descending order using a five-fold cross-validation as follows: LGBM > XGB > CatB > RF. The rank of MSEcv values was the opposite of that of the $R_{\text{cv}}^2$ values for the test datasets. Notably, LGBM achieved the highest performance metrics among all models, with an MSE of 0.23 and R2 of 0.92 for the test dataset. The hyperparameters of the optimal ML models are listed in Table S3 of the supporting information. Therefore, the missing pH values for 514 instances were predicted using the LGBM model, resulting in the establishment of Dataset III with 813 instances.

Table 3Full-screen View

Details of the features and instances of datasets

Show more

Figure 2 exhibits the data distribution and characteristics of the relationship between pH and each input feature. Blue and orange represent the data distributions of Dataset I and the imputed 514 instances, respectively. It clearly demonstrates a non-linear relationship between the pH and each input feature. The predicted pH values ranged from 5.0 to 9.0, exhibiting a Gaussian type distribution.

Fig. 2Full-screen View

(Color online) Data distribution of features and the relationship between pH and each input feature

pH is an important porewater parameter that influences both the radionuclide species and surface charge of clay [37]. Figure 3 shows the pH dependence on the external surface area and ion molar conductivity, which are associated with the bentonite and radionuclide properties, respectively. Dataset I exhibits that the pH values ranged from 3.0 to 13.4. The predicted pH values are concentrated in the range from 5.0 to 9.0, suggesting a close adherence to a normal distribution of porewater for Dataset III.

Fig. 3Full-screen View

(Color online) Analyzing the dependency of pH on the external surface area and ion molar conductivity

Model development for radionuclide diffusion

Ten ML models, namely LGBM-CatB, LGBM-XGB, LGBM-RF, LGBM, CatB, XGB, RF, ANN, DNN, and SVM, were used to predict the D_e values of radionuclides in compacted bentonite. Figure 4 shows the performance metrics of the ML models for the test datasets of Dataset II and III using the optimal hyperparameters tuned with PSO techniques (Table S4 in the supporting information). The performance metrics were assessed using five-fold cross-validation. The red lines represent the smooth kernel curve of the distribution of performance metrics. The black lines within and outside the box plots denote the mean values and standard deviations of the performance metrics, respectively, with a lower standard deviation indicating strong robustness of the ML models. The detailed performance metrics for the training, validation, and test datasets are listed in Table S5 of the supporting information.

Fig. 4Full-screen View

Mean performance metric values using five-fold cross-validation for machine learning models in the test datasets of Dataset II and III

As the number of instances increased from 316 (Dataset II) to 813 (Dataset III), the performance metrics of all ML models improved significantly, as evidenced by the higher $R_{\text{cv}}^2$ values, lower MSEcv, and reduced standard deviation. These findings indicate that expanding the dataset contributed to enhanced predictive performance and robustness of the ML models. It is noteworthy that the ensemble models were established by combining LGBM with other individual decision-tree models, primarily due to the relatively high training speed of the LGBM algorithm [38]. However, no significant difference is observed in the computational efficiencies of the ensemble and single models. The difference in running time was approximately five minutes. In the case of decision-tree algorithms, gradient boosting (GB) models (LGBM, CatB, and XGB) outperformed the RF models. The excellent predictive performance of GB models is consistent with previous findings in predicting the chloride diffusion coefficient in concrete [39]. In addition, the ensemble ML models (LGBM-CatB, LGBM-XGB, and LGBM-RF) and LGBM surpassed the other ML models, achieving an $R_{\text{cv}}^2$ above 0.90. This can be attributed to their ability to harness the strengths of various algorithms to thoroughly capture potentially complex patterns and errors within the data, thereby enhancing the prediction accuracy and robustness [30, 40]. For Dataset III, the $R_{\text{cv}}^2$ values of the ML models ranked in descending order as follows: LGBM-CatB ≈ LGBM-XGB > LGBM ≈ LGBM-RF > CatB ≈ XGB > ANN > DNN > RF > SVM. Notably, LGBM-CatB surpassed LGBM-XGB due to its lower standard deviation, indicating stronger robustness. SVM exhibited the lowest predictive performance based on Dataset III, with $R_{\text{cv}}^2\text{=}0.75$ and MSE_cv = 0.06. Compared with ensemble models, SVM is a relatively simple model. The ensemble models are designed to capture more complex patterns and relationships in the data through a combination of multiple decision trees. This lack of complexity in the SVM limits its ability to generalize across different data instances in the dataset. Notably, some studies have reported test R² values below 0.80, such as an R² of 0.74 for predicting the retention rate of Cd in biochar [41] and an R² of 0.76 for predicting alcohol space-time yield [42]. Therefore, the prediction accuracy of SVM remained satisfactory, despite exhibiting a lower predictive performance than the other models.

Figure 5 shows the regression plots comparing the experimental and predicted D_e values for the training (green triangle), validation (red circle), and test (purple square) datasets of Datasets II and III, using the LGBM-CatB, LGBM-XGB, LGBM, and LGBM-RF algorithms. These algorithms were selected owing to their excellent predictive accuracies. The plots reveal a close alignment between the experimental and predicted D_e values with the slope line, underscoring the effective simulation capability of these ML models for predicting radionuclide diffusion processes. The performance metrics of the best-performing models are shown in Fig. 5. Notably, the ML models applied to the test dataset of Dataset III outperformed those applied to Dataset II. This disparity can be attributed to the augmentation of instances in Dataset III, which facilitates the models’ ability to capture complex relationships within the data more effectively. For Dataset III, the ranking of models was as follows: LGBM-CatB (R² = 0.94) ≈ LGBM-XGB (R² = 0.94) > LGBM (R² = 0.92) ≈ LGBM-RF (R² = 0.92). These results indicate that both LGBM-CatB and LGBM-XGB exhibit high predictive accuracy.

Fig. 5Full-screen View

(Color online) Regression plots of experimental versus predicted effective diffusion coefficients based on Datasets II and III: (a, e) LGBM-CatB, (b, f) LGBM-XGB, (c, g) LGBM, and (d, h) LGBM-RF

Spearman and Shapley additive explanation analyses

ML models can uncover predictive principles through analysis techniques that rank the importance of influencing factors in predictions, such as feature importance and SHAP analysis [19, 21, 43, 44]. Additionally, Spearman analysis, a non-parametric statistical method, assesses the monotonic relationship between two variables by correlating ranked data. These approaches provided valuable insights into the consistency and strength of the relationships within a dataset. It worthy notes that the reliability of these analytical techniques is intrinsically linked to the quality of the data used. Increasing the dataset size enhances the depth, broadness, and reliability of the ML models.

Spearman correlation and SHAP analysis techniques were employed to analyze the correlation and importance of the input features, presenting intuitively global interpretations of the ML models (Fig. 6). The features were ranked from left to right according to their correlation and contribution to the prediction. The Spearman correlation analysis revealed that the most influential factor among the 16 input features was the ion diffusion coefficient in water for Dataset II, and the total porosity for Dataset III. This feature exhibited a positive correlation with D_e (Figs. 6a and b). This is consistent with the previous findings [19] and Archie’s law [31, 45].

Fig. 6Full-screen View

(Color online) (a, b) Spearman correlation analysis and global interpretations of ML models based on Dataset II and III: (c, d) LGBM-CatB, (e, f) LGBM-XGB, and (g, h) LGBM

In the case of Dataset II, the SHAP analysis revealed that the most important input features varied across different ML models: the compacted dry density for LGBM-CatB, ionic radius for LGBM-XGB, and ion diffusion coefficient in water for LGBM (Figs. 6c, e, and g). Notably, only the SHAP results for LGBM were consistent with the Spearman correlation analysis. This discrepancy can be attributed to differences in the feature importance assessment and prediction mechanisms inherent to each ML algorithm. As the number of instances increased from 316 (Dataset II) to 813 (Dataset III), both Spearman and SHAP analyses identified the total porosity as the primary contributor, which is consistent with Archie’s law [31, 45]. The total porosity for radionuclide diffusion in compacted bentonite blocks is expressed as a percentage of the total interconnected pore space within the blocks. A higher total porosity implies greater availability of transport pathways. These findings suggest that larger datasets may reduce the discrepancies between ML models in terms of feature importance assessment and prediction mechanisms.

Partial dependence plots

The dependence of D_e on the 16 input features has been discussed in our previous study [19]. However, some relationships may remain unclear due to the limited size of the dataset. To address this, PDP analysis was performed to visually represent the univariate correlations and examine the influence of the size of the dataset on these relationships (Fig. 7). The histograms and lines correspond to the data distribution and correlation with each input feature and the PDP. Generally, a more concentrated data distribution generally leads to more accurate analytical results. These findings indicate that Dataset III, which was larger than Dataset II, exhibited more continuous PDP curves, suggesting a more stable and clear relationship between the features and D_e.

Fig. 7Full-screen View

(Color online) Partial dependence plot for (a) the rock capacity factor, (b) distribution coefficient, (c) pH, (d) ionic charge, (e) ion molar conductivity, (f) external surface area, (g) montmorillonite stacking number, (h) grain density, (i) ionic strength, (j) total porosity, (k) dry density, (l) montmorillonite content, (m) ion diffusion coefficient in water, (n) ionic radius, (o) molecular weight, and (p) temperature

Figures 7a and b shows that both the rock capacity factor and distribution coefficient exhibit a clear positive correlation with the prediction for Dataset III. This finding is consistent with studies on radionuclides diffusion in crystalline rocks [46] and sodium montmorillonite [47]. Consistently, Fig. 7d illustrates the positive impact of ionic charge, where cations exhibit a higher D_e than neutral species, and anions display lower D_e values. This is consistent with previous studies, which attributed the differences in diffusion mechanisms to electrostatic interactions between the radionuclide species and charged bentonite surfaces [3]. Specifically, cation diffusion is controlled by surface diffusion effects, whereas anions diffusion is driven by anionic exclusion effects [47, 48].

pH values in the range from 6 to 9 negatively influence the prediction for Dataset III, whereas a peak was observed at approximately pH 8 for Dataset II (Fig. 7c). The negative effect of Dataset III might be more convincing because of its larger data size. Figure 7e shows a positive impact on the prediction when ion molar conductivity exceeded 0.01 m²⋅ S/mol for Dataset III. However, the relationships among the external surface area, montmorillonite stacking number, grain density, and ionic strength remained unclear for both Datasets II and III (Figs. 7f-i). This lack of clarity can be attributed to data dispersion, despite the larger dataset size.

In the case of remaining input features, such as the total porosity, ion diffusion coefficient in water, and temperature, exhibited positive impacts on the prediction, whereas the dry density, montmorillonite content, ionic radius, and molecular weight showed negative impacts (Figs. 7j-p). The positive influences of the total porosity and ion diffusion coefficient in water could be explained by Archie’s law [16, 45], whereas the positive impact of temperature followed the Arrhenius equations [49-51]. The detailed explanations are provided in our previous studies [19, 21]. It is worth mentioning that a negative influence of ionic radius was observed at Logr < -9.6 (2.5 Å). This positive relationship can be attributed to the limited data for species with ionic radius above 2.5 Å. Overall, the univariate correlation results visualized using the PDP technique align with the diffusion laws observed in the experiments and diffusion mechanisms derived from the numerical models. This consistency underscores the reliability of the interpretation capabilities of the ML models.

Determination of the diffusion parameters using diffusion experiments

Figure 8 shows the breakthrough curves of EuEDTA^- and the species distribution of Eu-EDTA complexes. A_cum denotes the accumulated mass of EuEDTA^- and HCrO₄^- that penetrated a 1.2 cm thick bentonite block to reach the sample reservoirs. The data show that the accumulated mass increased with decreasing dry density, which is consistent with the general understanding that lower dry density facilitates radionuclide diffusion through porous media [3, 5]. The pH was maintained at 5.3 ± 0.1 during the Eu(III) diffusion experiments. Simulations using Vision MINTEQ indicated that Eu(III) exists as a mixture of species, including Eu³⁺, EuHEDTA(aq), EuEDTA^-, and EuCl²⁺, in 0.6 mol/L NaCl solution (Fig. 8c). EuEDTA^- was the main species at pH above 2.0. It indicates that this study measured the diffusion parameters of EuEDTA^- in compacted Ba-bentonite.

Fig. 8Full-screen View

Relationship between the accumulated mass (A_cum) and time for (a) EuEDTA^- and (b) HCrO₄^- in saturated compacted bentonites. (c) Species distribution of Eu(III)-EDTA system in aqueous solution

Table 4 summarizes the diffusion parameters of HCrO₄^- and EuEDTA^-, including D_e, rock capacity factor, accessible porosity, total porosity, and distribution coefficient. Both D_e and distribution coefficient are important parameters in the safety assessment of repositories, whereas the other parameters play a crucial role in elucidating the diffusion mechanism. The error in the compacted dry density measurement was primarily attributed to a loss of approximately 0.3 g during the preparation of bentonite blocks. Both HCrO₄^- and EuEDTA^- are monovalent anions that can’t access the interlayer pores of compacted bentonite [17, 21]. The rock capacity factor of HCrO₄^- was lower than the total porosity, indicating that the accessible porosity was equal to the rock capacity factor. This suggests that the predominant diffusion path of HCrO₄^- was through the free pores of compacted bentonite. In contrast, EuEDTA^- exhibited an adsorptive behavior similar to that of simulated trivalent actinide complexes, such as AmEDTA^- and CmEDTA^-, with the rock capacity factor being higher than the total porosity. The distribution coefficient, K_d, of EuEDTA^- was calculated as follows: $K_{\text{d}}=\frac{\alpha -{\epsilon }_{\text{acc}}}{{\rho }_{\text{d}}}\mathrm{,}$ (7) where the accessible porosity, ε_acc, was obtained using the I^- diffusion experiments [19].

All diffusion parameters decreased with increasing dry density for both EuEDTA^- and HCrO₄^-. The distribution coefficient of EuEDTA^- ranged from 4.2 × 10^-4 m³/kg to 6.7 × 10^-4 m³/kg, which is lower than the range reported for EuEDTA^- in hard rock clay (1.3 × 10^-3-3.2 × 10^3- m³/kg) [52] and for CeEDTA^- in compacted Zhisin bentonite (0.8 × 10^-3-1.2 × 10^-3 m³/kg) [17]. The distribution coefficient of EuEDTA^- was lower than that of Eu³⁺, indicating that EDTA facilitated the diffusion of Eu(III), thereby reducing the retardation capacity of the bentonite barrier [52, 53]. This observation is consistent with the diffusion behavior of CeEDTA^- and CoEDTA^2- [17, 19, 31].

Model application

The LGBM-CatB and LGBM-XGB models were employed to predict the D_e of HCrO₄^- in compacted Wyoming bentonite and EuEDTA^- in compacted Ba-bentonite, which were compared with published diffusion experimental results for HCrO₄^- and the simulated actinides CeEDTA^- and CoEDTA2^- [17, 19, 21]. Additionally, both models were used to predict the D_e of radionuclide cation ¹³⁷Cs⁺ and neutral species HTO [8, 54, 55] (Fig. 9). It shows that D_e/D_w decreased with increasing compacted dry density, which is consistent with the result of previous studies [3, 5, 45]. In this study, the D_w value for metal-EDTA complexes was assumed to be 5.0 × 10^-10 m2/s [56]. The D_e of EuEDTA^- was observed to be higher than those of CeEDTA^- [17]. and CoEDTA2^- [19]. The LGBM-CatB and LGBM-XGB models successfully predict D_e, as evidenced by the good agreement with the experimental D_e values (Fig. 9a).

Fig. 9Full-screen View

Generalization ability validation of LGBM-CatB and LGBM-XGB: (a) M-EDTA(z-4)+ diffusion, (b) HCrO₄^- diffusion, (c) ¹³⁷Cs+ diffusion, and (d) HTO diffusion

Figure 9b shows that the D_e of HCrO₄^- in compacted Wyoming bentonite is lower than that in Anji bentonite [19] and GMZ bentonite [21], likely due to the higher montmorillonite content. LGBM-CatB slightly underestimated D_e for HCrO₄^- in Wyoming bentonite, with the predicted D_e values being 25%-47% lower than the experimental D_e. Although this discrepancy is less pronounced than the predictions for HCrO₄^- in GMZ and Anji bentonites using LGBM and PSO-LGBM, the difference was reported to be 9%-27% [19, 21]. This performance is significantly superior to that predicted using Archie’s law, according to which the predictive D_e values were 1.0 to 1.5 orders of magnitude higher than the experimental results [45].

Figure 9c shows that the predicted D_e values of ¹³⁷Cs⁺ are consistent with the experimental results at a compacted density of 1400 kg/m³. However, a significant underestimation was observed at a compacted density of 800 kg/m³, with a difference of approximately four times. This can be explained by the limited number of experimental data points available for this density in the dataset, which comprised only 58 instances, accounting for approximately 7% of the total dataset. It indicates that additional diffusion experiments for ¹³⁷Cs⁺ should be conducted at a compact density of approximately 800 kg/m³ to facilitate the identification of diffusion patterns using ML models. Figure 9d illustrates that both the LGBM-CatB and LGBM-XGB models accurately predict the D_e of HTO. Under similar experimental conditions, the D_e in Wyoming bentonite (red squares) was higher than that in FEBEX bentonite (blue pentagrams), primarily because of the lower montmorillonite content, with m = 0.85 for Wyoming bentonite and m = 0.92 for FEBEX bentonite [54, 55].

Notably, the experimental diffusion data from this study, as well as from the ¹³⁷Cs⁺ [8] and HTO [54, 55] diffusions, were not included in the test datasets, highlighting the strong generalization ability of both LGBM-CatB and LGBM-XGB models. The generalization ability of LGBM-XGB was superior to that of LGBM-CatB, indicating that model selection plays a crucial role in accurately predicting radionuclide diffusion in complex geological environments. Given that HLW repositories have been designed to operate for over 10,000 years, the prediction of radionuclide diffusion in bentonite barriers must consider the complex coupling effects among radionuclides, porewater, and bentonite under intrinsic disposal conditions. Current diffusion datasets remain insufficient for safety assessments of bentonite barriers owing to limitations in data size and dimensionality. Therefore, additional diffusion experiments should be conducted to enhance the dimensionality and scale of the datasets.