A triple porosity hydro-mechanical model for MX-80 bentonite pellet mixtures

Journal Pre-proof A triple porosity hydro-mechanical model for MX-80 bentonite pellet mixtures Vicente Navarro, Laura Asensio, Heidar Gharbieh, Gema De la Morena, Veli-Matti Pulkkanen

PII:

S0013-7952(19)30362-X

DOI:

https://doi.org/10.1016/j.enggeo.2019.105311

Reference:

ENGEO 105311

To appear in: Received Date:

25 February 2019

Revised Date:

10 July 2019

Accepted Date:

24 September 2019

Please cite this article as: Navarro V, Asensio L, Gharbieh H, De la Morena G, Pulkkanen V-Matti, A triple porosity hydro-mechanical model for MX-80 bentonite pellet mixtures, Engineering Geology (2019), doi: https://doi.org/10.1016/j.enggeo.2019.105311

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier.

A triple porosity hydro-mechanical model for MX-80 bentonite pellet mixtures

Vicente Navarro a*, Laura Asensio a, Heidar Gharbieh b, Gema De la Morena a, Veli-Matti Pulkkanen b a

Geoenvironmental Group, Universidad de Castilla-La Mancha. Avda. Camilo José Cela

s/n, 13071 Ciudad Real, Spain.

of

VTT Technical Research Centre of Finland Ltd. P.O. Box 1000, FI-02044 VTT, Finland.

ro

b

*: Corresponding author: Vicente Navarro. Tel.: +34 926295300. Email:

re

-p

[email protected].

E-mail addresses: [email protected] (V. Navarro), [email protected] (L.

lP

Asensio), [email protected] (H. Gharbieh), [email protected] (G. De la

Jo

ur na

Morena), [email protected] (V.M. Pulkkanen)

Abbreviations BBM Barcelona Basic Model BExM Barcelona Expansive Model CM Comsol Multiphysics DPM Double Porosity Model PTPM Pellets Triple Porosity Model 1

Highlights A triple porosity model for MX-80 bentonite pellets mixtures is proposed (PTPM) An overlapping continua approach is adopted to describe the three functional levels The model for compacted bentonite blocks is adapted to m and M1 The PTPM improves the time evolution prediction of double porosity models

of

The PTPM improves the simulation of the behaviour of open pellet mixtures

ro

Abstract

This paper describes the application of a triple porosity model to characterise the hydro-

-p

mechanical behaviour of MX-80 bentonite pellet mixtures. The macroscopic behaviour of the system was assumed to be fundamentally controlled by three functional levels, the

re

microstructure (associated with the pores existing inside the clay particle aggregates), the

lP

Macrostructure 1 (which integrates the space between the aggregates inside the pellets) and the Macrostructure 2 (related to the void space between the pellets). The role played by each functional level in the overall behaviour of the system was modelled by means of an

ur na

overlapping continua approach. The usual framework in the simulation of the behaviour of compacted bentonite blocks has been adapted to describe the behaviour of the microstructure and the Macrostructure 1. To model the behaviour of the Macrostructure 2, a

Jo

simplified formulation of granular materials was adopted. To illustrate the capacity of the proposed approach, two hydration tests were simulated in which significantly different pellet mixtures were used: an open mixture with a dry density of 953 kg/m3, and a closed mixture with 1520 kg/m3 of dry density. Despite the difference between the two mixtures, the model allowed very satisfactory fits to be obtained using in

2

both cases the same material parameters. This shows the consistency of the proposed conceptual framework. Since the response of the two Macrostructural levels is differentiated, the model provides a richer description of the mixture behaviour than that obtained using a double porosity model. However, the interest of the model lies not only in its descriptive richness. The differentiation of the effect of both Macrostructural levels improves the description of the

of

time evolution of the system compared with that obtained when using a double porosity

ro

model. However, if the Macrostructure 2 is small (as occurs in the second of the two tests analysed), or there is sufficient confinement for the inter-pellet void space to be

-p

considerably reduced when the microstructure swells, double porosity models can

re

satisfactorily characterise the performance of the bentonite pellet mixture.

lP

Keywords

Bentonite pellet mixtures, bentonite, hydro-mechanical model, pellets triple porosity model,

Jo

ur na

double porosity model, high level radioactive waste

3

1. Introduction Because of their low hydraulic conductivity, swelling and chemical retention capacity, bentonites are popular for components of the engineered barrier systems of different repository concepts for high level radioactive waste. These components can include bentonite in the forms of compacted blocks, pellets, granules, powder, or a combination of them (Sellin and Leupin, 2013). An intensive research effort has been made to characterise

of

the thermo-hydro-chemical-mechanic behaviour of bentonite to better assess its

ro

performance in the repositories. This is especially true for compacted bentonite blocks,

obtaining very satisfactory results when applying conceptual models based on a double

-p

porosity approach. In this regard, Guimarães et al. (2013), Nowamooz et al. (2009),

re

Sánchez et al. (2005) and Vilarrasa et al. (2016) present examples of the capacity of the Barcelona Expansive Model (BExM) (Gens and Alonso, 1992; Alonso et al., 1999) to

lP

reproduce the main characteristics of the behaviour of compacted bentonite blocks. BExM, probably the most widely used double porosity model (DPM) nowadays, is not based on

ur na

homogenization techniques (for more details on those techniques, see, for instance, Borja and Koliji, 2009; Mainka et al., 2014), but on the overlapping continua approach introduced by Barrenblatt et al. (1960) for fissured rock, extended to any soil by Khalili et al. (1999) and Wilson and Aifantis (1982), among others. These models assume that two different

Jo

continuous media overlap occupying the same spatial domain (Hueckel, 1992). For bentonites, the differentiation between the two media is usually not based on the microscopic topology of the material but on a functional criterion: the possibility of water flow occurring as a consequence of hydrodynamic gradients. One of the continua, commonly identified as “microstructure”, models the behaviour of a continuum equivalent to the porous medium formed by the pores in which such flow is negligible. The 4

“Macrostructure” models the continuum equivalent to the rest of the porous medium. The “restrictive” definition of the microstructural functional level (no hydrodynamic flow: adsorbed water) implies that, despite its functional definition, it has a specific structural support, since it is fundamentally associated with the aggregates, that is, the intra-aggregate void space. This provides consistency to the macroscopic magnitudes of the microstructure, because they characterise the state of a part of the domain that is expected to develop the

of

same behaviour trend. This is not the case with the Macrostructure, since its “wider”

ro

definition (“what is not microstructure”) means that it encompasses a relevant range of different pore sizes, see Fig.1a, in which different behaviour trends may occur.

-p

Nevertheless, the good results obtained when using DPMs do not provide evidence to

re

support this idea (Guimarães et al., 2013; Navarro et al., 2017; Nowamooz et al., 2009; Sánchez et al., 2005; Vilarrasa et al., 2016). On the contrary, they seem to indicate that,

lP

despite the differences, the same constitutive model can be used for the pores considered, which provides consistency to the hypothesis of an equivalent single continuum.

ur na

However, when bentonites with a multi-modal pore structure that have much larger pores are considered, such as pellet mixtures (Fig. 1b), the hypothesis of two functional levels is less clear. Some authors have maintained the same concept as in blocks, merging the interaggregate spaces and the inter-pellet spaces in the Macrostructural functional level (Gens et

Jo

al., 2011). But this, especially before the pellets swell and reduce the inter-pellet space, means integrating two very different behaviour trends into the Macrostructure. Regarding flow, this idea would imply that the flow in a gravel (inter-pellet) is considered to be similar to the flow in a compacted clay (intra-pellet). In others works, the Macrostructure has been identified with the larger inter-pellet pores and all the intra-pellet pore spaces are merged in the microstructural functional level (Alonso et al., 2011). In this case, the 5

microstructure comprises the Macrostructure and microstructure considered in compacted bentonite blocks, with remarkably different behaviours. Therefore, to adjust the considered functional levels to the expected behaviour trends, while maintaining the overlapping continua approach, this paper proposes to introduce a new functional level to characterise the behaviour of the inter-pellet porosity, adopting a triple porosity model to describe the macroscopic behaviour of pellet mixtures (PTPM: Pellets Triple Porosity Model). The

of

modelling framework is similar to the triple porosity models presented by Straughan (2017)

ro

and Svanadze (2018) but adapted for bentonites.

In the following sections, the constitutive basis of this proposal is described. The flow

-p

model, the stress-strain relationships proposed for each functional level, and the equations

re

that characterise the interaction (mass exchange and induced strains) between the three overlapping continua are defined. Afterwards, the simulation of two hydration tests is

lP

presented to demonstrate the scope of the proposed formulation. Before starting the description of the model, the next Section describes the material used and the procedures

Jo

ur na

followed to perform the modelled tests.

6

of

ro

Fig. 1. Pore size distributions obtained with mercury intrusion porosimetry of (a) a

compacted bentonite sample (dry density 1800 kg/m3; adapted from Lloret et al., 2003), and

-p

(b) a bentonite pellets sample (prepared using a gravity fall compaction, overall dry density

re

of 1150 kg/m3; adapted from Hoffmann et al. 2007). (PDF: Pore Distribution Frequency).

lP

2. Materials and Methods

The tests simulated in Section 4 were carried out by Clay Technology (Sweden) (Åkesson

ur na

et al., in prep.) and CEA (France) (Bernachy-Barbe et al., in prep.). The material used was MX-80 bentonite, a sodium bentonite from Wyoming (USA) with a montmorillonite content of 80% (Andersson and Sandén, 2012; Molinero Guerra et al., 2017). Although the material was labelled as MX-80 in both cases, it came from two different suppliers, and the

Jo

reported properties are not exactly equal. In the first test, the bentonite was supplied by HOSOKAWA Bepex GmbH. Its cation exchange capacity is of 75 meq/100 g, the main exchangeable cations are Na (55 meq/100 g) and Ca (12 meq/100 g), and the mineral density is 2780 kg/m3 (Karnland et al., 2006). In the second test, the bentonite was supplied by Laviosa-MPC. Its cation exchange capacity is of 98 meq/100 g, the main exchangeable

7

cations are Na (52 meq/100 g) and Ca (37 meq/100 g), and the mineral density is 2770 kg/m3 (Molinero Guerra et al., 2017). In the first test, the pellets had a pillow shape with a maximum dimension of 16 mm (Luterkort et al., 2017). They were produced by compaction between two counter-rotating rollers to a pellet dry density of 1841 kg/m3 and a water content of 15%. In the second test, the pellets had a cylindrical central part and two spherical poles, and their maximum

of

dimension was 32 mm (Bosgiraud and Foin, 2016). They were produced by compacting

ro

bentonite powder in a mould to a pellet dry density between 2010 and 2050 kg/m3 and a

water content of 4%. This test also contained crushed pellets (30% in mass) with a particle

-p

size lower than 2.5 mm.

re

In both cases, the sample was contained in a cylindrical cell. The cell in the first test had an inner diameter of 100 mm and a height of 500 mm. The cell was filled with a pellet mixture

lP

to a dry density of 953 kg/m3. In the second test, the cell had an inner diameter of 240 mm and a height of 105.15 mm. The cell was initially filled with a 70% in mass of pellets and a

ur na

30% of crushed pellets, poured alternatively in three layers each plus a final layer of split pellets. The initial dry density of the sample was 1520 kg/m3. In both cases, the cell was closed before the start of the test. Water was then supplied to the bottom face of the sample through a filter, pressurized to 1 kPa and 10 kPa above atmospheric pressure in the first and

Jo

second test, respectively. The tests were run in such conditions for approximately 200 days and 940 days, respectively.

3. Conceptual model 3.1 System abstraction

8

In this first application of the PTPM, a hydro-mechanical analysis is carried out, assuming isothermal conditions and a sufficiently reduced salinity to neglect its effect on the system's behaviour. In addition, the pressure of gas (which is assumed to consist of only two species: water vapour and “dry” air, or simply air, both perfect gases) is considered constant and equal to the atmospheric pressure, Patm, in the whole system. Therefore, only the liquid water flow and the water vapour transport are taken into account. Both processes are

of

considered to occur both in M1 (functional level associated with the void space that exists

ro

between the clay particle aggregates inside the pellets) and in M2 (functional level

associated with the inter-pellet void space). On the contrary, the water mass flow with

-p

respect to the bentonite skeleton in m (functional level associated with the microstructural

re

pores) is assumed null, since, as implicit in the functional definition of the microstructure, water is supposed to be linked to the bentonite skeleton. As usual (see for example Gens

lP

and Alonso, 1992; Mašín and Khalili, 2016; Yong, 1999), the microstructural pores are assumed saturated (therefore, without vapour).

ur na

The model takes into account the potential exchange of water mass between the three functional levels considered. However, although the three media simultaneously occupy the same spatial domain, the void space describing m is assumed to be surrounded by the void space of M1 (as stated by Alonso et al., 2011), and M1 is in turn surrounded by the void

Jo

space of M2. Therefore, a hierarchy in the interaction is considered, assuming that m cannot exchange mass with M2 directly. Any mass exchange must occur through M1. It is assumed that the rearrangement of m may cause strains in M1, but not in M2. In the same way as when modelling compacted bentonite blocks it is accepted that the rearrangement of the Macrostructure does not cause the rearrangement of the microstructure (Lloret and Villar, 2007), this PTPM assumes that the rearrangement of M1 9

does not cause strains in m, and that the rearrangement of M2 does not cause strains in M1 either.

3.2 Mathematical model The state of the system is therefore defined by four state variables: u, vector of the solid skeleton displacements; em, microstructural void ratio (volume of intra-aggregate voids per

of

volume of mineral), which defines the packing of the microstructure; PL1, liquid pressure in

ro

the functional level M1, which determines the flow conditions and the water content in M1; and PL2, which will do the same in M2. Their values are obtained by solving in a coupled

-p

way three mass balance equations (balance of water mass in m, M1 and M2) to obtain the

re

flow state variables (respectively, em, PL1 and PL2), and the mechanical equilibrium to

lP

obtain u (vectors and tensors are identified in bold throughout the article).

3.2.1 Water flow model. State functions, field and balance equations

[1]

ur na

The mass balance equations have the same structure as in DPM:

mk    l k  rk t

where “ /t” and “·” define the partial time derivative and the divergence operator,

Jo

respectively, and mk is the volumetric density of the water mass in each functional level (k: m, M1 and M2):

em 1 e

[2]

mm  m

[3]

mMi,L  mMi,V   W

eMi SrMi 1 e

 V,Mi

eMi (1  SrMi ) 1 e

10

; i  1, 2

where m and W are, respectively, the densities of micro and free water, SrMi (i=1, 2) defines the degree of saturation of Mi, eMi is its void ratio, and e is the total void ratio (e=eM1eM2em). In Eq [1], lk defines the water mass flow: [4]

l m  mm v

[5]

l Mi  mMi v  WqMi  jMi ; i  1, 2

of

where v is the solid skeleton velocity vector (time derivative of u), qMi is the liquid water specific discharge, and jMi is the vapour diffusion. In turn, rk represents a mass source term.

ro

If there are no external sources or sinks, this term is only associated with the mass

rk 



j m,M1 ,M 2

rk  j

re

[6]

-p

exchange between the functional levels:

where rkj defines the mass contribution to the functional level k from the functional level j.

lP

According to the description of the system abstraction (Section 3.1), only the exchanges rmM1 (=  rM1m) and rM1M2 (=  rM2M1) are not zero. The exchange rmM1 is modelled

[7]

ur na

using the non-linear formulation proposed by Navarro et al. (2015):

 p  sM1  rm-M1   W H m-M1     

Cm-M1

   p  s   M1

where the material parameters Hm-M1= 1.5×10−8 (kPa·s)−1 and Cm-M1 = 0.4 are assumed

Jo

(Table 1). It should be noted that these parameters were identified by modelling a double porosity system. Therefore, rmM1 defines mass exchange from M1 to m per unit volume if the porosity associated with M2 is not present. Consequently, the correction described in Appendix A should be applied. To calculate rM2-M1, the approach proposed by Huyakorn et al. (1983) for a fractured rock is used:

11

[8]

rM2-M1  W H M2-M1  PL1  PL2 

where HM2-M1=1.0×10−11 (kPa·s)−1, adapted from Gens et al. (2011), is adopted (Table 1).

As in DPM (see for instance Olivella et al., 1994, and Rutqvist et al., 2001), and as Svanadze (2018) did when applying a triple porosity model to analyse the thermoelastic

[9]

q Mi  

kI,Mi kr,Mi

W

( PLi +  W g z ) ; i  1, 2

ro

discharge qMi using the advective approach by Pollock (1986):

of

behaviour of a porous medium, it is proposed to calculate the liquid water specific

-p

where g is the gravitational acceleration, z is the vertical coordinate upwards oriented, “”

viscosity (Ewen and Thomas, 1989):

W  661.2 103 T  229 

1.562

lP

[10]

re

expresses the gradient differential operator, and W defines the liquid water dynamic

where T is the absolute temperature and W is obtained in Pa·s.

ur na

The model introduced by Gens et al. (2011) is adopted to compute the liquid intrinsic permeability kI,M1 of M1 as a function of its porosity M1: [11]

kI,M1  kO1 exp  b1 M1  O1  

Jo

As occurs when rm-M1 is calculated using the parameters deduced from a DPM, if the material parameters kO1, b1 and O1 (identified when the flow in compacted bentonite blocks is simulated, see Table 1) are used to calculate kI,M1, it should be noted that the specific discharge is obtained without taking into account the presence of the inter-pellet void space. Therefore, the correction outlined in Appendix A should be also applied. The liquid relative permeability kr,M1 can be modelled using the Brooks and Corey (1964) 12

formulation, assuming an exponent value of 3 for the Macrostructural degree of saturation for water SrM1 (volume of water in M1 per volume of voids in M1), as it is done in DPMs (Gens et al., 2011): [12]

kr,M1   SrM1 

3

A van Genuchten (1980) retention curve is used to calculate SrM1:

SrM 1 

1  

1

1 sM1  1 n



m1

of

[13]

ro

where the Macrostructural suction is defined as sM1 = PG – PL1, PG being the gas pressure,

-p

equal to Patm in the present paper. The relationship n1 = 1/(1m1) is assumed, and the approach of Navarro et al. (2015) is applied to obtain differentiated models of the intra-

re

aggregate and inter-aggregate water content. Table 1 includes the material parameters 1

lP

and m1 (identified for a main wetting path), as well as the other hydraulic parameters of M1 used to perform the simulations presented in Section 4. The intrinsic permeability of M2, kI,M2, was deduced from the saturated hydraulic

ur na

conductivity proposed by Ergun (1952) for granular solids based on the classic KozenyCarman model:

[14]

kI,M2

2 e  4  Req  M2  150 1 e

3

Jo

In this formulation, an equivalent radius Req is assigned to the pellets, which can be derived from the evolution of the inter-pellet void ratio, em+eM1: [15]

Req  Req,ini

3

1  em  eM1 1  em,ini  eM1,ini

13

Both the relative permeability, kr,M2, and the degree of saturation, SrM2, of M2 are defined according to the formulation of Brooks and Corey (1954): [16]

[17]

kr,M2   SrM2 

SrM2

23 



 if sM2  sA , SrM2  1     sA   otherwise, SrM2    sM2  

of

Consequently, the air entry pressure sA, the shape parameter  and Req,ini are the flow

ro

parameters of M2. Req,ini depends on the geometry of the pellets arranged. Values of 0.5 kPa

-p

and 1.1 (Table 1) are taken for sA and  respectively, consistently with the values for NW5 gravel, Pea gravel and Sublayer (1.0-2.0 mm) in Benson et al. (2014). Like sM1, suction in

re

M2, sM2, is also identified with the capillary suction: sM2 = PG – PL2.

As in DPM, the vapour flow in M1 and M2 is calculated adopting an advective-diffusive

lP

model. Since the gas pressure is assumed constant, the advective flow of the gas phase can be considered null, and the hydrodynamic dispersion can be neglected. Therefore, the

[18]

ur na

vapour diffusion jMi (i: 1 and 2) is modelled as the vapour molecular diffusion:

jMi    Mi Mi (1  SrMi ) DV  V, Mi  ; i  1, 2

The binary diffusion coefficient of water vapour in gas DV is calculated as proposed by

Jo

Pollock (1986) and Philip and De Vries (1957): [19]

DV  5.9 10

3

T 2.3 Patm

where Patm should be introduced in kPa to obtain DV in m2/s. Tentatively, as Olivella and Gens (2000) did when simulating a double porosity soil, the soil tortuosity to vapour flow

14

Mi is assumed equal to 1 both in M1 and M2. To determine the density of water vapour V,Mi, the psychrometric law is used (Edlefsen and Anderson, 1943): [20]

 WMM  sMi  ; i  1, 2  W R T 

V,Mi  VO exp  

where WMM is the molar mass of water and R is the universal gas constant. The formulation in Ewen and Thomas (1989) is applied to calculate the density of saturated

VO 

exp  0.06374 (T  273.15)  0.1634 103 (T  273.15)2 

ro

[21]

of

water vapour VO (in kg/m3) as a function of the absolute temperature T:

194.4

-p

The density of free water W is assumed to be constant and equal to 1000 kg/m3 for the

re

isothermal problems considered. In addition, although the density of adsorbed water m may be greater than W (see Jacinto et al., 2012), m  W is assumed for modelling

lP

purposes (Tournassat and Appelo, 2011).

To conclude the description of the flow problem analysis, it is important to note that the

ur na

value of the void ratios will be updated as the simulation progresses over time. The value of em is automatically updated, since, as introduced at the beginning of this Section, it is the state variable associated with solving the equation of the water mass balance of the microstructure. Therefore, as it is done in some DPM (see, for example, Navarro et al.,

Jo

2017), its updating is inherent to the simulation process. On the other hand, if there are not erosion processes or external contribution of mineral mass, the total void ratio e (total volume of voids per volume of mineral) can be obtained explicitly considering the mineral mass balance equation (Navarro et al., 2016): [22]

mmineral     mineral v   0 t 15

where mmineral (= mineral/[1+e]) is the mineral mass per unit volume, and mineral is the mineral density. Using the sign convention in Soil Mechanics, the volumetric strain ratio dV/dt can be expressed as: [23]

d V    v dt

Therefore, from the above two equations it is deduced that:

e  1  eO  exp   V 

of

[24]

ro

where eO is the initial value of the total porosity e, and V is obtained by the time

integration of the increment of total volumetric strain dV, computed from the total strain

d V  mt  dε

re

[25]

-p

increment vector from d (defined in the next Section) as:

transpose operator.

lP

where m is the vector form of the Kronecker delta, and the superscript “t” indicates the

If M2 did not exist, the value of eM1 could be determined from e and em, since e = em + eM1

ur na

would be met. This is done in DPM. However, it is not true in PTPM, since the relationship that holds is e = em + eM1 + eM2. Therefore, even knowing em and e, it is necessary to calculate the value of the void ratio in one of the two Macrostructural levels to obtain the void ratio in the other. In this work, eM2 is calculated by the time integration of its

Jo

increment deM2, obtained from dV,M2, i.e., the volumetric strain increment associated with M2:

[26]

deM2   1  e  d V,M2

where dV,M2 is given by dM2 (strain increment associated with M2), which, as d, is defined in the next Section. 16

Table 1. Hydraulic and mass-transfer parameters. Definition

Value

b1

Material parameter of the intrinsic permeability kI,M1 9.91 a

Cm-M1

Constitutive parameter of the mass exchange rm-M1

0.4 b

Hm-M1, 1/(kPa·s)

Constitutive parameter of the mass exchange rm-M1

1.5×10−8 b 1.0×10−11 c

kO1, m2

Reference intrinsic permeability of M1

2.34×10-21 a

m1

Parameter of the water retention curve of M1

0.733 a

sA, kPa

Parameter of the water retention curve of M2

0.5 d

α1, 1/kPa

Parameter of the water retention curve of M1

1.15×10−4 a



Parameter of the water retention curve of M2

1.1 d

O1

Reference porosity of M1 associated with kO1

0.047 a

re

-p

ro

of

HM2-M1, 1/(kPa·s) Constitutive parameter of the mass exchange rM2-M1

lP

a

Parameter

Data from Navarro et al. (2017). b Data from Navarro et al. (2015). c Data adapted from

ur na

Gens et al. (2011). d Data adapted from Benson et al. (2014).

3.2.2 Mechanical problem. Equilibrium equation and constitutive model

Jo

Even if no account is taken of cracking processes that wetting can cause in pellets, in addition to the increment in strain caused by the rearrangement of the functional levels m (dm) and M1 (dM1), both considered in DPM, the increment in strain due to the rearrangement of M2 (dM2) should also be considered in the PTPM proposed here: [27]

dε  dε m  dε M1  dε M2

17

The same conceptual approach used when the behaviour of compacted bentonite blocks is simulated using a DPM can be used to calculate dm and dM1. Therefore, dM1 can be computed as: [28]

dε M1  dε M1 M1  dεM1 m

where dM1M1 defines the increment of strain that the rearrangement of M1 causes if the

of

volume of the microstructure remains constant, and dM1m defines the rearrangement caused by dm in M1. For the isothermal conditions assumed in this work, dM1M1 can be

ro

modelled using the Barcelona Basic Model (BBM; Alonso et al., 1990). The model is

-p

formulated in terms of two significant stresses, net stress and Macrostructural suction. The net stress is defined as =   PG·m, where  is the total stress vector (the

re

engineering or Voigt notation is used for both stress and strain tensors throughout the paper). The BBM has become the reference model in modelling unsaturated soils (see

lP

Sheng et al., 2004; Sołowski et al., 2012). In this work, it is assumed that plastic strains only occur if the ellipsoidal generalisation of the load-collapse (LC) yield surface (Alonso

[29]

ur na

et al., 1990) is achieved. Therefore, dM1M1 can be expressed as:

dεM1M1  dεe M1 σ  dεe M1 s  dε p M1 LC

where deM1- defines the elastic strain that occurs when changing  in M1, deM1s

Jo

introduces the increase in elastic strain caused by the changes of sM1 in M1, and dpM1LC describes the plastic strain that occurs when the LC yield surface is reached. According to Navarro et al. (2017), the coupling term dM1m is caused by two differentiated mechanism, which in volumetric terms can be expressed as: [30]

d V, M1 m  d

 d V, M1 m, 2

V, M1 m,1

18

The first mechanism (first term on the right) is a consequence of changes in packing that the Macrostructure undergoes when the volume of the aggregates changes. It can be described using the BExM by means of two interaction functions: fC, when em decreases (microstructural compression), and fS when em increases (microstructural swelling). The shape of both functions is indicated in Fig. 2a. The value of the functions depends on the degree of openness of the Macrostructure relative to the applied stress state, expressed as

of

pR/pO (Sánchez et al., 2005), where p is the net mean stress (p = pTOT − PG, pTOT being the

ro

mean stress) and pR is defined in Fig. 2b. In normal consolidated conditions, pR is equal to pO, the net mean yield stress at the current suction (see Fig. 2b). Once fi (i= C or S) is

d V,M1 m,1  fi d V, m

re

[31]

-p

defined, dV,M1m,1 can be obtained as:

[32]

d  V,m  

dem 1 e

lP

where dV,m, the volumetric strain associated with dm, can be calculated as:

ur na

The value of em, and consequently of dem, is defined by the state surface of Fig. 3. In this figure, based on the analysis of the change in microstructural water content (or microstructural void ratio) for a large number of MX-80 samples (Navarro et al., 2015),  defines the thermodynamic swelling pressure (Navarro et al., 2018), and emR is the

Jo

microstructural void ratio remaining in dry conditions.

19

of ro -p re

Fig. 2. (a) Shape of the BExM interaction functions fC and fS. (b) Yield surface of the BBM

Jo

ur na

from Sánchez et al. (2005).

lP

in the p-q space and definition of the reference stress pR (CSL: critical state line). Adapted

20

Fig. 3. State surface that relates the microstructural void ratio em with the thermodynamic swelling pressure  Adapted from Navarro et al. (2015).

The second term on the right of Eq. 30 describes the free swelling that the Macrostructure experiences when the destructuration of the microstructure occurs under low confinement

of

conditions. Saiyouri et al. (2004) suggested that free swelling is a process more closely linked to the subdivision of mineral clay particles than a homogeneous increase in the

ro

distance between layers. Aggregates undergo similar subdivision processes (Salles et al., 2009). This subdivision inside the aggregates generates a new porosity M1 introduced in

-p

the model by dV,M1-m,2. It is proposed to model this term using the formulation presented

re

by Navarro et al. (2017), which allowed to obtain satisfactory results when simulating free swelling processes under different salinity conditions.

lP

Once defined how to obtain dM1 and dm, the way in which dM2 is calculated should be determined to complete the constitutive formulation of the mechanical behaviour of the

ur na

PTPM. For simplicity, in this first approach to the use of three functional levels, a nonlinear isotropic elastic model is adopted, such as the one proposed by Lade and Nelson (1987) to describe granular materials. The formulation is adapted to partial saturation conditions assuming in M2 the definition of effective stress ’ proposed by Alonso et al. (2010), ’=

Jo

 SrM2·sM2·m, where, as stated by Alonso et al. (2010), the term multiplying sM2 (i.e., SrM2), defines all the water available in the modelled media, M2. For the axisymmetric conditions in the tests considered in the present paper (Section 2), dM2 is calculated as:

21

 d  V,M2σ   d  V,M2s  dε M2  dε M2σ  dε M2s      d  S,M2σ   0  [33]

dε M2σ

1 K M2 e  C 2  dσ    0

dSrM2  SrM2  dsM2    dp  ; dε dsM2 V,M2 s     1 dq  K M2  3GM2  0

where q is the von Misses stress, S,M2 is the deviatoric strain associated with M2, and Ce2 is

of

the elastic compliance matrix associated with M2. The Poisson’s ratioM2 is assumed to be constant, considered in a first approximation equal to that of M1, M1, Table 2. The shear

GM2 

3 1  2 M2  K M2 2 1  M2 

-p

[34]

ro

modulus GM2 is obtained as a function of the bulk modulus KM2 as:

[35]

K M2 

p

 M2

lP

porosity clayfills is proposed to be used:

re

To determine the bulk modulus, the equation introduced by Mašín et al. (2005) for double

[36]

ur na

where the following expression is used to calculate the elastic stiffness M2:

 M2  cM2 exp  dM2 eM2 

adapted from Najser et al. (2012). Both cM2 and dM2 are material parameters. In this first approximation, cM2 = p,M1/10 and dM2 = Ln(3×p,M1/cM2)/eM2,ini are taken. Thus, for initial

Jo

conditions (eM2 = eM2,ini) M2 will be three times p,M1 (elastic stiffness of M1 for changes in net stress; Table 2), and when eM2 becomes very low, since M2 tends to cM2, it will be an order of magnitude lower than p,M1. The formulation introduced (Eqs. from [28] to [36]) allows to write Eq. [27] as: [37]

dε  dε m  dεe M1 σ  dεe M1 s  dε p M1 LC  dε M1 m  dε M2 σ  dε M2 s 22

where it follows that: [38]

dεe M1 σ  dεM2 σ   Ce1  Ce 2   dσ  dε   dε m  dεe M1 s  dε p M1 LC  dε M1 m  dε M2 s 

Ce1 being the compliance matrix associated with M1, defined analogous to Ce2, but taking the bulk modulus Kp,M1 of M1 as:

K p,M1 

1  e  p  p,M1

of

[39]

The shear modulus GM1 is calculated as a function of Kp,M1 in the same way as GM2 is

ro

calculated as a function of KM2 (Eq.[34]). If the elastic matrix De of the whole system is defined as: De   Ce1  Ce 2 

1

-p

[40]

re

(where the subscript “1” indicates the inverse matrix), the increment of net stress d, equal in M1 and M2 since the value of PG is the same in both functional levels, is given by:



dσ  De  dε   dε m  dε e M1 s  dε p M1 LC  dε M1 m  dε M2 s 

lP

[41]



ur na

The use of a unique elastic matrix for M1 and M2 allows to apply a formulation analogous to that used in DPM (see, for instance, Sánchez et al., 2005) to calculate dpM1LC at each computational time step. For this purpose, it is enough to take into account the additional effect of dM2s. Thus, using Eq. [41], d is obtained, and consequently TOT, as a function

Jo

of dem (which allows to compute dm), dPL1, dPL2 and du. The mechanical problem will be solved in that time step if, besides meeting the boundary conditions of displacement and load, the mechanical equilibrium equation is met: [42]

 σTOT   g z  0

where  is the soil bulk density.

23

Table 2. Macrostructural stress-strain parameters. Parameter Definition Increase in cohesion with suction

0.1 a

pC, kPa

Reference stress of the net mean yield stress (BBM in M1)

10 a

r

Macrostructural soil compressibility parameter (BBM in M1)

0.8 a

, 1/kPa

Macrostructural soil compressibility parameter (BBM in M1)

2.0×105 a

p,M1

Elastic stiffness of M1 for changes in net mean stress

0.1 a

S,M1

Elastic stiffness of M1 for changes in suction

(0)

Saturated slope of the virgin compression curve (BBM in M1) 0.3 b



Slope of the critical state line (BBM in M1)

M1

Poisson's ratio of M1

M2

Poisson's ratio of M2

ro

-p

re

lP

of

k

0.05 a

1.07 a 0.35 a 0.35 a

Data from Navarro et al. (2017). b Data from Toprak et al. (2018).

ur na

a

Value

3.3 Numerical model

To solve the boundary problem defined by the water mass conservation in M2, M1 and m

Jo

(Eq. [1]), as well as the mechanical equilibrium equation (Eq. [42]), the implementation platform provided by the program Comsol Multiphysics (CM) (Comsol, 2015), a partial differential equation solver based on the application of the finite element method with Lagrange multipliers (Babuška, 1973), was used. The program includes prebuilt models that allow to solve the problem without carrying out any numerical development. However, if this option is adopted, a previously programmed 24

constitutive formulation must be adopted, and therefore the desired conceptual model cannot be introduced. In addition, the structure of the differential equation to be solved is also predefined. Nevertheless, CM also includes modules that do allow to define both the differential equation to be solved and the constitutive formulation to be adopted. The researcher focuses on the definition of the physical problem and the software takes automatic control of assembling the system of equations and solving it. When this solving

of

strategy is adopted, CM should be understood as an implementation platform in which the

ro

researcher develops a new calculation program. This was the option adopted in the present paper, solving Eqs. [1] and [42], and implementing the mechanical model, the equations of

-p

motion and the state functions defined in the previous Sections.

re

The multiphysics logic of CM presents an important advantage. Without any additional programming effort, any of the balance equations can be disregarded in the calculation. If it

lP

is decided not to solve a balance equation, the associated state variable remains constant, equal to its initial value throughout the problem. Thus, for example, if it is selected not to

ur na

solve the water mass balance (M1, M2 and m), the evolution of the deformations would be calculated assuming that PL1, PL2 and em are constantly equal to the initial value introduced. Similarly, it would suffice not to solve the equilibrium equation to solve the flow in a nondeformable porous medium. Therefore, the numerical model can be simplified to match the

Jo

desired conceptual model. This fact, for example, allows to change from a DPM to a PTPM effortlessly, without compromising the numerical efficiency. But perhaps more relevant in this case is that, in the same way, the model can become more complex by adding new physical processes that, if solved, are automatically coupled to the rest. This was the strategy followed in the present work, i.e., the development of a PTPM from a previous DPM. 25

CM includes calculation processes based on the application of symbolic algebra. Thus, automatic symbolic differentiation techniques are applied (see, for instance, Martins and Hwang, 2013) to obtain the derivatives needed to define the iteration matrix. This significantly improves CM’s computational performance (Gobbert et al., 2009). The symbolic algebra and the symbolic differentiation present an additional advantage. Thanks to them, it is not necessary to explicitly define the numerous spatial and time derivatives

of

used in the formulation. Simply by indicating them, they are automatically calculated by

ro

CM. This makes the code clearer and simpler, facilitating its use by other researchers and reducing the implementation errors.

-p

However, the automatic symbolic differentiation of CM presents an important

re

disadvantage. If there is a state function that is defined implicitly (as occurs, for example, with the stress when an elastoplastic stress-strain model is adopted, or in many cases of

lP

non-linear elasticity), the program is not able to define its derivatives, resulting in an error that blocks calculation. Navarro et al. (2014) solved this problem using a mixed method

ur na

(Malkus and Hughes, 1978) considering stresses and plastic variables as new state variables. Satisfactory results were obtained (Alonso et al., 2012; Navarro et al., 2014, 2016, 2017, 2018), without finding the stability problems that occur in some cases of

Jo

application of mixed methods (Cervera et al., 2010).

4. Results and Discussion For the PTPM used to analyse the tests in Section 2, the material parameters in Tables 1 and 2 were used. Those associated with the flow of M1, since they were identified for compacted bentonite blocks (Navarro et al., 2017), should be adapted following the Appendix A to be used in the PTPM. For the sake of comparison, the tests were also 26

simulated using a DPM, merging M2 and M1 into a single Macrostructural “M” with the same strategy as in Gens et al. (2011). The BExM approach used by Navarro et al. (2017) was applied. The first test analysed, carried out by Clay Technology (Sweden) (Åkesson et al., in prep.), was focused on the hydraulic behaviour of an open pellet mixture (initial dry density of 953

of

kg/m3), characterising the time evolution of the water uptake and the change of the relative humidity RH in three different points along time (located at 3, 8 and 15.5 cm from the wet

ro

end, sensors 1, 2 and 3).

The initial water inflow in this test was very rapid, Fig. 4, showing an almost instant

-p

flooding of the sample to a certain height as a response to the water pressure applied at the

re

bottom of the sample. For this reason, in the simulation, the largest sample voids were taken initially flooded to a height consistent with the initial water inflow. Both using the

lP

PTPM and the DPM, the computed evolution of the water inflow with time matched closely the experimental results. However, to obtain this fit using DPM the value of the elastic

ur na

stiffness for suction changes had to be modified, taking S,M = 0.001 (as did Alonso et al., 2011), and the reference porosity O (Table 4), as in Toprak et al. (2018), was taken equal

Jo

to 0.319.

27

of ro

-p

Fig. 4. Water inflow evolution with time in the first test.

Although the identification of O allows to reproduce satisfactorily the “total” water inflow

re

using the DPM, the characterisation of the evolution of the distributed behaviour is not so

lP

satisfactory. This can be seen in Fig. 5, in which the deviation of the prediction of RH on sensors 2 and 3 (located above the initially flooded area) obtained with the DPM is not

ur na

negligible. The deviation from the sampled final water content obtained using the DPM is also appreciable (Fig. 6). The space available for flow is the same in both models (those voids different to the microstructural void ratio). However, in the PTPM that space is differentiated into M1 and M2. The M2 void ratio at the beginning of the test is large

Jo

(greater than 1, Fig. 7). It should be noted that the initial values of eM2 and eM1em were determined from the experimental value of the bulk density of the pellets and the overall dry density of the mixture after the samples were prepared. To differentiate the values of eM1 and em, the state surface in Fig. 3 was used to calculate em, obtaining eM1 later. Since the overall dry density is low, the inter-pellet space has not been completely occupied at the end of the test (Fig. 8). Then, all along the test, the Macrostructural space (voids available 28

for water flow) is comprised of M1 and M2 voids. The inter-pellets space works as the flow path that contains the initial rapid inflow. However, because of its low air entry pressure (0.5 kPa, Table 4), and despite its higher saturated hydraulic conductivity, M2 has a very low hydraulic conductivity in the non-flooded area at the start of the test. Therefore, the liquid water flow mainly occurs through M1 along the test. Nevertheless, the DPM assumes that the flow occurs in all M, which includes M1 and M2. Thus, the space available for

of

water transport is overestimated. Therefore, to simulate the water inflow, it is necessary to

ro

reduce the hydraulic conductivity, which was achieved by increasing O. However, the

retention capacity of M was not changed, which remained equal to that of M1. Since the air

-p

entry pressure of M1 (8.7 MPa = 1/1, see Table 1) is much higher than that of M2, a lower

re

RH than the real one should be predicted to fit the degree of saturation. In addition, the degree of saturation “applies to” a higher void ratio (since in the DPM the void ratio of the

lP

Macrostructure includes eM1 and eM2), and a greater mass of water is thus calculated. Therefore, in the DPM an even lower value of the RH is predicted to reproduce the water

Jo

ur na

content that the sample experimentally contents.

29

of ro

Jo

ur na

lP

re

-p

Fig. 5. RH evolution with time in the first test (s1: sensor 1, s2: sensor 2, s3: sensor 3).

Fig. 6. Experimental and numerical water content distribution at the end of the first test.

30

of ro

-p

Fig. 7. Computed void ratio profiles at the beginning and at the end of the first test obtained

Jo

ur na

lP

re

in the PTPM.

Fig. 8. First test. Computed evolution of the total, M2, M1 and m void ratios with time at the bottom (A), the centre (B) and top (C) of the sample obtained in the PTPM.

31

This problem could be avoided by defining a modified retention curve model for M. However, throughout the sample, the modification should be adapted to the relative variation of eM1 and eM2 shown in Fig. 8, being lower as the pellets swell and the weight of M1 is greater. The simplification of merging M1 and M2 requires loading the weight of this simplification on the material parameters. The DPM parameters should be understood, rather than as material parameters, as fitting parameters whose value should vary along

of

space and time if a correct spatial and temporal prediction of the behaviour of the system is

ro

sought. The weighting will no longer be necessary if eMeM1 (eM20). However, this does not occur in the first test analysed, since the pellet mixture is very open.

-p

As shown in Fig. 9, in the second test, by CEA (France) (Bernachy-Barbe et al., in prep.),

re

eM2 is almost null. In this test, in addition to the evolution of the hydration process, the evolution of the axial stress and of the radial stress at heights of 60 and 80 mm was

lP

measured. Therefore, it is a qualification exercise not only of the hydraulic formulation, but also of the mechanical simulation capacity of the model. Starting from a considerably more

ur na

closed initial situation (initial dry density equal to 1520 kg/m3), the initial value of eM2 is significantly lower (0.363, see Fig. 8) than in the previous test. Thus, it was not necessary to reduce the intrinsic permeability of M, which was taken equal to that of M1. However, the value of the elastic stiffness for changes in suction was still different, assuming again a

Jo

value of S,M=0.001.

32

of ro

-p

Fig. 9. Second test. Computed void ratio profiles at the beginning and at the end of the test

re

obtained in the PTPM.

lP

The reduced initial value of eM2 causes this void space to close when em swells, Fig. 10. Therefore, in this type of mixture, the DPM provides better predictions. It allows not only

ur na

to simulate the integrated behaviour of the whole sample, as it is shown in the prediction of the water inflow in Fig. 11 (in which the initial water inflow was more progressive). It also provides good predictions of the of the final state of the system, as shown in Fig 12, where the final water content and dry density distributions are plotted, and in Fig. 13, in which the

Jo

value of the radial stresses are represented. However, although acceptable predictions of the final values are achieved using the DPM, the results obtained of the evolution of the radial stress in the transient stage in which eM2 and eM1 change are of lower quality than those obtained using the PTPM (Fig 11). To end this Section, it is interesting to note that, although the fits obtained with the PTPM are satisfactory, they could be improved by changing the tentative estimation of the 33

parameters of M2 indicated in Section 3.2. However, since the results obtained allow to illustrate the scope of the conceptual model, a special fitting effort was not considered

lP

re

-p

ro

of

necessary.

Fig. 10. Second test. Computed evolution of the total, M2, M1 and m void ratios with time

Jo

ur na

at the bottom (A), centre (B) and top (C) of the sample obtained in the PTPM.

34

ur na

lP

re

-p

ro

of

Fig. 11. Water inflow evolution with time in the second test.

Jo

Fig. 12. Final water content (a) and dry density distribution (b) at the end of the second test.

35

of ro

-p

Fig. 13. Evolution with time of axial stress and radial stress at heights of 60 and 80 mm in

re

the second test.

lP

5. Conclusions

The paper presents a triple porosity model of the hydro-mechanical behaviour of MX-80

ur na

bentonite pellet mixtures. The model is based on the assumption that the macroscopic behaviour of the system is conditioned by three functional levels, associated with the three principal porosity levels that form the multi-modal pore structure of pellet mixtures: (i) the space inside the clay particle aggregates (microstructure), (ii) the voids between the

Jo

aggregates inside the pellets (Macrostructure 1), and (iii) the space between the pellets (Macrostructure 2). An overlapping continua approach was adopted. The macroscopic behaviour of each functional level was described by means of an equivalent continuous medium, assuming that the three continuous media are overlapped occupying the same spatial domain. The hydro-mechanical behaviour of the microstructure and the Macrostructure 1, as well as 36

their interaction, were modelled adapting the formulation usual in the simulation of the behaviour of compacted blocks. In this first approximation, a simplified model was adopted for the Macrostructure 2, using a formulation developed for granular materials to describe the flow and a non-linear isotropic elastic approach to describe its mechanical behaviour. The model was applied to simulate two hydration tests. The first test, focused on the characterisation of the flow in an open system (initial dry density of 953 kg/m3), allowed

of

the qualification of the hydraulic formulation. With the second test, in which the system

ro

was initially more closed (initial dry density of 1520 kg/m3), the mechanical simulation

capacity was also checked. Despite the different structure of both systems, and despite the

-p

unrefined model used to characterise the behaviour of the Macrostructure 2, encouraging

re

fits were obtained in both cases using the same material parameters. This shows the consistency of the proposed conceptual framework.

lP

For the sake of comparison, the tests were also simulated using a double porosity model, in which a single Macrostructural level that includes the Macrostructures 1 and 2 was

ur na

assumed. Since both Macrostructural levels were not differentiated, the descriptive potential of this model is lower (information such as that shown in Figs. 6, 7, 8 and 9 could not be obtained). In addition, since the two Macrostructural levels were merged in a single functional level, the double porosity model obtained a worse characterisation of the time

Jo

evolution of the behaviour of the system when the volume of voids associated with the inter-pellets space was not negligible. However, if the pellet mixture is initially closed or the confinement makes the inter-pellets space disappear when the microstructure swells, the simulation capacity of double porosity models could be sufficient to characterise the performance of bentonite pellet mixtures.

37

Acknowledgements This study was funded in part by Posiva Oy project 2113799, and by FPU Grant FPU15/02655 from the Spanish Ministry of Education awarded to Ms. De la Morena. In addition, we would like to thank Clay Technology and CEA for allowing us to use the

Jo

ur na

lP

re

-p

ro

of

results of their experimental work.

38

References Åkesson, M., Goudarzi, R., Börgesson, L., In preparation. EBS TF – THM Modelling. Water transport in pellets-filled slots. Laboratory tests and task descriptions. SKB P-19-06. Svensk Kärnbränslehantering AB. Alonso, E.E., Gens, A., Josa, A., 1990. A constitutive model for partially saturated soils.

of

Géotechnique 40 (3), 405-430. https://doi.org/10.1680/geot.1990.40.3.405.

ro

Alonso, E.E., Vaunat, J., Gens, A., 1999. Modelling the mechanical behaviour of expansive

-p

clays. Eng. Geol. 54 (1-2), 173-183. https://doi.org/10.1016/S0013-7952(99)00079-4. Alonso E.E., Pereira, J.-M., Vaunat, J., Olivella, S., 2010. A microstructurally based

lP

https://doi.org/10.1680/geot.8.P.002.

re

effective stress for unsaturated soils. Géotechnique 60 (12), 913-925.

Alonso, E.E., Romero, E., Hoffmann, C., 2011. Hydromechanical behaviour of compacted

ur na

granular expansive mixtures: experimental and constitutive study. Géotechnique 61 (4), 329-344. https://doi.org/10.1680/geot.2011.61.4.329. Alonso, J., Navarro, V., Calvo, B., Asensio, L., 2012. Hydro-mechanical analysis of Co2 storage in porous rocks using a critical state model. Int. J. Rock Mech. Min. Sci. 54, 19-26.

Jo

https://doi.org/10.1016/j.ijrmms.2012.05.016. Andersson, L., Sandén, T., 2012. Optimization of backfill pellet properties. ÅSKAR DP2. Laboratory tests. SKB Report R-12-18. Svensk Kärnbränslehantering AB, Swedish Nuclear Fuel and Waste Management Co. See: www.skb.se/publikation/2571606/R-12-18.pdf, last accessed July 2019. 39

Babuška, I., 1973. The finite element method with Lagrangian multipliers. Numer. Math. 20 (3), 179-192. https://doi.org/10.1007/BF01436561. Barrenblatt, G.I., Zheltov, I.P., Kochina, I.N., 1960. Basic concepts in the theory of seepage of homogeneous liquids in fissured rocks [strata]. J. Appl. Math. Mech. 24 (5), 1286-1303. https://doi.org/10.1016/0021-8928(60)90107-6.

ro

pore-size distribution and porosity. Can. Geotech. J. 50 (4), 1-16.

of

Beckett, C.T.S., Augarde, C.E., 2013. Prediction of soil water retention properties using

-p

https://doi.org/10.1139/cgj-2012-0320.

Benson, C.H., Chiang, I., Chalermyanont, T., Sawangsuriya, A., 2014. Estimating van

re

Genuchten parameters α and n for clean sands from particle size distribution data. Geotechnical Special Publication (ASCE), Geo-Congress 2014, 410-427.

lP

https://doi.org/10.1061/9780784413265.033.

ur na

Bernachy-Barbe, F., Imbert, C., Conil, N., Guillot, W., Gatabin, C., Talandier, J., In preparation. Heterogeneities in the final state of MX-80 bentonite under granular form in laboratory scale resaturation tests.

Borja, R.I., Koliji, A., 2009. On the effective stress in unsaturated porous continua with

Jo

double porosity. J. Mech. Phys. Solid. 57 (8), 1182-1193. https://doi.org/10.1016/j.jmps.2009.04.014. Bosgiraud, J.M., Foin, R., 2016. DOPAS Work Package 3 – Deliverable 3.1. WP3 FSS Construction Summary Report. DOPAS Project. See:

40

www.posiva.fi/files/4491/DOPAS_D3_1_WP3_FSS_Construction_Summary_Report_final _2016.pdf, last accessed July 2019. Brooks, R.M., Corey, A.T., 1964. Hydraulic properties of porous media. Colorado State University, Fort Collins, Colorado. Cervera, M., Chiumenti, M., Codina, R., 2010. Mixed stabilized finite element methods in

ro

(37-40), 2559-2570. https://doi.org/10.1016/j.cma.2010.04.006.

of

nonlinear solid mechanics. Part I: Formulation. Comput. Method. Appl. Mech. Eng. 199

-p

COMSOL (2015). COMSOL Multiphysics Reference Manual, version 5.1.

Dueck, A., Nilsson, U., 2010. Thermo-hydro-mechanical properties of MX-80. Results

re

from advanced laboratory tests. SKB Technical Report TR-10-55. Svensk

lP

Kärnbränslehantering AB, Swedish Nuclear Fuel and Waste Management Co. See: www.skb.se/publication/2223073/TR-10-55.pdf, last accessed July 2019.

ur na

Edlefsen, N.E., Anderson, A.B.C., 1943. Thermodynamics of soil moisture. Hilgardia 15 (2), 31-298. https://doi.org/10.3733/hilg.v15n02p031. Ergun, S., 1952. Fluid flow through packed columns. Chem. Eng. Prog. 48 (2), 89-94.

Jo

Ewen, J., Thomas, H.R., 1989. Heating unsaturated medium sand. Géotechnique 39 (3), 455-470. https://doi.org/10.1680/geot.1989.39.3.455. Gens, A., Alonso, E.E., 1992. A framework for the behaviour of unsaturated expansive clays. Can. Geotech. J. 29 (6), 1013-1032. https://doi.org/10.1139/t92-120.

41

Gens, A., Valleján, B., Sánchez, M., Imbert, C., Villar, M.V., van Geet, M., 2011. Hydromechanical behaviour of a heterogeneous compacted soil: Experimental observations and modelling. Géotechnique 61 (5), 367-386. https://doi.org/10.1680/geot.SIP11.P.015. Gobbert, M.K., Churchill, A., Wang, G., Seidman, T.I., 2009. COMSOL Multiphysics for efficient solution of a transient reaction-diffusion system with fast reaction. In: Y. Rao (Ed.)

of

Proceedings of the COMSOL Conference 2009. Boston.

ro

Guimarães, L.D.N., Gens, A., Sánchez, M., Olivella, S., 2013. A chemo-mechanical

constitutive model accounting for cation exchange in expansive clays. Géotechnique 63 (3),

-p

221-234. https://doi.org/10.1680/geot.SIP13.P.012.

re

Hoffmann, C., Alonso, E.E., Romero, E., 2007. Hydro-mechanical behaviour of pellet

lP

mixtures. Phys. Chem. Earth 32, 832-849. https://doi.org/10.1016/j.pce.2006.04.037. Hueckel, T.A., 1992. Water-mineral interaction in hygromechanics of clays exposed to

ur na

environmental loads: a mixture-theory approach. Can. Geotech. J. 29 (6), 1071-1086. https://doi.org/10.1139/t92-124.

Huyakorn, B., Lester, A., Faust, C., 1983. Finite element techniques for modeling groundwater flow in fractured aquifers. Water Resour. Res. 19 (4), 1019-1035.

Jo

https://doi.org/10.1029/WR019i004p01019. Jacinto, A.C., Villar, M.V., Ledesma, A., 2012. Influence of water density on the waterretention curve of expansive clays. Géotechnique 62 (8), 657-667. https://doi.org/10.1680/geot.7.00127.

42

Kahr, G., Kraehenbuehl, F., Stoeckli, H.F., Muller-von Moos, M., 1990. Study of the water–bentonite system by vapour adsorption, immersion calorimetry and X-ray techniques. II. Heats of immersion, swelling pressures and thermodynamic properties. Clay Miner. 25 (4), 499-506. https://doi.org/10.1180/claymin.1990.025.4.08. Karnland, O., Olsson, S., Nilsson, U., 2006. Mineralogy and sealing properties of various

of

bentonites and smectite-rich clay minerals. SKB Technical Report TR-06-30. Svensk Kärnbränslehantering AB, Swedish Nuclear Fuel and Waste Management Co. See:

ro

www.skb.com/publication/1419144/TR-06-30.pdf, last accessed July 2019.

-p

Khalili, N., Valliappan, S., Wan, C., 1999. Consolidation of fissured clay. Géotechnique 49

re

(1), 75-89. https://doi.org/10.1680/geot.1999.49.1.75.

Lade P.V., Nelson R.B., 1987. Modelling the elastic behaviour of granular materials. Int. J.

lP

Numer. Anal. Method. Geomech. 11 (5), 521-542. https://doi.org/10.1002/nag.1610110507.

ur na

Lloret, A., Villar, M.V., Sánchez, M., Gens, A., Pintado, X., Alonso, E.E., 2003. Mechanical behaviour of heavily compacted bentonite under high suction changes. Géotechnique 53 (1), 27-40. https://doi.org/10.1680/geot.2003.53.1.27. Lloret, A., Villar, M.V., 2007. Advances on the knowledge of the thermo-hydro-

Jo

mechanical behaviour of heavily compacted "FEBEX" bentonite. Phys. Chem. Earth 32 (814), 701-715. https://doi.org/10.1016/j.pce.2006.03.002. Luterkort, D., Johannesson, L.E., Eriksson, P., 2017. Buffer design and installation method. Installation report. SKB Technical Report TR-17-06. Svensk Kärnbränslehantering AB,

43

Swedish Nuclear Fuel and Waste Management Co. See: www.skb.com/publication/2490718/TR-17-06.pdf, last accessed July 2019. Mainka, J., Murad, M.A., Moyne, C., Lima, S.A., 2014. A modified effective stress principle for unsaturated swelling clays derived from microstructure. Vadose Zone J. 13 (5). https://doi.org/10.2136/vzj2013.06.0107.

of

Malkus, D.S., Hughes, T.J.R., 1978. Mixed finite element methods - Reduced and selective

-p

(1), 63-81. https://doi.org/10.1016/0045-7825(78)90005-1.

ro

integration techniques: A unification of concepts. Comput. Method. Appl. Mech. Eng. 15

Martins, J.R.R.A., Hwang, J.T., 2013. Review and unification of methods for computing

lP

https://doi.org/10.2514/1.J052184.

re

derivatives of multidisciplinary computational models. AIAA J. 51 (11), 2582-2599.

Mašín, D., Herbstová, V., Boháč, J., 2005. Properties of double porosity clayfills and

ur na

suitable constitutive models. Proceedings of the 16th International Conference of ICSMGE, Osaka, Japan, 12-14 September 2005, Vol. 2, pp. 827-830. Mašín, D., Khalili, N., 2016. Swelling phenomena and effective stress in compacted

Jo

expansive clays. Can. Geotech. J. 53 (1), 134-147. https://doi.org/10.1139/cgj-2014-0479. Molinero Guerra, A., Mokni, N., Delage, P., Cui, Y.J., Tang, A.M., Aimedieu, P., Bernier, F., Bornert, M., 2017. In-depth characterisation of a mixture composed of powder/pellets MX80 bentonite. Appl. Clay Sci. 135, 538-546. https://doi.org/10.1016/j.clay.2016.10.030. Najser, J., Mašín, D., Boháč, J., 2012. Numerical modelling of lumpy clay landfill. Int. J. Numer. Anal. Method. Geomech. 36 (1), 17-35. https://doi.org/10.1002/nag.990. 44

Navarro, V., Asensio, L., Alonso, J., Yustres, Á., Pintado, X., 2014. Multiphysics implementation of advanced soil mechanics models. Comput. Geotech. 60, 20-28. https://doi.org/10.1016/j.compgeo.2014.03.012. Navarro, V., Asensio, L., De la Morena, G., Pintado, X., Yustres, Á., 2015. Differentiated intra-and inter-aggregate water content models of mx-80 bentonite. Appl. Clay Sci. 118,

of

325-336. https://doi.org/10.1016/j.clay.2015.10.015.

ro

Navarro, V., Asensio, L., Yustres, Á., De la Morena, G., Pintado, X., 2016. Swelling and

https://doi.org/10.1016/j.enggeo.2016.01.005.

-p

mechanical erosion of MX-80 bentonite: Pinhole test simulation. Eng. Geol. 202, 99-113.

re

Navarro, V., Yustres, Á., Asensio, L., De la Morena, G., González-Arteaga, J., Laurila, T., Pintado, X., 2017. Modelling of compacted bentonite swelling accounting for salinity

lP

effects. Eng. Geol. 223, 48-58. https://doi.org/10.1016/j.enggeo.2017.04.016.

ur na

Navarro, V., De la Morena, G., González-Arteaga, J., Yustres, Á., Asensio, L., 2018. A microstructural effective stress definition for compacted active clays. Geomech. Energy Env. 15, 47-53. https://doi.org/10.1016/j.gete.2017.11.003. Nowamooz, H., Mrad, M., Abdallah, A., Masrouri, F., 2009. Experimental and numerical

Jo

studies of the hydromechanical behaviour of a natural unsaturated swelling soil. Can. Geotech. J. 46 (4), 393-410. https://doi.org/10.1139/T08-127. Olivella, S., Carrera, J., Gens, A., Alonso, E.E., 1994. Nonisothermal multiphase flow of brine and gas through saline media. Transp. Porous Media 15 (3), 271-293. https://doi.org/10.1007/BF00613282. 45

Olivella, S., Gens, A., 2000. Vapour transport in low permeability unsaturated soils with capillary effects. Transp. Porous Media 40 (2), 219-241. https://doi.org/10.1023/A:1006749505937. Philip, J.R., De Vries, D.A., 1957. Moisture movement in porous materials under temperature gradients. Eos, Trans. Am. Geophys. Union 38 (2), 222-232.

of

https://doi.org/10.1029/TR038i002p00222.

ro

Pollock, D.W., 1986. Simulation of fluid flow and energy transport processes associated

with high‐ level radioactive waste disposal in unsaturated alluvium. Water Resour. Res. 22

-p

(5), 765-775. https://doi.org/10.1029/WR022i005p00765.

re

Rutqvist, J., Börgesson, L., Chijimatsu, M., Kobayashi, A., Jing, L., Nguyen, T.S., Noorishad, J., Tsang, C.F., 2001. Thermohydromechanics of partially saturated geological

lP

media: Governing equations and formulation of four finite element models. Int. J. Rock

ur na

Mech. Min. Sci. 38 (1), 105-127. https://doi.org/10.1016/S1365-1609(00)00068-X. Saiyouri, N., Tessier, D., Hicher, P.Y., 2004. Experimental study of swelling in unsaturated compacted clays. Clay Miner. 39 (4), 469-479. https://doi.org/10.1180/0009855043940148. Salles, F., Douillard, J. M., Denoyel, R., Bildstein, O., Jullien, M., Beurroies, I., Van

Jo

Damme, H., 2009. Hydration sequence of swelling clays: Evolutions of specific surface area and hydration energy. J. Colloid Interface Sci. 333 (2), 510-522. https://doi.org/10.1016/j.jcis.2009.02.018.

46

Sánchez, M., Gens, A., Guimarães, L.D.N., Olivella, S., 2005. A double structure generalized plasticity model for expansive materials. Int. J. Numer. Anal. Method. Geomech. 29 (8), 751-787. https://doi.org/10.1002/nag.434. Sane, P., Laurila, T., Olin, M., Koskinen, K., 2013. Current Status of Mechanical Erosion Studies of Bentonite Buffer. Posiva Report 2012-45. Posiva Oy. See:

of

www.posiva.fi/files/3349/POSIVA_2012-45.pdf, last accessed July 2019.

https://doi.org/10.1346/CCMN.2013.0610601.

-p

management - A review. Clay. Clay Miner. 61 (6), 477-498.

ro

Sellin, P., Leupin, O.X., 2013. The use of clay as an engineered barrier in radioactive-waste

re

Sheng, D., Sloan, S.W., Gens, A., 2004. A constitutive model for unsaturated soils: Thermomechanical and computational aspects. Comput. Mech. 33 (6), 453-465.

lP

https://doi.org/10.1007/s00466-003-0545-x.

ur na

Sołowski, W.T., Hofmann, M., Hofstetter, G., Sheng, D., Sloan, S.W., 2012. A comparative study of stress integration methods for the Barcelona Basic Model. Comput. Geotech. 44, 22-33. https://doi.org/10.1016/j.compgeo.2012.03.007. Straughan, B., 2017. Mathematical Aspects of Multi–Porosity Continua. Advances in

Jo

Mechanics and Mathematics. Vol. 38. Springer, 208 pp. Svanadze, N., 2018. Fundamental solutions in the linear theory of thermoelasticity for solids with triple porosity. Math. Mech. Solid., In Press. https://doi.org/10.1177/1081286518761183.

47

Toprak, E., Olivella, S., Pintado, X., 2018. Modelling engineered barriers for spent nuclear fuel repository using a double-structure model for pellets. Environ. Geotech, In Press. https://doi.org/10.1680/jenge.17.00086. Tournassat, C., Appelo, C.A.J., 2011. Modelling approaches for anion-exclusion in compacted Na-bentonite. Geochim. Cosmochim. Acta 75 (13), 3698-3710.

of

https://doi.org/10.1016/j.gca.2011.04.001.

ro

van Genuchten, M.T., 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44 (5), 892-898.

-p

https://doi.org/10.2136/sssaj1980.03615995004400050002x.

re

Vilarrasa, V., Rutqvist, J., Blanco Martin, L., Brikholzer, J., 2016. Use of a dual-structure constitutive model for predicting the long-term behavior of an expansive clay buffer in a

lP

nuclear waste repository. Int. J. Geomech. 16 (6), D4015005.

ur na

https://doi.org/10.1061/(ASCE)GM.1943-5622.0000603. Wadsö, L., Svennberg, K., Dueck, A., 2004. An experimentally simple method for measuring sorption isotherms. Dry. Technol. 22 (10), 2427-2440. https://doi.org/10.1081/LDRT-200032898.

Jo

Wilson, R., Aifantis, E., 1982. On the theory of consolidation with double porosity. Int. J. Eng. Sci. 20 (9), 1019-1035. https://doi.org/10.1016/0020-7225(82)90036-2. Yong, R.N., 1999. Overview of modeling of clay microstructure and interactions for prediction of waste isolation barrier performance. Eng. Geol. 54 (1-2), 83-91. https://doi.org/ 10.1016/S0013-7952(99)00064-2. 48

Appendix A The specific discharge vector qM1 defines the volume of liquid water in M1 that, taking the displacement of the solid skeleton as a reference, crosses a unit of domain area per unit time. If the parameters identified for a block of compacted bentonite using a DPM are considered, the specific discharge disregarding the inter-pellet void space is calculated. Identifying this flow as qM1*, the following can be written:

dVWL,M1

of

qM1* dt 

AM1  Am  Amin

ro

[A1]

where qM1* defines the flow that crosses a surface A in which the area of M1 is equal to

-p

AM1, the area of m is equal to Am, and the area of mineral is equal to Amin. The equation considers a time differential dt in which a liquid water volume dVWL,M1 flows through AM1

qM1* dt 

dVWL,M1 AM2  AM1  Am  Amin

AM2  AM1  Am  Amin AM1  Am  Amin

lP

[A2]

re

(always in relation to the movement of the solid). Eq. [A1] can be formulated as:

where AM2 is the area of M2 in A. The first term on the right of this equation can be written

ur na

as qM1·dt, where qM1 is the normal flow to the total surface A (= AM2  AM1  Am  Amin). Therefore:

qM1* dt   qM1 dt 

AM2  AM1  Am  Amin AM1  Am  Amin

Jo

[A3]

Since the equations above are valid for any point in the domain and for any surface A, it follows that: [A4]

q M1* 

AM2  AM1  Am  Amin q M1 AM1  Am  Amin

49

In an isotropic media, and considering that the dimensions of the area are comparable to those of a representative elementary volume, then: [A5]

AM2  AM1  Am  Amin VM2  VM1  Vm  Vmin 1 e   AM1  Am  Amin VM1  Vm  Vmin 1  eM1  em

where VM2, VM1 and Vm are, respectively, the volumes of porosity associated with the structural levels M2 (inter-pellet), M1 (intra-pellet, not intra-aggregate) and m (intra-

of

aggregate), and Vmin is the volume of mineral present in the representative elementary

q M1 

1  eM1  em q M1* 1 e

-p

[A6]

ro

volume considered. Consequently:

If the parameters identified in a DPM are used in the calculation process, the specific

re

discharge obtained should be multiplied by the reducing coefficient (1 + eM1 + em) / (1 + e) to obtain the specific discharge to be used in the PTPM.

rmM1 dt 

dVW,mM1

VM2  VM1  Vm  Vmin

ur na

[A7]

lP

Something similar occurs with the mass exchange rmM1, which can be formally defined as:

i.e., the exchange of liquid water mass dVWmM1 from M1 to m per unit of total volume of soil in a time differential dt. If the parameters identified in a DPM are used, rmM1* is actually calculated, i.e., the exchange per unit volume if the porosity associated with M2

Jo

were not present: [A8]

rmM1* dt 

dVW,mM1 VM1  Vm  Vmin

Therefore, from the above two equations it follows that: [A9]

rmM1 

1  eM1  em rmM1* 1 e 50

This equation shows that, to convert rmM1* to rmM1, the same reducing coefficient should

Jo

ur na

lP

re

-p

ro

of

be used as to convert qM1* to qM1 (Eq. [A6]).

51