1 Introduction

Compacted active clays are used for the construction of engineered barriers in underground nuclear waste repositories, where they are subjected to complex thermohydromechanical–chemical processes. This includes exposure to wetting from the saturated host rock and drying from the hot nuclear spent fuel [48], which makes the characterisation of the water retention behaviour of these materials particularly important. The water retention behaviour of clays is also relevant to other applications including, for example, the prediction of settlements under superficial foundations, the study of groundwater flow and the design of agricultural irrigation systems [3, 44, 69].

Early retention models have proposed a unique relationship between degree of saturation and suction, thus disregarding the effect of porosity on soil water content (e.g. [3, 15, 70]). Nevertheless, experiments published in the literature (e.g. [51,52,53, 62, 74]) have shown that changes in pore volume can significantly influence the water retention behaviour of the soil. In this sense, the effect of stress level and volume change and, in particular, the influence of the change in soil pore size distribution on the water retention of soils have also been studied [54, 68]. Subsequent models have, therefore, introduced an additional dependency of degree of saturation on void ratio without, however, considering the different structural levels of double porosity clays [16,17,18, 27, 37, 45, 57, 65, 67]. For instance, in the study developed by Huyghe et al. [22], the soil water retention equation proposed takes into account the soil volume deformation, establishing its relationship with the mechanical behaviour (i.e. the effective stress) from a mixture theory standpoint.

Two distinct porosity structures [19] have been experimentally observed in compacted clays, namely a microstructural porosity corresponding to the intra-aggregate voids and a macrostructural porosity corresponding to the inter-aggregate voids [6, 7, 52, 53]. Consistent with these observations, several models have been proposed to incorporate the effects of double porosity on the mechanical behaviour of compacted clays (e.g. [1, 20, 38, 42, 59]).

Only recently, some authors have, however, started to analyse the effects of porosity structure on soil–water retention and have proposed suitable modelling strategies that can account for the evolution of pore size distribution (e.g. [8,9,10,11, 23]). In this context, the present paper proposes a new approach to describe the influence of double porosity on soil–water retention. This approach combines the formulations of Navarro et al. [40] and De la Morena et al. [5] for separating microstructural and macrostructural water contents with a retention law that describes the effect of macrostructural volume changes on degree of saturation. This macrostructural retention law consists of an adaptation of a previous formulation proposed by Gallipoli et al. [17, 18] for single porosity soils.

The resulting model has been validated against a comprehensive set of experimental data from free swelling, constant load and confined (constant volume) wetting tests on three different active clays, i.e. MX-80 bentonite, FEBEX bentonite, and Boom clay. The main characteristics of these three clays and the relevant experimental procedures are summarised in the first part of the paper, followed by the description of the proposed retention model. Comparison between experiments and model simulations confirms that the incorporation of volumetric deformation inside a double porosity framework can improve the prediction of degree of saturation in compacted clays.

2 Materials and methods

Table 1 illustrates the main properties of MX-80 bentonite, FEBEX bentonite, and Boom clay. Retention tests on compacted samples of these three materials have been performed by a number of authors as listed in Tables 2, 3, and 4. In these tests, the samples were compacted and subsequently subjected to monotonic wetting with deionised water under free swelling, constant load or confined (constant volume) conditions. Since the initial water content was very low, the main wetting curve was assumed to be followed. The vapour equilibrium technique, consisting in the imposition of known levels of relative humidity to the sample, was generally employed to control suctions higher than 2–3 MPa (e.g. Romero [51]). Membrane cells or pressure plates) were usually employed to control suctions below 2–3 MPa. The samples analysed were not fully saturated in any case during the tests.

Table 1 Properties of MX-80 bentonite [28], FEBEX bentonite [14] and Boom clay [51, 53, 60]
Table 2 Free swelling tests on MX-80 bentonite, FEBEX bentonite, and Boom clay
Table 3 Constant load wetting tests on Boom clay
Table 4 Confined (constant volume) wetting tests on MX-80 bentonite, FEBEX bentonite, and Boom clay

Figure 1 shows the changes in water content recorded during free swelling and constant load wetting tests on MX-80 bentonite, FEBEX bentonite, and Boom clay (Tables 2 and 3). In these tests, volumetric strains were recorded throughout wetting, which allowed the measurement of void ratio in addition to water content. This is important for the calibration and validation of the proposed model, which, in contrast to the model previously proposed by Navarro et al. [40], includes the effect of the void ratio variation on the retention behaviour.

Fig. 1
figure 1

Free swelling and constant load wetting tests on a MX-80 bentonite, b FEBEX bentonite and c Boom clay

Figure 2 shows the changes in water content recorded during confined (constant volume) wetting tests on MX-80 bentonite, FEBEX bentonite, and Boom clay (Table 4). In these tests, wetting occurred at constant void ratio as the deformation of the specimen was prevented.

Fig. 2
figure 2

Confined (constant volume) wetting tests on a MX-80 bentonite, b FEBEX bentonite and c Boom clay

3 Model description

The proposed retention model is based on the double porosity concept by Gens and Alonso [19]. In this framework, as mentioned in Introduction, two structural levels are differentiated: the microstructure, which is associated with the voids inside the aggregates of clay particles, and the macrostructure, related to the space between the aggregates ( [30, 61]; among others). Different authors (e.g. [46, 52, 53, 74]) have suggested that, at high suction levels, virtually all water is adsorbed inside the micropores, and the presence of free water in the soil macropores is negligible. Therefore, over this suction level, defined as sMlim, the microstructural water content wm is assumed to coincide with the total water content w of the soil [52, 74]. For this reason, if the microstructural water is considered to occupy the whole microstructural pore space (see, for example, [1, 7, 77]), the microstructural void ratio em (microstructural pore volume per mineral volume) can be calculated over this suction range as:

$$ e_{{\text{m}}} = G_{{\text{S}}} \, w_{{\text{m}}} = G_{{\text{S}}} \, w $$
(1)

where GS is the specific gravity of the soil particles. Navarro et al. [40] and De la Morena et al. [5] applied this equation to obtain the values of em associated with different values of the thermodynamic swelling pressure π [35], which, in equilibrium between microstructural water and a low-salinity macrostructural aqueous solution can be defined as [42]:

$$ \pi = p + s_{{\text{M}}} $$
(2)

where p is the mean net stress, defined as the difference between the external mean stress pTOT and the pore gas pressure PG, and sM is the macrostructural matric suction, identified with the capillary suction: PG minus the liquid pressure PL. Navarro et al. [41] demonstrated that π controls the change in em, and therefore π can be understood as a microstructural effective stress. In addition, this expression agrees with the expression found by Borja and Choo [2] to define the microstructural effective stress. For monotonic hydration paths, em and π are related through a state surface as that defined by Navarro et al. [40] and De la Morena et al. [5]. The problem when characterising this surface arises when Eq. (1) is no longer valid, and the experimental values of w are associated with both wm and the macrostructural water content wM. This happens when sM is close to sMlim, and water starts to fill the macropores. The value of sMlim has been studied for different active clays. For example, Romero et al. [53] performed Mercury Intrusion Porosimetry (MIP) tests on Boom clay concluding that a suction of 2 MPa marks the onset of water flooding in the inter-aggregate pores. Similarly, Or and Tuller [46] indicated that the amount of free capillary water becomes negligible at suctions greater than 10 MPa for different soils. In another work on MX-80 bentonites, Jacinto et al. [25] found that dry density influences the retention behaviour only for suctions smaller than 30 MPa, which corresponds to the lower limit of the suction range dominated by adsorptive storage mechanisms [52]. In a similar work on MX-80 bentonite, Villar [74] found instead that water occupies only the soil micropores for suctions larger than 10 MPa. The results by Jacinto et al. [25] and Villar [74] indicate that, for MX-80 bentonites, the suction threshold marking flooding of macropores can vary between 10 and 30 MPa. In any case, these values are below 50 MPa, which, according to the psychrometric law [13], at 20 ºC, is associated with a relative humidity of 70%. This value is proposed by Cases et al. [4] and Woodruff and Revil [76] to differentiate crystalline hydration from hydration associated with double-layer processes. Therefore, in wetting processes, until sMlim is reached, there will be a nonnegligible range of water retention behaviour associated with double-layer processes. If, in that range of suctions, a logarithmic law is adopted to define it [12], the increase in microstructural void ratio, Δem, can be calculated as:

$$ \Delta e_{{\text{m}}} = - \kappa_{{\text{m}}} {\text{Ln}}\left( {\frac{{\pi + P_{{{\text{ATM}}}} }}{{\pi_{{{\text{ref}}}} + P_{{{\text{ATM}}}} }}} \right) $$
(3)

where πm can be understood as the microstructural stiffness and PATM is the atmospheric pressure. Any value of π can be considered as a reference πref, in particular, the 50 MPa mentioned above. Navarro et al. [40] and De la Morena et al. [5] verified the capacity of Eq. (3) to approximate the state surface emπ for π values lower than 50 MPa, which is applied in the present work. Experimental values associated with decreasing suction are progressively incorporated into the analysis, estimating the value of κm. When wM becomes relevant, Eq. (1) is no longer valid, and em cannot be estimated directly from w, since w = wm + wM. Therefore, κm can no longer directly relate π to water content, since an increment in water content in both em and eM occurs, and the quality of the estimation is reduced. Thus, sMlim is identified as the value of suction from which macrostructural water starts to be relevant.

For lower values of suction, wM and wm are differentiated assuming that the κm value identified for higher suction allows to estimate em (Eq. (3)) and, consequently, wm (Eq. (1)). In this way, the experimental water content data, which include both macrostructural and microstructural water content, are higher than the estimated wm values. Therefore, from the experimental values of w, wM is calculated as:

$$ w_{{\text{M}}} = w - w_{{\text{m}}} $$
(4)

Navarro et al. [40] and De la Morena et al. [5] showed the consistency of this procedure for MX-80 bentonite. However, it is interesting to assess its scope for a wider range of clays. This is done in the following section. In addition, as explained in Introduction, it is of great interest to include in the analysis the effect of the variation in the void ratio on the soil retention behaviour. For this purpose, the formulation proposed by Gallipoli et al. [17, 18] and Gallipoli [16] for a single porosity retention model (SPRM) is taken as a reference model:

$$ Sr_{{{\text{SPRM}}}} = \left[ {1 + \left( {s \, \frac{{e^{{{1 \mathord{\left/ {\vphantom {1 {\lambda_{{\text{s}}} }}} \right. \kern-\nulldelimiterspace} {\lambda_{{\text{s}}} }}}} }}{\beta }} \right)^{{{{\lambda_{{\text{s}}} } \mathord{\left/ {\vphantom {{\lambda_{{\text{s}}} } m}} \right. \kern-\nulldelimiterspace} m}}} } \right]^{ - m} $$
(5)

where SrSPRM is the degree of saturation computed with the model, s is the suction, e is the total void ratio (e = eM + em, being eM the macrostructural void ratio or macrostructural pore volume per mineral volume), and λs, β and m are material parameters. It is important to note that this model is thermodynamically consistent, as explained by Song and Borja [63]. If the model is expanded to a double porosity system, it seems reasonable to introduce the effect of the variation in the void ratio in SrM (macrostructural degree of saturation), and to do so using the macrostructural void ratio:

$$ Sr_{{\text{M}}} = \left[ {1 + \left( {s_{{\text{M}}} \, \frac{{e_{{\text{M}}}^{{{1 \mathord{\left/ {\vphantom {1 {\lambda_{{\text{s,M}}} }}} \right. \kern-\nulldelimiterspace} {\lambda_{{\text{s,M}}} }}}} }}{{\beta_{{\text{M}}} }}} \right)^{{{{\lambda_{{\text{s,M}}} } \mathord{\left/ {\vphantom {{\lambda_{{\text{s,M}}} } {m_{{\text{M}}} }}} \right. \kern-\nulldelimiterspace} {m_{{\text{M}}} }}}} } \right]^{{ - m_{{\text{M}}} }} $$
(6)

where λs, M, βM, and mM are macrostructural material parameters. Since SrM (i.e. wM) is relevant for nonnegligible suction (sMlim not lower than 2 MPa), a nonnegligible variation in the void ratio can be expected in unconfined conditions. This variation will be mainly due to eM, since, as shown by Navarro et al. [40] and De la Morena et al. [5] (and will be illustrated in the next section), the variation in em is not relevant for these suctions. Therefore, Eq. (6) represents a significant improvement with respect to the previously cited works. In addition, unlike other proposals (e.g. [8,9,10, 23]), the presence of air in the intra-aggregate is neglected (a hypothesis accepted and used by, e.g., [1, 7, 77]). Therefore, the microstructural strain model is directly related to the evolution of the microstructural water content, and additional parameters are not necessary to describe the microstructural water retention model.

Finally, the full double porosity retention model (DPRM) combines the above microstructural and macrostructural models to calculate the overall degree of saturation SrDPRM of the soil according to the following equation:

$$ Sr_{{{\text{DPRM}}}} = \frac{{Sr_{{\text{M}}} \, e_{M} + e_{{\text{m}}} }}{e} $$
(7)

4 Model calibration

This section presents the calibration of the proposed retention model against wetting tests on MX-80 bentonite, FEBEX bentonite, and Boom clay. The limiting value of macrostructural suction sMlim, which marks the onset of macrostructural flooding, was first determined from free swelling and confined wetting tests. To this end, the experimental values of water ratio ew (volume of water per volume of mineral) against the corresponding values of macrostructural suction sM for different tests at different dry densities have been plotted. These curves should overlap when sM is greater than sMlim since, in this case, all pore water exists inside the soil microstructure, whose retention behaviour depends only on the microstructural effective stress. The same curves should, however, start to diverge when sM becomes smaller than sMlim as water starts to flood the soil macrostructure, whose retention behaviour depends also on dry density due to the link between the macrostructural degree of saturation SrM and void ratio eM, see Eq. (6).

Figure 3a shows the experimental values of water ratio ew plotted against the macrostructural suction sM for free swelling and confined wetting tests at different dry densities on the three materials considered in the present study (Tables 2 and 4). Inspection of Fig. 3a indicates that, given the large scatter of data, it is not easy to identify the point where the experimental curves start to diverge. For that reason, the mean curve of each experiment was calculated using a centred moving average method [39], see Fig. 3b. To calculate the suction in which the mean curves start to diverge, a statistical analysis of the data in Fig. 3b was carried out. In this way, the values of sMlim were obtained as those for which the water ratios start to differ more than a 5%. This analysis leads to sMlim values equal to 20, 10 and 2 MPa for MX-80 bentonite, FEBEX bentonite, and Boom clay, respectively. These results are consistent with the values obtained in previous studies of the same materials (see Sect. 3).

Fig. 3
figure 3

a Experimental data from free swelling and confined wetting tests at different dry densities. b Mean curves of free swelling and confined wetting tests at different dry densities

The limits obtained from Fig. 3 can be transformed in terms of microstructural effective stress by means of Eq. (2). For free swelling tests, since the confinement is null during the tests, the mean net stress p can be assumed to be almost zero. Therefore, without salinity effects, the value of the macrostructural suction sMlim is equal to the microstructural effective stress πlim. For confined tests, the swelling pressure developed for the higher suctions (far from saturation) can be considered low and therefore negligible compared with sM. In this way, it can be assumed that πlim is roughly the same as sMlim in the three active clays studied. The values of πlim are reported in Table 5. It should be noted that, as expected, the value of πlim depends on the clay nature and, in particular, on the clay mineralogy [43, 50, 64], but also on the type of exchangeable cations [4, 32, 58].

Table 5 Values of πlim for the three clays considered in the present study

Once the value of πlim was determined, the state surfaces defined by Navarro et al. [40] and De la Morena et al. [5] were obtained using the results from free swelling tests in Table 2 limited to the range where π is greater than πlim. As stated in Sect. 3, for the range in which the hydration associated with double-layer processes is not negligible, the increase in microstructural void ratio is calculated using Eq. (3). As noted, in free swelling tests π coincides with sM. The values of πm and the obtained state surfaces are shown in Table 6 and in Fig. 4, respectively, for the three active clays studied. As expected, inspection of Fig. 4 indicates that, when the microstructural effective stress π is lower than πlim, the water content is higher than the value predicted by the microstructural model due to the presence of macrostructural water.

Table 6 Values of πm for MX-80 bentonite, FEBEX bentonite, and Boom clay
Fig. 4
figure 4

Fitting of microstructural model to data for a MX-80 bentonite, b FEBEX bentonite and c Boom clay. The dashed lines represent the corresponding values of πlim

Next, the parameter values of the macrostructural retention model were determined by fitting Eq. (6) to data from: (a) free swelling tests by Dieudonné et al. [11], Likos and Lu [31], Likos and Wayllace [32] and Saiyouri et al. [55] for MX-80 bentonite, (b) free swelling tests by Villar [52] for FEBEX bentonite, and (c) free swelling tests by Salager et al. [57] and constant load test (at constant vertical loads of 0.026, 0.55 and 0.6 MPa) by Romero [51] for Boom clay. Confined tests were not considered in this stage since the experimental values of the net mean stress developed during the tests were not available, and π could not be calculated. In each test, the value of em was computed with the microstructural model using the parameter values in Table 6. Moreover, Eq. (7) was employed to determine the experimental values of the macrostructural degree of saturation SrM. To simplify graphical representation of the data, Eq. (6) is rewritten in the following form [17]:

$$ Sr_{{\text{M}}} = \left[ {1 + \left( {\frac{{\overline{s}_{{\text{M}}} }}{{\beta_{{\text{M}}} }}} \right)^{{{{\lambda_{{\text{s,M}}} } \mathord{\left/ {\vphantom {{\lambda_{{\text{s,M}}} } {m_{{\text{M}}} }}} \right. \kern-\nulldelimiterspace} {m_{{\text{M}}} }}}} } \right]^{{ - m_{{\text{M}}} }} $$
(8)

where \(\stackrel{-}{s}\) M is the scaled macrostructural suction defined as:

$$ \overline{s}_{{\text{M}}} = s_{{\text{M }}} e_{{\text{M}}}^{{^{{1/\lambda_{{{\text{sM}}}} }} }} $$
(9)

In this way, since SrM is considered to depend on two independent variables (i.e. macrostructural suction and void ratio), the three-dimensional fit of Eq. (6) in the (SrM − sMeM) space is recast as a two-dimensional fit of Eq. (8) in the (SrM-\(\stackrel{-}{s}\) M) plane, see Fig. 5a. Therefore, the model can be represented as a 2D curve.

Fig. 5
figure 5

a Fitting of macrostructural retention model to experimental data of MX-80 bentonite, FEBEX bentonite, and Boom clay. b Fitting of single porosity retention model to experimental data of MX-80 bentonite, FEBEX bentonite, and Boom clay

For comparison, the single porosity retention model of Gallipoli et al. [17], Eq. (5), was also calibrated against the same experimental data, in Fig. 5b. Also in this case, the calibration is presented as a two-dimensional fit in the (SrSPRM-\(\stackrel{-}{s}\)) plane where \(\stackrel{-}{s}\) is the scaled suction defined according to Gallipoli et al. [17] as:

$$ \overline{s} = s \, e^{{1/\lambda_{{\text{s}}} }} $$
(10)

where e is obtained from the experimental data. The resulting parameter values are listed in Table 7 for both the macrostructural retention model and the single porosity retention model, which were obtained applying the least squared method.

Table 7 Parameter values of the macrostructural retention model and single porosity retention model for MX-80 bentonite, FEBEX bentonite, and Boom clay

To better visualise the results obtained, the comparison of both models with experimental data is also shown in the conventional Srs space in Fig. 6. In this case, 2D curves cannot be represented, since, for the same value of s, different values of e are possible, and vice versa. For Boom clay, since multiple tests were carried out with the same vertical load and initial dry density, only one test for each condition is represented for clarity.

Fig. 6
figure 6

Fitting of a double porosity model and b single porosity model to experimental data of MX-80 bentonite, FEBEX bentonite (tests from Villar [72]) and Boom clay in the Srs space

The accuracies of the double and single porosity retention models were initially compared by calculating the root-mean-squared errors (RMSE) of the respective predictions of degree of saturation in the previous calibration tests. For MX-80 bentonite, the RMSE of the double porosity retention model was two times lower than that of the single porosity retention model, i.e. 0.024 against 0.047. For FEBEX bentonite, the double porosity retention model yielded a RMSE of 0.015, while the single porosity retention model produced a RMSE of 0.024, i.e. a gain in accuracy of 1.6 times in the former model. Finally, for Boom clay, the RMSE was 0.019 for the double porosity retention model against 0.072 for the single porosity retention model, corresponding to a gain of accuracy of 3.7 times in the former model. The accuracy of both models is graphically shown in Fig. 7, which shows the higher correlation of the double porosity model with the experimental degree of saturation data.

Fig. 7
figure 7

Experimental versus predicted values of degree of saturation for a the proposed double porosity retention model and b the original single porosity retention model (calibration data)

5 Model validation

Although the results in the previous section are satisfactory, they do not demonstrate the improved predictive capabilities of the double porosity model, as they are based on the same experimental data employed during calibration of parameter values. For this reason, a validation exercise was carried out where additional data not used during calibration were employed for FEBEX bentonite and Boom clay. Therefore, the same parameters obtained in the calibration were used here. In particular, this additional validation used the free swelling tests on FEBEX bentonite by Lloret and Villar [33], Lloret et al. [34], Pintado [47] and Villar [73] and the constant load wetting tests (with vertical loads of 0.085 and 0.3 MPa) on Boom clay by Romero [51]. Figure 8 compares the predicted degree of saturation versus scaled suction for both the double and single porosity models against the corresponding experimental data.

Fig. 8
figure 8

Validation of a macrostructural retention model and b single porosity retention model against data of FEBEX bentonite and Boom clay

Figure 9 plots the same results obtained with the double and the single porosity models in the Srs space to make the comparison easier. Only one test for each vertical load is represented again for Boom clay. Inspection of Figs. 8 and 9 indicates that the double porosity retention model provides a more accurate prediction of the degree of saturation compared to the single porosity one. For the FEBEX bentonite, the double porosity retention model yields a RMSE of 0.016, while the single porosity retention model produces a RMSE of 0.027, i.e. a gain in accuracy of 1.7 times in the former model. The gain in accuracy associated with the use of the double porosity model is similar, even a little greater, to that observed for the calibration data. Similar results were obtained in the case of Boom clay with a RMSE of 0.016 for the double porosity model compared to 0.065 for the single porosity model, i.e. a gain in accuracy of 4.1 times in the former model. It is interesting to note that the model is particularly efficient for not very low suction values, where the sample is not completely saturated and the differences between both models are higher. These can be clearly seen for Boom clay. On the other hand, when the macrostructure is saturated, both models give similar results.

Fig. 9
figure 9

Validation of a double porosity retention model and b single porosity retention model against experimental data of MX-80 bentonite, FEBEX bentonite, and Boom clay in the Srs space

Figure 10 provides a complementary representation of the accuracy of the double and single porosity retention models by plotting their respective predictions of degree of saturation against the corresponding experimental values. This figure illustrates that the double porosity predictions lie significantly closer to the 1:1 line compared to single porosity predictions for both FEBEX bentonite and Boom clay. This is confirmed by the higher R2 values obtained.

Fig. 10
figure 10

Experimental versus predicted values of degree of saturation for a the proposed double porosity model and b the original single porosity model

6 Conclusions

This paper has presented a model to predict the retention behaviour of compacted clays with a double porosity structure consisting of micropores inside clay aggregates and macropores between clay aggregates. The model calculates the degree of saturation by explicitly considering the different retention mechanisms that occur at the microstructural and macrostructural scales. The approach relies on a new formulation that describes the effect of macrostructural volume changes on the variation in the degree of saturation of the soil. The formulation integrates the simulations of the distinct retention mechanisms within macropores and micropores while accounting for the effect of the volumetric deformation at separate porosity scales. The model requires the definition of three main components: (a) a swelling law for the saturated microstructure, (b) a criterion to assess the presence of water inside the macrostructure and (c) a retention law for the deformable unsaturated macrostructure. In particular, the retention behaviour of the deformable macrostructure was modelled by adapting a previously proposed law for predicting the variation in saturation in single porosity deformable soils.

A review of the existing literature was conducted to collect a considerable amount of data from wetting tests on compacted MX-80 bentonite, FEBEX bentonite, and Boom clay under free swelling, constant load and confined conditions. The collected data were used to formulate, calibrate and validate the proposed double porosity retention model. In particular, part of the experimental data were used to calibrate the model parameters, while the rest of the data were employed to validate the calibrated model. Predictions from the proposed double porosity retention model were also compared with those of an existing single porosity retention model for deformable soils. This comparison confirmed that the explicit consideration of volumetric deformations at different structural scales significantly improves the accuracy of the prediction of degree of saturation in compacted clays.