Z Michalec

R Blaheta
email: radim.blaheta@ugn.cas.cz

M Hasal

T Ligurský

Fully coupled thermo-hydro-mechanical model with oversaturation and its validation to experimental data from FEBEX experiment

Keywords: modelling of coupled thermo-hydro-mechanical (THM) processes, specific aspects of THM processes in bentonite, oversaturation, implementation

The paper is devoted to the formulation of a mathematical model which enables analysis of coupled thermo-hydro-mechanical (THM) processes with a special focus on processes in bentonite and bentonite-based sealing barriers which are component of designs of deep geological repositories for a high-level radioactive waste including the spent nuclear fuel. Such modelling is important as the analysis of coupled THM processes in the near field is crucial for the understanding of the long-term behaviour of the barriers. The developed model uses generally available data from standard tests but attempts to incorporate special characteristics of the bentonite behaviour, let us mention the oversaturation, swelling and strong couplings of the processes. The model is implemented via COMSOL Multiphysics and validated by the simulation of THM processes which were monitored during operation of the FEBEX experiment and analysed during two dismantling phases of the experiment.

Introduction

The paper is motivated by the need for analysis of coupled thermo-hydro-mechanical (THM) processes which is necessary for the assessment of performance and safety of deep geological high-level radioactive waste storage facilities, see e.g. [START_REF] Pusch | Geological Storage of Highly Radioactive Waste. Current Concepts and Plans for Radioactive Waste Disposal[END_REF]. The need for coupled THM analysis is especially important for processes in the near field, close to the deposited canisters with the nuclear waste (the spent nuclear fuel), i.e. in the bentonite-based sealing barriers and the surrounding host rock. This near field is in the initial stage of repository operation strongly influenced by the heat produced by the spent nuclear fuel; the bentonite barrier is slowly water-saturated from the surrounding (host) rock mass and influenced by the forces arising from swelling of the bentonite.

The processes in the host rock can be relatively well described by classical models for flow in porous media, but a special interest should be devoted to the modelling of the processes in bentonite barriers. This paper is devoted to the development of a reasonably simple THM model, which however considers the specific properties of the bentonite and provides the capability of simulation of THM processes in the bentonite sealing barriers. The specific features of the bentonite, as swelling and oversaturation induced by the complicated microstructure, will be taken into consideration. Moreover, these specific features will be described on the base of data from relatively standard laboratory tests.

The capability of simulating the coupled processes will be demonstrated by simulation of processes occurring in the FEBEX experiment, i.e. Full-scale Engineered Barriers Experiment performed in the Grimsel test site in Switzerland.

Data available from monitoring THM processes during more than 18-year operation of the FEBEX experiment as well as the data from two phases of dismantling the experiment (after 5 and 18 years), will serve for validation of the described model. A selection of required data was available due to the participation in the international DECOVALEX project. The necessary data, as well as methodology for validation, were prepared by A. Gens, who was the coordinator of the relevant task of the project. The DECOVALEX project also offers the possibility of comparing the modelling technique described in this paper with modelling used by other teams with their approaches and different software, see [START_REF] Gens | DECOVALEX-2019[END_REF], [START_REF] Gens | HM and THM interactions in bentonite engineered barriers for nuclear waste disposal[END_REF]. Note that modelling related to FEBEX experiment was also investigated in an earlier phase of the DECOVALEX project, see [START_REF] Alonso | The FEBEX benchmark test: Case definition and comparison of modelling approaches[END_REF], [START_REF] Stephansson | Coupled Thermo-Hydro-Mechanical-Chemical Processes in Geo-Systems: Fundamentals, Modelling, Experiments and Applications[END_REF].

Many papers provide a description of the FEBEX experiment and the material behaviour of the FEBEX bentonite and the granite host rock for the FEBEX experiment [START_REF] Fuentes-Cantillana | FEBEX: Full-scale engineered barriers experiment in crytalline host rock: Final design and installation of the 'in situ' test at Grimsel[END_REF], [START_REF] Alberdi | FEBEX project: full-scale engineered barriers experiment for a deep geological repository for high level radioactive waste in crystalline host rock[END_REF], [START_REF] Bárcena | Dismantling of heater No. 1 at the FEBEX in situ test. Description of operations[END_REF], [START_REF] Lloret | FEBEX II Project Final report on thermo-hydro-mechanical laboratory tests[END_REF], [START_REF] Huertas | Full-scale engineered barrier experiment[END_REF], [START_REF] García-Siñeriz | FEBEX DP: Dismantling of heater 2 at the FEBEX "in situ" test. Description of operations[END_REF], [START_REF] Villar | FEBEX-DP Post-mortem THM/THG Analysis Report[END_REF] and [START_REF] Villar | Report on Thermo-Hydro-Mechanical Laboratory Tests Performed by CIEMAT on Febex Bentonite 2004-2008[END_REF]. The development of mathematical models for the simulation of THM processes occurring during the experiment operation is described e.g. in [START_REF] Alonso | The FEBEX benchmark test: Case definition and comparison of modelling approaches[END_REF], [START_REF] Stephansson | Coupled Thermo-Hydro-Mechanical-Chemical Processes in Geo-Systems: Fundamentals, Modelling, Experiments and Applications[END_REF], [START_REF] Rutquist | A Fully Coupled Three-Dimensional THM Analysis of the Febex in Situ Test with the Rocmas Code: Prediction of THM Behavior in a Bentonite Barrier[END_REF], [START_REF] Gens | A full-scale in situ heating test for high-level nuclear waste disposal: Observations, analysis and interpretation[END_REF] and [START_REF] Papafotiou | Predismantling THM modelling of the FEBEX in situ experiment[END_REF], the realisation of such models uses different software tools. The modelling, which was done during the phase DECOVALEX 2019 by using the software COMSOL Multiphysics, THAMES, TOUGH2+FLAC and HYDROGEOMECH, is described in [START_REF] Gens | DECOVALEX-2019[END_REF], [START_REF] Gens | HM and THM interactions in bentonite engineered barriers for nuclear waste disposal[END_REF]. The model developed in this paper is implemented with the aid of COMSOL Multiphysics software, which provides enough tools and capabilities for simulation of complex coupled problems.

The specific aims of this paper are (1) treating of bentonite oversaturation by using the available data, [START_REF] Gens | DECOVALEX-2019[END_REF] showing a significant influence of different expressions for vapour diffusion, [START_REF] Gens | HM and THM interactions in bentonite engineered barriers for nuclear waste disposal[END_REF] demonstration of the possibility of numerical implementation using the COMSOL Multiphysics software, (4) simulation of THM processes during FEBEX experiment and comparison of the computed results with data from monitoring and dismantling.

The rest of the paper is organised as follows. Section 2 describes the basic principles for deriving separated THM models and their couplings which are considered in this paper. Section 3 is devoted to the explanation of oversaturation and swelling, construction of retention function suitable for FEBEX bentonite, discussion of the choice of the vapour diffusion coefficient and construction of other constitutive relations using the available data. Section 4 shortly describes the numerical implementation via COMSOL Multiphysics and the solvers offered there. The FEBEX experiment is briefly described in Section 5, the procedure of modelling the processes which occur in FEBEX in situ test is described in Section 6 including several phases of the test with changing geometry, necessary transfer of data between the following phases etc. Section 7 describes the outputs from the modelling and discusses the results, including the effect of different expressions for vapour diffusion and the correspondence with the measured data. The paper ends with general conclusions in Section 8.

Conservation equations and coupled THM model

This section aims to develop a THM model as simple as possible, but suitable for simulation of processes in bentonite sealing with various applications as it is demonstrated later in the paper. The focus is on the saturation and swelling of the bentonite which consists of solid part (s), liquid water (w), water vapour (v) and air (a). The pressure in gas (vapour and air) is assumed constant and consequently a Richards type, one-phase model is used for the description of flow in variably saturated bentonite. This section introduces the conservation equations for separate THM processes, then the interconnections among THM processes are mainly considered in Section 3. To keep the model relatively simple, we neglect some couplings which we consider weaker, not very significant. As general references for some details related to the provided description of the THM model, we can mention [START_REF] Bear | Modeling Groundwater Flow and Contaminant Transport[END_REF], [START_REF] Lewis | The Finite Element Method in the Static and Dynamic Deformation and Consolidation of Porous Media[END_REF] and also [START_REF] Abed | A study on how to couple thermo-hydro-mechanical behaviour of unsaturated soils: Physical equations, numerical implementation and examples[END_REF] and the references therein.

The model is described by equations which hold in the space-time domain Ω 𝑡 = Ω × (0, 𝑡 𝑚𝑎𝑥 ).

For brevity of this section, we shall not discuss the boundary and initial conditions, but they must be included in the implementation of the model.

Water and vapour flow model

We consider a situation when water appears in two forms of liquid water and vapour. The mass conservation will be described by one equation [START_REF] Pusch | Geological Storage of Highly Radioactive Waste. Current Concepts and Plans for Radioactive Waste Disposal[END_REF], which can be derived as a sum of two mass Using the above equations, the mass conservation equation for the liquid water and vapour is considered in the following form

𝜕 𝜕𝑡 (𝑤𝜌 𝑑 ) = -∇ ⋅ (𝜌 𝑤 𝑣 𝑤 + 𝑞 𝑣 ) + 𝑄 𝑤 . (1)
The equation ( 1) represents a generalised Richards equation including non-isothermal vapour diffusion. A constant (atmospheric) gas pressure is assumed, 𝑝 𝑔 = 𝑝 𝑎𝑡𝑚 ≈ 0. This assumption corresponds with vapour transport only through diffusion and enables to formulate one-phase model without separate mass balance equation for the gas flow. For more details concerning the Richards equation, see e.g. [START_REF] Bear | Modeling Groundwater Flow and Contaminant Transport[END_REF]. The Richards and two-phase flow approaches are compared in [START_REF] Wang | Non-isothermal flow in low permeable porous media: a comparison of 'Richards' and two-phase flow approaches[END_REF].

In [START_REF] Pusch | Geological Storage of Highly Radioactive Waste. Current Concepts and Plans for Radioactive Waste Disposal[END_REF], 𝑣 𝑤 is the Darcy velocity

𝑣 𝑤 = - 𝑘 𝑠 𝜇 𝑘 𝑟 (𝑆 𝑒 )∇𝑝 𝑤 , (2)
where 𝑝 𝑤 is the pore water pressure, 𝑝 𝑤 ≥ 0 and 𝑝 𝑤 < 0 for the water-saturated and unsaturated case, respectively, 𝑘 𝑠 is the (absolute, intrinsic) permeability in the fully saturated case, 𝑘 𝑟 = 𝑘 𝑟 (𝑆 𝑒 )

is the relative permeability depending on the effective saturation with water 𝑆 𝑒 , 𝜇 is the dynamic viscosity of water.

We shall use definitions of the water saturation (𝑆 𝑤 ) and effective saturation (𝑆 𝑒 ) which are not related to volumes and porosity but which are expressed through the mass of water as follows

𝑆 𝑤 = 𝑤 𝑤 𝑚𝑎𝑥 , 𝑆 𝑒 = 𝑤 -𝑤 𝑟𝑒𝑠 𝑤 𝑚𝑎𝑥 -𝑤 𝑟𝑒𝑠 , (3)
where 𝑤 𝑚𝑎𝑥 and 𝑤 𝑟𝑒𝑠 are maximum and residual (minimum) water mass content. The definition of 𝑆 𝑤 and 𝑆 𝑒 given in (3) enables to include the oversaturation into the developed model. More detailed discussion on this topic and relation to saturation defined in a more usual way through (bulk) porosity is provided in the next section, see the equation [START_REF] Villar | Report on Thermo-Hydro-Mechanical Laboratory Tests Performed by CIEMAT on Febex Bentonite 2004-2008[END_REF].

Note that the saturated permeability is considered as depending on the dry density, 𝑘 𝑠 = 𝑘 𝑠 (𝜌 𝑑 ), and through dry density also on the temperature, see Section 3.7. The dynamic viscosity also depends on the temperature 𝑇, 𝜇 = 𝜇(𝑇). These topics are also discussed in the next section.

Assuming that 𝑤 depends on 𝑝 𝑤 , 𝜌 𝑑 and 𝑇, the change of the water mass can be expressed as

In the unsaturated case, the terms , 𝜀 𝑒𝑙,𝑣 is defined in Section 2.3 and more details can be found in Section 3.7.

Following [START_REF] Philip | Moisture Movement in Porous Materials under Temperature Gradients[END_REF], the vapour mass flux 𝑞 𝑣 can be primarily expressed as

𝑞 𝑣 = -𝐷 𝑣 ∇𝜌 𝑣 , (5)
where the flux is driven by the gradient of the mass density of vapour and 𝐷 𝑣 is the diffusion coefficient. The vapour density 𝜌 𝑣 can be expressed from the Kelvin-Laplace equation

𝜌 𝑣 = 𝜑𝜌 𝑣𝑠𝑎𝑡 (𝑇), 𝜑 = exp ( 𝑝 𝑤 𝜌 𝑤 𝑟 𝑣 𝑇
), see e.g. [START_REF] Philip | Moisture Movement in Porous Materials under Temperature Gradients[END_REF] and [START_REF] Lewis | The Finite Element Method in the Static and Dynamic Deformation and Consolidation of Porous Media[END_REF], where 𝜑 is relative humidity, 𝜌 𝑣𝑠𝑎𝑡 (𝑇) is the density for full vapour saturation at the temperature 𝑇, 𝑟 𝑣 = 𝑅/𝑀 𝑚𝑤 is the specific gas constant for water vapour, 𝑅 is the universal gas constant and 𝑀 𝑚𝑤 is the molecular weight of water. From the above formula, we get

Thus, the vapour mass flux has a part driven by the pore pressure gradient and a part driven by the temperature gradient. The determination of the diffusion coefficient is discussed in Section 3.3. Above, 𝜂 = 1 gives the expression (5) and will be used in our model; 𝜂 > 1 is a possibility (enhancement factor) introduced already in [START_REF] Philip | Moisture Movement in Porous Materials under Temperature Gradients[END_REF].

Heat transport

We shall assume that the heat propagation is solely done by the heat conduction; the heat from phase changes of water and the heat transport by convection are neglected. Therefore, the heat balance is described by the following equation,

𝑐𝜌 𝜕𝑇 𝜕𝑡 = ∇ ⋅ (𝜆 ∇𝑇) + 𝑄 𝑇 , (7)
where c is the heat capacity, 𝜌 is the mass density, 𝜆 is the heat conduction coefficient and 𝑄 𝑇 denotes the density of heat sources.

The heat capacity and the mass density in [START_REF] Alberdi | FEBEX project: full-scale engineered barriers experiment for a deep geological repository for high level radioactive waste in crystalline host rock[END_REF] concern the mixture of the solid matrix and water (contribution of the gaseous phase can be neglected), so that they are computed as averages of the corresponding quantities, i.e. where the already introduced notation is used.

𝑐 = 𝑐

Mechanics

The deformation of the porous body Ω is assumed to be described by the displacement 𝑢 and the small strain tensor 𝜀 = 𝜀(𝑢). Moreover, 𝜎 will denote the Cauchy stress tensor. The sign convention of continuum mechanics is used, i.e. expansive strains and tensile stresses are positive.

It is assumed that the triple (𝑢, 𝜀, 𝜎) evolves in time in balance with the acting volume forces 𝑓. The time evolution of the system is due to time-dependent swelling and temperature changes. This evolution is slow so that inertia forces can be neglected, and therefore the momentum conservation has the form -∇ • 𝜎(𝑥, 𝑡) = 𝑓(𝑥) for (𝑥, 𝑡) ∈ Ω 𝑡 = Ω × (0, 𝑡 𝑚𝑎𝑥 ).

We shall assume that the total stress has three components

𝜎 = 𝜎 𝑡𝑜𝑡𝑎𝑙 = 𝜎 𝑒𝑓𝑓𝑒𝑐𝑡 -𝜎 𝑤𝑎𝑡𝑒𝑟 -𝜎 𝑠𝑤𝑒𝑙𝑙 , (9)
where 𝜎 𝑒𝑓𝑓𝑒𝑐𝑡 is the (effective) stress related to elastic deformation of the solid matrix, 𝜎 𝑤𝑎𝑡𝑒𝑟 is hydrostatic stress from the pore pressure and 𝜎 𝑠𝑤𝑒𝑙𝑙 is the stress due to swelling. The (total) strain is composed of two components:

𝜀 = 𝜀 𝑡𝑜𝑡𝑎𝑙 = 𝜀(𝑢) = 𝜀 𝑒𝑙 + 𝜀 𝑡ℎ , (10)
which are the elastic strain 𝜀 𝑒𝑙 and the isotropic thermal expansion 𝜀 𝑡ℎ = 𝛼 𝑣 3 (𝑇 -𝑇 𝑟𝑒𝑓 )𝐼. Here, 𝛼 𝑣 is the coefficient of volumetric thermal expansion, 𝐼 is the unit tensor, 𝑇 𝑟𝑒𝑓 is the reference (initial) temperature. Note that small strains are assumed, 𝜀 = (𝜀 𝑖𝑗 ), 𝜀 𝑖𝑗 = ( Note also that the swelling is expressed through the stress.

We shall assume the strain-stress relation with the elastic tensor 𝐶, 𝜎 𝑒𝑓𝑓𝑒𝑐𝑡 (𝑥, 𝑡) = 𝐶(𝑥, 𝑡): 𝜀 𝑒𝑙 (𝑥, 𝑡), [START_REF] García-Siñeriz | FEBEX DP: Dismantling of heater 2 at the FEBEX "in situ" test. Description of operations[END_REF] where, in the case of isotropic material, the inverse relation to [START_REF] García-Siñeriz | FEBEX DP: Dismantling of heater 2 at the FEBEX "in situ" test. Description of operations[END_REF] can be described in the form

𝜀 𝑒𝑙,𝑖𝑗 = 1+𝜈 𝐸 𝜎 𝑖𝑗 - 𝜈 𝐸 σ 𝑒𝑓,𝑣 𝛿 𝑖𝑗 , σ 𝑒𝑓,𝑣 = ∑ 𝜎 𝑖𝑖 𝑖 ,
where 𝜎 𝑒𝑓𝑓𝑒𝑐𝑡 = (𝜎 𝑖𝑗 ), 𝜀 𝑒𝑙 = ( 𝜀 𝑒𝑙,𝑖𝑗 ), 𝐸 is the elasticity modulus, 𝜈 is the Poisson ratio, 𝜎 𝑒𝑓,𝑣 denotes the volumetric component of the effective stress and 𝛿 is the Kronecker delta. In the case of some kinds of material nonlinearity, a variation of 𝐶 = 𝐶(𝑥, 𝑡, 𝜀(𝑥, 𝑡)) is admitted, see e.g. [START_REF] Duncan | Modern techniques for the analysis of engineering problems in soft clay[END_REF], [START_REF] Lewis | The Finite Element Method in the Static and Dynamic Deformation and Consolidation of Porous Media[END_REF]. Later in ( 24) and ( 25), we consider the evolution of elastic modulus in time, in dependence on dry density (involving strains) and suction (water saturation). We assume that plasticity effects can arise only as an effect of an increase of loading and therefore, the deformation theory of plasticity (materially nonlinear elasticity) is sufficient. Nevertheless, the model could be generalised for adopting loading/unloading behaviour by involving the Basic Barcelona Model (BBM) and its generalisation to non-isothermal states, see e.g. [START_REF] Alonso | A constitutive model for partially saturated soils[END_REF], [START_REF] Gens | Constitutive laws. In Modern Issues[END_REF].

The pore pressure contribution is as follows,

𝜎 𝑤𝑎𝑡𝑒𝑟 = 𝛼𝜒𝑝 𝑤 𝐼, ( 12
)
where 𝛼 is the Biot-Willis constant, which is assumed to be equal to one in this paper, and 𝜒 is 'Bishop's function depending on the saturation. In the fully saturated case 𝜒 = 1. In the unsaturated case, a standard 𝜒(𝑆 𝑤 ) = 𝑆 𝑤 , see [START_REF] Bear | Modeling Groundwater Flow and Contaminant Transport[END_REF], [START_REF] Lewis | The Finite Element Method in the Static and Dynamic Deformation and Consolidation of Porous Media[END_REF], is used for granite and concrete, but for the bentonite, we use 𝜒 = 0, see e.g. [START_REF] Chijimatsu | Model development and calibration for the coupled thermal, hydraulic and mechanical phenomena of the bentonite[END_REF]. Correspondingly, the source term 𝑄 𝑤 = 𝛼𝜒 𝜕𝜀 𝑒𝑙,𝑣 𝜕𝑡 , see also Section 3.7. The swelling stress 𝜎 𝑠𝑤𝑒𝑙𝑙 is described in the next section.

Constitutive relations and THM interactions

This section deals with the constitutive relations, which also describe the most significant couplings among the THM processes in the bentonite. Some more details can be found e.g. in [START_REF] Lewis | The Finite Element Method in the Static and Dynamic Deformation and Consolidation of Porous Media[END_REF] and [START_REF] Abed | A study on how to couple thermo-hydro-mechanical behaviour of unsaturated soils: Physical equations, numerical implementation and examples[END_REF].

Inner structure and specific properties of the bentonite

The term of oversaturation was already mentioned in Section 1 for the phenomenon that the bentonite is capable of storing more water than the amount corresponding to its standard (bulk) pore space 𝑉 𝑏 = 𝑉 -𝑉 𝑠 where 𝑉 𝑠 is the volume of solid part within 𝑉. The corresponding porosity 𝜙 can be defined as

𝜙 = 𝑉 -𝑉 𝑠 𝑉 = 1 - 𝜌 𝑑 𝜌 𝑠 .
The effect of oversaturation is a consequence of a complicated inner structure, see Fig. 1, and chemistry acting in the bentonite or, especially, in its component, the montmorillonite. The influence of this inner structure with (1) bulk pores, (2) interlayers and (3) double layers on the flow and diffusion could be explained by a triple porosity model with specific flow (diffusion) mechanism for each component -bulk pore water, interlayer water and double layer water with proper exchange rules for these water components. Also, chemistry, charges and ion system could be considered. For more details see e.g. [START_REF] Kröhn | Modelling the re-saturation of bentonite in final repositories in crystalline rock[END_REF], [START_REF] Kröhn | New evidence for the dominance of vapour diffusion during the re-saturation of compacted bentonite[END_REF], [START_REF] Idiart | Models for Diffusion in Compacted Bentonite[END_REF] and the references therein. However, such description requires much knowledge, experimental data and research which is hardly available for the considered applications.

Therefore, we suggest to replace the quantities connected with porosity with quantities related to the mass of water like the water mass content 𝑤, see the definition of saturation (3). This definition of saturation avoids the necessity of analysis of more complicated pore structure and also the deformation of pores during the saturation process.

Using this idea, we construct a model using the available macroscopic data, which are, of course, affected by the mentioned inner structure. The saturation enters the model many times, explicitly or implicitly, see the description of the retention function, the relative permeability, the vapour diffusion coefficient, the stress-strain relation and the swelling stress in the following subsections.

Fig. 1: Free (bulk) pore water, interlayer water and double layer water, see [START_REF] Bradbury | Porewater chemistry in compacted re-saturated MX-80 bentonite[END_REF].

An important value 𝑤 𝑚𝑎𝑥 can be derived from the data provided for the construction of the retention function, see the next subsection. The accuracy of determination of 𝑤 𝑚𝑎𝑥 is influenced by the fact if data are available for the states close to full saturation (states with very low suction).

Note that many porous materials have less complicated inner structure and therefore can work with saturation 𝑆 𝑟 related to the bulk pore space and defined as follows

𝑆 𝑟 = 𝑤 𝑤 𝑝𝑜𝑟𝑒 , 𝑤 𝑝𝑜𝑟𝑒 = 𝜌 𝑤 𝜙 𝜌 𝑑 , 𝑆 𝑟 = 𝑤 𝑚𝑎𝑥 𝑤 𝑝𝑜𝑟𝑒 𝑆 𝑤 ≥ 𝑆 𝑤 . (13)
Above, the value 𝑤 𝑝𝑜𝑟𝑒 corresponds to the full saturation of the bulk pores and obviously 𝑤 𝑝𝑜𝑟𝑒 ≤ 𝑤 𝑚𝑎𝑥 , but a sharp inequality is typical for the bentonite. The third identity in [START_REF] Villar | Report on Thermo-Hydro-Mechanical Laboratory Tests Performed by CIEMAT on Febex Bentonite 2004-2008[END_REF] provides a relation between 𝑆 𝑟 and 𝑆 𝑤 , which indicates that 𝑆 𝑟 can be greater than one. Values of 𝑆 𝑟 > 1 are confirmed by measurements, which can be seen e.g. in Fig. 5 and Fig. 13.

Retention function

We fit the data provided for the retention relation with a modified van Genuchten formula of the form

w = 𝑓 𝑟 (𝑠, 𝜌 𝑑 , 𝑇) = 𝑤 𝑚𝑎𝑥 (𝜌 𝑑 , 𝑇) ⋅ (1 + ( 𝑠 𝑃 0 (𝜌 𝑑 ) ) 1 1- 1 ) - 1 ⋅ (1 - 𝑠 𝑃 00 )  0 for 0 ≤ 𝑠 ≤ 𝑃 00 , (14)
where 𝜆 1 , 𝜆 0 are dimensionless parameters and 𝑃 0 , 𝑃 00 are parameters expressed in MPa. 𝑃 00 is the closing pressure, 𝑠 = 𝑝 𝑔 -𝑝 𝑤 ≈ -𝑝 𝑤 is the suction in MPa. The suction s= 0 provides 𝑤 = 𝑤 𝑚𝑎𝑥 , s= 𝑃 00 provides 𝑤 = 0. More precisely, w = 𝑓 𝑟 (𝑠, 𝜌 𝑑 , 𝑇) = 𝑤 𝑚𝑎𝑥 (𝜌 𝑑 , 𝑇) for full saturation 𝑠 ≤ 0, w = 𝑓 𝑟 (𝑠, 𝜌 𝑑 , 𝑇) = 𝑤 𝑟𝑒𝑠 ≈ 0 for 𝑠 > 𝑃 00 .

Note that 𝑤 𝑟𝑒𝑠 , appearing in the definition of the effective saturation (3), is assumed to be 𝑤 𝑟𝑒𝑠 = 0.001 ≈ 0.

The expression [START_REF] Rutquist | A Fully Coupled Three-Dimensional THM Analysis of the Febex in Situ Test with the Rocmas Code: Prediction of THM Behavior in a Bentonite Barrier[END_REF] was fitted to the measured data for the FEBEX bentonite provided for different dry densities and different temperatures. It provides the following expression

𝑤 𝑚𝑎𝑥 (𝜌 𝑑 , 𝑇) = 𝑒 -0.015⋅(𝑇-𝑇 0 )  𝑑 3 100 + 𝜌 𝑤 𝜙  𝑑
, 𝑃 0 (𝜌 𝑑 ) = 10 -5 ⋅ 𝑒 (8,2⋅ 𝑑 ) , 𝑃 00 = 1000,

𝜆 1 = 0.14 and 𝜆 0 = 1.5, the dry density and specific mass density of water are in the unit Mg/m 3 . The term 𝜌 𝑤 𝜙/ 𝑑 represents the mass water content in the bulk pores, see ( 13), 𝑒 -0.015⋅(𝑇-𝑇 0 )  𝑑 3 /100 is an expression for the mass water content in extra pores, 𝑇 0 = 20 ℃ is the initial (reference) temperature.

A comparison of the fitted expression ( 14) with ( 16) with experimental data can be seen in Fig. 2. Note that for the porous material like granite or concrete, the value 𝑤 𝑚𝑎𝑥 is fully determined by the bulk porosity 𝜙, 𝑤 𝑚𝑎𝑥 = 𝑤 𝑝𝑜𝑟𝑒 = 𝜌 𝑤 𝜙/𝜌 𝑑 , 𝜌 𝑤 ≈ 1 𝑀𝑔/𝑚 3 . For 𝑤 𝑚𝑎𝑥 in bentonite, it is important to have water mass content measurements for saturation in the state of low suction.

The difference between 𝑤 𝑚𝑎𝑥 and 𝑤 𝑝𝑜𝑟𝑒 in the bentonite can be balanced by introducing a higher mass density of the water 𝜌 𝑤 + = 𝜌 𝑤 (𝑤 𝑚𝑎𝑥 /𝑤 𝑝𝑜𝑟𝑒 ). Fig. 3 shows both the values 𝜌 𝑤 + derived from the experiments on the FEBEX bentonite and 𝜌 𝑤 + calculated from the expression ( 16) for 𝑤 𝑚𝑎𝑥 .

In the case of experiments, the water density 𝜌 𝑤 + has been calculated from 𝑤 𝑚𝑎𝑥 , which was the final water mass content averaged over the whole tested samples. The correlation of both types of values provides additional evidence of the suitability of the expression of oversaturation given in [START_REF] Papafotiou | Predismantling THM modelling of the FEBEX in situ experiment[END_REF].

Fig. 3: Comparison of data derived from tests [START_REF] Villar | Influence of dry density and water content on the swelling of a compacted bentonite[END_REF] (diamond and square symbols) and data derived from w max (ρ d , T = 20°C) defined in ( 16) (the red line in the picture).

Permeabilities and vapour diffusion

We assume that the saturated permeability depends on the dry density (similarly like 𝑤 𝑚𝑎𝑥 ) and the following expressions ( 17) and ( 18) for the FEBEX bentonite were derived from the measured values [START_REF] Villar | FEBEX-DP Post-mortem THM/THG Analysis Report[END_REF],

𝑘 𝑠 = 𝑘 𝑠 (𝜌 𝑑 ) = 6.46 ⋅ 10 -17 ⋅  𝑑 -22.5 [m 2 ], (17)
𝑘 𝑠 = 𝑘 𝑠 (𝜌 𝑑 ) = 2.0 ⋅ 10 -17 ⋅  𝑑 -18.95 [m 2 ], (18)
see Fig. 4. The corresponding conductivities are computed through the standard formula 𝐾 𝑠 = 𝑘 𝑠 𝑔𝜌 𝑤 /𝜇, where 𝑔 is the gravitational acceleration and 𝜇 = 𝜇(𝑇) is the dynamic viscosity of water, which depends on the temperature 𝑇, see e.g. [START_REF] Bear | Modeling Groundwater Flow and Contaminant Transport[END_REF]. Note that 𝜇 at 𝑇 = 20° is used for obtaining the graphs in Fig. 4. The relative permeability depends on the effective saturation (3) and the following dependence is used for the FEBEX bentonite

𝑘 𝑟 = 𝑘 𝑟 (𝑆 𝑒 ) = 𝑆 𝑒 3 , (19)
The water vapour diffusion coefficient 𝐷 𝑣 , which appears in ( 5) and ( 6), should be determined experimentally, see e.g. [START_REF] Moldrup | Modeling diffusion and reaction in soils.7. Predicting gas and ion diffusivity in undisturbed and sieved soils[END_REF], [START_REF] Moldrup | Predicting the Gas Diffusion Coefficient in Repacked Soil: Water-Induced Linear Reduction Model[END_REF]. Generally, the expression for this coefficient has the form 𝐷 𝑣 = 2.16 • 10 -5 ( 𝑇 273. [START_REF] Gens | A full-scale in situ heating test for high-level nuclear waste disposal: Observations, analysis and interpretation[END_REF] )

1.8 𝐷 𝑣𝑟 , (20)
where the first part provides the diffusion coefficient through open air and 𝐷 𝑣𝑟 is the relative diffusion coefficient depending on the air-filled porous space, see [START_REF] Moldrup | Predicting the Gas Diffusion Coefficient in Repacked Soil: Water-Induced Linear Reduction Model[END_REF].

It can be assumed that 𝐷 𝑣𝑟 = 𝑐𝜙(1 -𝑆 𝑤 ), where 𝑐 is a proportionality factor and the expression 𝜙(1 -𝑆 𝑤 ) stands for approximation of the volume fraction Φ 𝑣 , where the diffusion can occur. In this way, we shall consider

𝐷 𝑣𝑟 = 𝜏𝜙(1 -𝑆 𝑤 ), ( 21
)
where 𝜏 is the tortuosity. The first choice is the value 𝜏 = 0.8, which is used for the FEBEX bentonite e.g. in [START_REF] Rutquist | LBNL Modelling of Task A-1 using ROCMAS[END_REF]. Another choice is the Penman-Millington-Quirk (PMQ) model

𝐷 𝑣𝑟 = 0.66 𝜙(1 -𝑆 𝑤 ) 2 , (22)
see [START_REF] Moldrup | Modeling diffusion and reaction in soils.7. Predicting gas and ion diffusivity in undisturbed and sieved soils[END_REF]. It can also be interpreted as [START_REF] Wang | Non-isothermal flow in low permeable porous media: a comparison of 'Richards' and two-phase flow approaches[END_REF] with 𝜏 = 0.66(1 -𝑆 𝑤 ).

The above choices of 𝜏 are used in the definition of Model 1 and Model 2 in Section 6.2, and the application of these models in the modelling of FEBEX in situ test is discussed later in Section 7.

Heat

The heat conduction coefficient significantly depends on water saturation, i.e.

𝜆 = 𝜆(𝑆 𝑤 ).

The available data, see Fig. 5, show the dependency of 𝜆 on 𝑆 𝑟 , which is the relative saturation related to the bulk pore space defined in [START_REF] Villar | Report on Thermo-Hydro-Mechanical Laboratory Tests Performed by CIEMAT on Febex Bentonite 2004-2008[END_REF], see [START_REF] Alonso | The FEBEX benchmark test: Case definition and comparison of modelling approaches[END_REF]. The dependence can be approximated by the expression 𝜆(𝑆 𝑟 ) = 1.28 -0.71/(1 + 𝑒 10⋅(𝑆 𝑟 -0.65) ).

The dependence 𝜆 = 𝜆(𝑆 𝑤 ) is then obtained from [START_REF] Duncan | Modern techniques for the analysis of engineering problems in soft clay[END_REF] and the relation between 𝑆 𝑟 and 𝑆 𝑤 given in [START_REF] Villar | Report on Thermo-Hydro-Mechanical Laboratory Tests Performed by CIEMAT on Febex Bentonite 2004-2008[END_REF]. The evaporation or condensation heat, see [START_REF] Abed | A study on how to couple thermo-hydro-mechanical behaviour of unsaturated soils: Physical equations, numerical implementation and examples[END_REF], is neglected in the developed THM model.

Fig. 5: Experimental data for FEBEX and S2 bentonites from [START_REF] Alonso | The FEBEX benchmark test: Case definition and comparison of modelling approaches[END_REF]. FEBEX correlation (solid black line) corresponds to the formula (23).

Nonlinear elasticity

The elastic modulus is obtained from the oedometric experiments showing dependency of vertical (in this case also volumetric) strain 𝜀 𝑣𝑒𝑟𝑡 to the uniaxial (vertical) pressure (stress) 𝜎 𝑣𝑒𝑟𝑡 , see Fig. 6. These oedometric tests can be simulated with hydro-mechanical model using the described relations plus information about the mechanical behaviour of bentonite. If isotropic behaviour is assumed and the Poisson ratio is estimated as 𝜈 = 0.3, then we need to consider elastic modulus 𝐸 depending on dry density, deformations and level of saturation.

With respect to the linear dependence of strains on the logarithm of the load, which is typical for soils, the elastic modulus for the saturated state is searched in the form of exponential dependence on dry density and strains. Because the dry density also depends on strains, the searched modulus can be expressed as depending only on dry density.

The inverse analysis of the data from oedometer tests, see Fig. 6, provides the relation

𝐸 𝑠𝑎𝑡 = 𝐸 𝑠𝑎𝑡 (𝜌 𝑑 ) = exp(8.2652 𝑑 -10.62) [MPa]. (24)
The above expression depends on the compaction expressed by the dry density 𝜌 𝑑 [Mg/m 3 ].

The dependence of elastic modulus on the level of saturation or equivalently on suction is studied e.g. in [START_REF] Man | Thermal, Hydraulic and Mechanical Properties of Sealing Materials[END_REF], and in the developed model it is simplified to the relation

𝐸 = 𝐸(𝑠, 𝜌 𝑑 ) = { 𝐸 𝑠𝑎𝑡 if 𝑠 ≤ 0, 𝐸 𝑠𝑎𝑡 + 𝑠 50 (𝐸 𝑚𝑎𝑥 -𝐸 𝑠𝑎𝑡 ) if 0 < 𝑠 ≤ 50 MPa, 𝐸 𝑚𝑎𝑥 if 𝑠 > 50 MPa, (25)
where 𝑠 denotes suction in MPa, 𝐸 𝑠𝑎𝑡 is defined in [START_REF] Alonso | A constitutive model for partially saturated soils[END_REF] and 𝐸 𝑚𝑎𝑥 = 100 MPa. This choice reflects the application in the analysis of FEBEX experiment with a significant number of gaps among blocks and the corresponding decrease of the dry density to 𝜌 𝑑 = 1.6 Mg/m 3 (in FEBEX the dry density of isolated block is 𝜌 𝑑 = 1.7 Mg/m 3 ).

Swelling

The effect of swelling is again a consequence of the inner microstructure of the bentonite. We use the expression

𝜎 𝑠𝑤𝑒𝑙𝑙 = 𝜎 𝑠𝑤𝑒𝑙𝑙 (𝑤, 𝜌 𝑑 ) = 𝑤 -𝑤 𝑖𝑛𝑖𝑡 𝑤 𝑚𝑎𝑥 -𝑤 𝑖𝑛𝑖𝑡 𝜎 𝑠𝑤𝑒𝑙𝑙,𝑚𝑎𝑥 𝐼, (26)
where 𝑤 𝑚𝑎𝑥 is taken from ( 16), 𝑤 𝑖𝑛𝑖𝑡 is the water mass content in the time when we start to consider swelling stress and 𝜎 𝑠𝑤𝑒𝑙𝑙,𝑚𝑎𝑥 depends on the dry density 𝜌 𝑑 expressed in Mg/m 3 and was determined from the fitting of experimental data from [START_REF] Alberdi | FEBEX project: full-scale engineered barriers experiment for a deep geological repository for high level radioactive waste in crystalline host rock[END_REF] as

𝜎 𝑠𝑤𝑒𝑙𝑙,𝑚𝑎𝑥 = 𝜎 𝑠𝑤𝑒𝑙𝑙,𝑚𝑎𝑥 (𝜌 𝑑 ) = exp(6.77𝜌 𝑑 -9.07) [𝑀𝑃𝑎]. (27)

THM induced changes of the densities and the pore space

Let us assume that the solid grains, which occupy the volume 𝑉 𝑠 , are relatively stiff so that their mechanical deformation is negligible. On the other hand, the thermal expansion of the solid is not neglected. Then the mass density 𝜌 𝑠 = 𝜌 𝑠 (𝑇) and the volumetric thermal expansion coefficient 𝛼 𝑣 = 𝛼 𝑣𝑠 can be expressed as

𝛼 𝑣𝑠 = 1 𝑉 𝑠 𝜕𝑉 𝑠 𝜕𝑇 = - 1 𝜌 𝑠 𝜕𝜌 𝑠 𝜕𝑇 .
Integration of this differential equation with the initial condition 𝜌 𝑠 (𝑇 𝑟𝑒𝑓 ) = 𝜌 𝑠0 then provides

𝜌 𝑠 = 𝜌 𝑠 (𝑇) = 𝜌 𝑠0 𝑒 -𝛼 𝑣𝑠 (𝑇-𝑇 𝑟𝑒𝑓 ) . (28)
Similarly, assuming incompressibility of water,

𝜌 𝑤 = 𝜌 𝑤 (𝑇) = 𝜌 𝑤0 𝑒 -𝛼 𝑣𝑤 (𝑇-𝑇 𝑟𝑒𝑓 ) , (29)
but here 𝛼 𝑣𝑤 still depends on the temperature.

Under the assumption of stiff solid grains and correspondingly the Biot-Willis coefficient 𝛼 = 1, the change of porosity 𝜙 can be expressed by the following differential equation derived from the mass balance of the solid component under the small strain assumption

𝜕𝜙 𝜕𝑡 = (1 -𝜙) ( 𝜕𝜀 𝑣 𝜕𝑡 -𝛼 𝑣𝑠 𝜕𝑇 𝜕𝑡 ) = (1 -𝜙) 𝜕𝜀 𝑒𝑙,𝑣 𝜕𝑡 , (30)
where 𝜀 𝑣 = 𝑡𝑟(𝜀) and 𝜀 𝑒𝑙,𝑣 = 𝑡𝑟(𝜀 𝑒𝑙 ). For more details, see e.g. [START_REF] Abed | A study on how to couple thermo-hydro-mechanical behaviour of unsaturated soils: Physical equations, numerical implementation and examples[END_REF].

The solution of the equation ( 30) has the form

𝜙(𝑡) = 1 -(1 -𝜙 0 ) exp(𝜀 𝑒𝑙,𝑣,0 -𝜀 𝑒𝑙,𝑣 (𝑡)), (31)
where 𝜙 0 and 𝜀 𝑒𝑙,𝑣,0 describe an initial state. The use of a simple Padé approximation 𝑒 𝑥 ≈ 1 1-𝑥 provides an approximate expression used in the implementation 𝜙 ≈ 𝜙 0 + 𝜀 𝑒𝑙,𝑣 -𝜀 𝑒𝑙,𝑣,0 1 + 𝜀 𝑒𝑙,𝑣 -𝜀 𝑒𝑙,𝑣,0 .

Finally, note that 𝜌 𝑑 = 𝜌 𝑠 (1 -𝜙) depends on temperature and the volumetric strain 𝜀 𝑒𝑙,𝑣 .

Numerical solution and implementation

The model is implemented via the finite element software COMSOL Multiphysics [START_REF]COMSOL Multiphysics® v. 5.3a. www.comsol[END_REF]. The space discretisation is done through weak formulation and the finite element (FEM) method. A locally refined mesh of finite elements is used for all THM processes, see Fig. 8. Finite elements with quadratic polynomials are used for discretisation of heat and fluid flow, finite elements with cubic polynomials are used for the analysis of stresses and strains. The time discretisation uses the backward Euler method, which leads to a time-stepping procedure with the requirement of solving complex nonlinear system in each time step. This system arises from the finite element discretisation of all coupled THM processes, i.e. thermal-hydraulic and mechanical parts.

COMSOL Multiphysics offers two tools for solving the system, monolithic (full coupled) and segregated solvers as well as their combinations, see [START_REF]COMSOL Multiphysics® v. 5.3a. www.comsol[END_REF]. So, the segregation is used between mechanical part and thermal-hydraulic part, and a monolithic damped Newton solver treats the thermal-hydraulic part. The adaptive time step control watches both discretisation error and convergence of the solver for the systems arising within the time steps.

COMSOL Multiphysics allows to abstract from details of the implementation of numerical methods and considers it as a black box. Thus, we can only mention that the default accuracy was used during the computations. Some difficulties with the convergence arose in the case of inconsistent boundary and initial conditions which occur e.g. during the transition from one modelling phase of the FEBEX experiment to another one with possibly changed geometry. To overcome such convergence problems very short artificial transition phases were used.

The FEBEX in situ test

The full-scale engineered barriers experiment (FEBEX) was a research and demonstration project that was initiated by the "Spanish Agency for Radioactive Waste Disposal" (ENRESA). It followed from the Spanish concept for the isolation of high-level radioactive waste. The subject of study in this article is the FEBEX in situ test, which was located in a crystalline host rock of the Aare Massif in central Switzerland in the Grimsel Test Site [START_REF] Lloret | FEBEX II Project Final report on thermo-hydro-mechanical laboratory tests[END_REF].

The primary goal of the FEBEX test was to study the behaviour of near-field components in a repository for high-level radioactive waste in granite formations, demonstration of the feasibility of constructing the engineered barrier system in a horizontal configuration, and analysis of the technical problems of this type of disposal method. The secondary aim was to develop modelling tools for the thermo-hydro-mechanical (THM) and thermo-hydrogeochemical (THG) processes and provide validation of models against experimental results [START_REF] Lloret | FEBEX II Project Final report on thermo-hydro-mechanical laboratory tests[END_REF].

The FEBEX in situ test involves the installation of two cylindrical heaters (4.54 m long and 0.9 m in diameter) disposed in a cylindrical steel liner with holes (diameter 0.97 m) centred in a tunnel (diameter 2.27 m) excavated by TBM (a tunnel boring machine). The rest of the space was filled by a barrier (thickness 0.64 m) made of highly compacted unsaturated blocks of the FEBEX bentonite. The backfilled area was sealed with a plain concrete plug placed into a recess excavated in the rock and having a length of 2.70 m (Fig. 7).

A total of 632 instruments were placed in the host rock and buffer to monitor relevant values of temperature, humidity, total and pore pressure, and displacements. A more detailed description of the spatial dimensions, bentonite placement, sensor placement, and experiment installation is provided in [START_REF] Fuentes-Cantillana | FEBEX: Full-scale engineered barriers experiment in crytalline host rock: Final design and installation of the 'in situ' test at Grimsel[END_REF].

Fig. 7: The FEBEX experiment.

The FEBEX project started in 1994. One year later started the FEBEX in situ test with rock-mass characterisation and tunnel excavation. The construction of the concrete plug was completed in October 1996, and the heating began on February 27, 1997. The power of heaters was first increased stepwise and manually to reach a targeted maximum heater temperature of 100 °C for 53 days. After that, the system was switched to automatic mode for keeping the constant temperature at 100 °C in one measurement point in the middle of each heater [START_REF] Lloret | FEBEX II Project Final report on thermo-hydro-mechanical laboratory tests[END_REF].

Partial dismantling of the concrete plug and buffer around the heater closer to the access tunnel was carried out during the summer of 2002, after five years of continuous heating. The removed area was up to a distance of 2 meters from the second heater to minimise disturbance of the nondismantled area. During this dismantling process, the second heater stayed operational. Afterwards, additional sensors were installed and a steel cylinder with a length of 1 m was inserted in the void left by the removed heater in the centre of the buffer, and the experiment was again sealed. All details of the dismantling process and measurements of dismantled data are described in [START_REF] Bárcena | Dismantling of heater No. 1 at the FEBEX in situ test. Description of operations[END_REF].

The second heater stayed operational for 18.4 years and the rest of the FEBEX test was dismantled in 2015, including the second heater. A number of parts of the FEBEX test were studied, like bentonite, rock, relevant interfaces, sensors, metallic components, and tracers. Detail description of the second dismantling process is in [START_REF] García-Siñeriz | FEBEX DP: Dismantling of heater 2 at the FEBEX "in situ" test. Description of operations[END_REF]. The monitoring data from the FEBEX in situ test are evaluated in [START_REF] Martínez | FEBEXe: Final Sensor Data Report (FEBEX "in situ" Experiment)[END_REF], and the final report from the dismantling of the FEBEX test is provided in [START_REF] García-Siñeriz | FEBEX DP: Dismantling of heater 2 at the FEBEX "in situ" test. Description of operations[END_REF].

Numerical model of the FEBEX experiment

Due to the horizontal deposition of heaters and axisymmetric geometry of the experiment, the model of the FEBEX experiment is treated as 2D-axisymmetric with axis placed in the middle of the circular tunnel with heaters. This axisymmetric simplification is natural for the saturation and heat flow and acceptable for the mechanics with swelling dominating over gravitational forces and the effect of possibly anisotropic initial stress. The cross-section shown in Fig. 8 has 50 m in radial and 120 m in axial directions, see Figs. 8 and9. These dimensions are supposed as sufficient; the model does not show a change of temperature against initial temperature for 18.4 years at the outer boundary of the computational domain. The finite element mesh is depicted in Fig. 8. Quadratic finite elements are used for space discretisation of heat and water flow, and cubic elements are used for the mechanical part. The mesh for space discretisation includes 8 659 triangular and 3 736 quadrilateral elements, i.e., in total maximum of 12 395 elements. The number of elements depends on the modelling phases, see Section 6.2.

Fig. 8: Finite element meshes -the whole computational domains and the details.

The computational domain is divided into five subdomains, see Fig. 9, where each subdomain represents the same type of material. Most of the material properties and initial conditions used in the developed THM models are summarised in Table 1. In more details, the initial conditions for each subdomain are as follows.

Crystalline host rock has initial water pressure 𝑝 𝑤 = 0.7 MPa and isotropic horizontal initial Labelled notes to the data in Table 1:

(a) Average value considering the dry density of blocks and the presence of gaps. (b) According to [START_REF] Schneider | Determining hydraulic properties of concrete and mortar by inverse modelling[END_REF]. (c) Maximal water content fully determined by the (bulk) pore space. (d) Retention function from [START_REF] Gens | A full-scale in situ heating test for high-level nuclear waste disposal: Observations, analysis and interpretation[END_REF] for Grimsel granite. (e) Retention function from [START_REF] Gens | A full-scale in situ heating test for high-level nuclear waste disposal: Observations, analysis and interpretation[END_REF] for bentonite with dry density 1.58-1.59 Mg/m 3 .

Further comments on the data in Table 1:

THM model described in Sections 2 and 3 is considered for all materials with the exception of the heater. For the heater, a thermo-mechanical behaviour is assumed. A simplification is used for rock and concrete, where a classical retention function depending only on suction and without the closing term is used. All other parameters are kept constant. The liner is supposed to have thermo-mechanical properties as steel, but due to the holes, the flow through liner is possible. Therefore, THM model is also used for the liner with modified parameters.

Computational phases

The FEBEX experiment evolves in time, and consequently, different phases must be considered in the modelling. FEBEX experiment formally started 25/09/1995 by tunnel excavation and ended 20/07/2015 by the end of the final dismantling. The computational phases are described in detail in [START_REF] Lloret | FEBEX II Project Final report on thermo-hydro-mechanical laboratory tests[END_REF] as well as in [START_REF] Gens | DECOVALEX-2019[END_REF] and [START_REF] Gens | HM and THM interactions in bentonite engineered barriers for nuclear waste disposal[END_REF]. It is essential to realise that the computational domain changes in time, e.g. by removing a part of the domain after the first dismantling or by the construction of the second plug. Phase 10 (14 days) is the period of cooling the FEBEX test after turning off the second heater.

In this phase zero heat flow is prescribed on the whole boundary of the second heater.

Phase 11 (114 days) starts with removing of the second plug [START_REF] García-Siñeriz | FEBEX DP: Dismantling of heater 2 at the FEBEX "in situ" test. Description of operations[END_REF]. Data from the second dismantling process [START_REF] Villar | FEBEX-DP Post-mortem THM/THG Analysis Report[END_REF] are obtained at the end of this phase. A comparison of the model outputs against the dismantling data is provided in the next section, see also [START_REF] Gens | DECOVALEX-2019[END_REF] and [START_REF] Gens | HM and THM interactions in bentonite engineered barriers for nuclear waste disposal[END_REF].

The simulation of the FEBEX in situ test included all eleven phases. The modelling results consider mainly two stages, the first one considers outputs after performing Phases 1 to 8, the second one after Phases 8 to 11. The detailed description of all settings is in [START_REF] Gens | DECOVALEX-2019[END_REF] (Annex F).

Variants of the THM model for bentonite

The modelling of the processes in FEBEX bentonite was performed by two THM models which use the relations and data described in Sections 2 and 3, and parameters described in Table 1. The models differ only slightly in the adopted expressions for saturated permeability. The main difference is in the choice of expression for vapour diffusion. As we shall see in the next section, this difference in model setting can lead to a more significant difference in some outputs. The differences in the variants of the model are just the following.

Model 1 uses the saturated permeability [START_REF] Lewis | The Finite Element Method in the Static and Dynamic Deformation and Consolidation of Porous Media[END_REF] and the relative vapour diffusion coefficient using the formula [START_REF] Wang | Non-isothermal flow in low permeable porous media: a comparison of 'Richards' and two-phase flow approaches[END_REF] with the tortuosity 𝜏 = 0.8. Model 2 uses the saturated permeability [START_REF] Bear | Modeling Groundwater Flow and Contaminant Transport[END_REF] and the relative vapour diffusion coefficient using the formula [START_REF] Philip | Moisture Movement in Porous Materials under Temperature Gradients[END_REF].

Note that the outputs from Model 1 are similar as outputs presented in [START_REF] Gens | A full-scale in situ heating test for high-level nuclear waste disposal: Observations, analysis and interpretation[END_REF] for the first fiveyear period, where outputs from this model fit very well the measured data from the first dismantling (dry density and water content). Model 1 was consequently used for modelling the saturation up to the final dismantling. But then it provided a significantly lower value of the mass water content than shown by the experimental data, see Fig. 13. Similar observations are presented in a blind study provided within [START_REF] Papafotiou | Predismantling THM modelling of the FEBEX in situ experiment[END_REF] with modelling performed by the CIMNE (International Center for Numerical Methods in Engineering, Spain) using the CodeBright software.

This inconsistency between the model outputs and measured data in the water content motivated searching for possible reason and change of Model 1. As the model results provide a good agreement with in situ data after first dismantling, a simple change of parameters was not acceptable.

A suitable change was found in the modification of the vapour diffusion coefficient, which played a significant role in the saturation, especially in locations close to heaters. Therefore, the Penman-Millington -Quirk (PMQ) model, see [START_REF] Moldrup | Modeling diffusion and reaction in soils.7. Predicting gas and ion diffusivity in undisturbed and sieved soils[END_REF], was adopted and used in Model 2. The hydraulic permeability should be slightly decreased to fit the data after the first dismantling. A further discussion concerning the introduced models can be found in the next section.

Comparison of computed results with experimental data

This section describes the simulation of the THM processes in the bentonite barrier during the operation of the FEBEX in situ test. The simulation follows all phases described in the previous section and uses the mathematical model described in Sections 2 and 3 with two model variants (Model 1 and

2) which differ mainly in the adopted vapour diffusion term. Some outputs were selected which characterise the fitting of computed values with experimental data in different space positions (different cross-sections and different distances from heaters) with time evolution to 5 years (about 1800 days) and 18 years (about 6600 days). The outputs also characterise sensitivity to some relations included in the model. Note that more outputs are presented in [START_REF] Gens | DECOVALEX-2019[END_REF].

Fig. 12 shows the comparison of experimental data from the monitoring of time evolution of temperatures, relative humidity and total radial and axial stress against computed results using both models 1 and 2. The main difference between the models 1 and 2 can be seen in the time evolution of the humidity. Both models fit this evolution for five years. For the evolution up to 18 years, there is less high-quality experimental data, especially for the relative humidity close to the remaining heater.

Thus, it is more difficult to evaluate the fit of the models, but the difference between the results provided by the models 1 and 2 is significant. A better comparison of the models is possible if the data from dismantling are considered, see Fig. 13.

Note that the measurement data in Fig. 12 also show some difference between measurements from different sensors in the same cross-sections. These differences show some heterogeneity and the difficulty of the measurements. As concerns the fit of radial stresses, note that real measurements were provided by pressure sensors placed into gaps between bentonite blocks and therefore they indicate real stress only after some time necessary for filling the gaps by swelling. Further, there was an asymmetric gap between the rock and bentonite at the top of the tunnel. This gap also had to be closed prior to accurate measurements. These construction facts, as well as asymmetric gravity forces, are not included in the model. Obviously, the model then gives an immediate increase of stresses. Also, the use of elasto-plastic behaviour for the bentonite could increase the accuracy of the model.

The comparison of computed data with data from dismantling after 5 and 18.4 years of experiment operation is depicted in Fig. 13. This figure focuses on the areas of two cross-sections passing through the middle of heaters with high-temperature gradients, where the thermal vapour diffusion plays a significant role. The comparison shows that both Models 1 and 2 are comparable and provide a good fit for dismantling after five years, but the Model 2 is better in fitting the experimental data from dismantling after 18.4 years, especially in water content and degree of saturation.

The obtained results suggest a question of how further optimise the relations and parameters in THM model. One possibility is e.g. 3D model; another possibility is the introduction of heterogeneity, namely the use of different parameter settings in different parts of the computational domain. This requires to get more input information, but otherwise, the model will not be able to be accurate in some aspects. This is illustrated below.

The graph in Fig. 14 shows the power of both heaters recorded during the experiment operation [START_REF] Gens | Additional data for FEBEX experiment which were kindly provided on request[END_REF]. In the first five-year period both heaters were active and the power of both was controlled by temperature sensors to be 100°C in the middle of the heater length. The available information about the FEBEX test does not indicate any source of difference in the heat power between both heaters as they should be surrounded by the same arrangement of bentonite blocks and surrounded by the same granitic rock, the water supply at the bentonite-granite interface should be sufficient, see [START_REF] Gens | A full-scale in situ heating test for high-level nuclear waste disposal: Observations, analysis and interpretation[END_REF].

However, in the experiment, there is about a 10% difference in the power of heaters. The power is lower for heater #1 located in the rear part of the FEBEX experiment, which was dismantled first, after five years of operation. It indicates that some difference (heterogeneity) in physical reality should exist. The different heat power of heaters can be explained by different heat dissipation due to higher water saturation around the heaters. Such a situation could be caused by higher permeability of the bentonite block around heater #2 or slower closing of technical gaps between the bentonite blocks and therefore better water conduction and faster bentonite saturation than in the vicinity of the first heater. Logically, without introducing some heterogeneity, the model is not able to explain this effect. Comparison with data from the first partial dismantling (5 years) and the complete dismantling (18.4 years) of the experiment. The description of the saturation uses the definition related to bulk pores [START_REF] Villar | Report on Thermo-Hydro-Mechanical Laboratory Tests Performed by CIEMAT on Febex Bentonite 2004-2008[END_REF] and corresponds to the state after dismantling. Consequently, the saturation goes over 100% (oversaturation occurs). Fig. 14: Measured and computed power of the heaters (Heater 1 is closer to the access tunnel), see [START_REF] Gens | Additional data for FEBEX experiment which were kindly provided on request[END_REF].

Conclusions

This paper formulates a model and its variants suitable for simulation of the THM processes occurring in the bentonite barriers around canisters with the spent nuclear fuel in a deep geological repository and similar applications. The main focus is on the processes in bentonite including the special character of bentonite behaviour with oversaturation and swelling.

The introduced THM model includes (1) flow in a variably saturated bentonite described by Richards type equation with vapour diffusion, (2) nonlinear elasticity and (3) heat propagation by conduction. A novelty of the model is the introduction of saturation which is not related to filling the pores but described by using the ratio of the current and maximal mass water content. This definition allows to include the oversaturation effect, which was measured in FEBEX in situ test. Subsequently, the model uses the water retention function providing the relation between the mass water content on one side and the suction, dry density and temperature on the second side. Such dependence is approximated by a modified multivariable van Genuchten function which includes the oversaturation term and fits reasonably well the provided measurement data. Special care is also devoted to vapour diffusion, which has a significant role in the model, as it was shown in the previous section.

Implementation of the model with proper discretisation in space and time and solvers for coupled multiphysics, including nonlinear processes, is a challenging task that was successfully done by using the tools provided by COMSOL Multiphysics software. Still, some effort in further development will be required in the future to fully control the discretisation error and enable the efficient solution to problems in 3D space (in the modelling of the FEBEX experiment with horizontal deposition of the heaters, it was natural to use the 2D axisymmetric simplification of the situation).

2 .

 2 Note that the volumetric parts of the strain tensors are defined through the trace, 𝜀 𝑣 = 𝑡𝑟(𝜀) = ∑ 𝜀 𝑖𝑖 𝑖 , 𝜀 𝑣 = 𝑡𝑟(𝜀 𝑒𝑙 ) + 𝑡𝑟(𝜀 𝑡ℎ ) = 𝜀 𝑒𝑙,𝑣 + 𝜀 𝑡ℎ,𝑣 , 𝜀 𝑡ℎ,𝑣 = ∑ 𝜀 𝑡ℎ,𝑖𝑖 𝑖 = 𝛼 𝑣 (𝑇 -𝑇 𝑟𝑒𝑓 ).

Fig. 2 :

 2 Fig. 2: Retention curves (bold and dashed lines) for different dry densities ρ d = 1.55; 1.6, ; 1.65; 1.7 Mg/m 3 (red, green, yellow and blue curves, respectively) and temperatures T = 20, 60 °C and a comparison with the measured data (isolated symbols) from [13].

Fig. 4 :

 4 Fig. 4: Comparison of the hydraulic conductivities derived from the expressions (17) and (18) (see dashed and solid red lines) with measured data (isolated symbols).

Fig. 6 :

 6 Fig. 6: Strains with elastic component for different vertical loads applied to fully saturated blocks with different values of initial dry density in g/cm 3 = Mg/m 3 , see [7].

  The thermal and hydraulic parts of the model are implemented in PDE (partial differential equation) module. It was necessary due to complicated links including thermal vapour term in the Richards equation. For axisymmetric model, the transformation into cylindrical coordinates should be additionally implemented as the PDE module directly offers only Cartesian coordinates. The mechanical part is implemented via the Structural mechanics module. Two of the basic equations describing THM processes, namely heat conduction and Richards equation, are treated as strongly coupled and involve time derivative of the basic unknowns. The momentum balance equation does not involve the time derivative of the displacement (no inertia forces); displacements just instantly react to changes of loading.

stress 28

 28 MPa. Bentonite blocks have initial dry density 𝜌 𝑑 = 1.6 Mg/m 3 . This value corresponds to the average between compacted blocks of FEBEX bentonite with initial 𝜌 𝑑 = 1.7 Mg/m 3 and technological gaps. The initial water content of bentonite blocks is 𝑤 𝑖𝑛𝑖𝑡 = 0.14. The corresponding suction is 135 MPa. Concrete plugs (original and second one constructed after the first dismantling) have initial suction 10 MPa. The initial temperature in all subdomains is 12 °C.

Fig. 9 :

 9 Fig. 9: Computational domains. 1. Crystalline host rock, 2. Bentonite, 3a. Original concrete plug, 3b. New concrete plug, 4. Heater, 5. Liner.

Fig. 10 Phase 1 and 2 -Phase 3 , 4 , 5 , 6 -Phase 7 -Phase 8 -Fig. 11 :

 10234567811 Fig. 10 Computational phases with different degrees of freedom (DoF) -active DoFs are inside the coloured parts of the model.

Fig. 12 :

 12 Fig. 12: Comparison of time-dependent evolution up to 5 and 18 years of measured and calculated values of temperature, relative humidity and total radial and axial stress in selected cross-sections and distances from the heaters. Circles denote measurements by several sensors in each cross-section, dashed and solid lines show output from Model 1 and Model 2, respectively. Colours indicate different distances from the tunnel axis. Note that positive stress in this figure means compression.

Fig. 13 :

 13 Fig. 13: Comparison of measured and calculated values of dry density, water mass content and degree of saturation in the area with a large temperature gradient (sections 27 and 49 through the heaters).

  conservation equations for separate liquid water and vapour flow. If 𝑀 𝑤 and 𝑀 𝑣 denote the mass of water and vapour in a volume 𝑉, then the changes of 𝑀 𝑤 and 𝑀 𝑣 in time 𝑡 can be expressed as where 𝜌 𝑤 is the mass density of water, 𝑣 𝑤 is the Darcy velocity described later in (2), 𝑞 𝑣 is the diffusive mass flux of the vapour, see[START_REF] Stephansson | Coupled Thermo-Hydro-Mechanical-Chemical Processes in Geo-Systems: Fundamentals, Modelling, Experiments and Applications[END_REF], and 𝑄 𝑤 is the mass source term for water. In the considered model, 𝑄 𝑤 is determined by the deformation and the level of water saturation of volume V, see Section 2.3.The term 𝐸 denotes possible mass exchange between the liquid water and vapour, which in the sum disappears, see the classical paper[START_REF] Vries | Simultaneous transfer of heat and moisture in porous media[END_REF]. Moreover,

	𝜕𝑀 𝑤 𝜕𝑡	𝑉 = ∫ [-∇ ⋅ (𝜌 𝑤 𝑣 𝑤 ) -𝐸 + 𝑄 𝑤 ]	d𝑥,
		𝜕𝑀 𝑣 𝜕𝑡	𝑉 = ∫ [-∇ ⋅ 𝑞 𝑣 + 𝐸]	d𝑥,
	𝜕(𝑀 𝑤 + 𝑀 𝑣 ) 𝜕𝑡	=	𝜕 𝜕𝑡	∫ (𝑤𝜌 𝑑 𝑉	+ 𝑚 𝑣 ) d𝑥 ≈	𝜕 𝜕𝑡	𝑉 ∫ 𝑤𝜌 𝑑	d𝑥,

where 𝑤 is the mass water content, 𝜌 𝑑 is the dry density and 𝑤𝜌 𝑑 expresses (density of) mass of the water. If 𝑀 𝑤 and 𝑀 𝑠 denote the mass of water and solid in the volume 𝑉, then 𝑤 = 𝑀 𝑤 /𝑀 𝑠 , 𝜌 𝑑 = 𝑀 𝑠 /𝑉 and 𝑤𝜌 𝑑 = 𝑀 𝑤 /𝑉. The term 𝑚 𝑣 = 𝜌 𝑣 Φ 𝑣 stands for the mass of vapour, which can be expressed through 𝜌 𝑣 and Φ 𝑣 being the mass density and volume fraction for vapour, respectively. The mass density 𝜌 𝑣 can be expressed from the Kelvin-Laplace equation, see later in this subsection; the volume fraction Φ 𝑣 plays a role in the vapour diffusion, see Section 3.3. The term 𝑚 𝑣 can regularise the equation for 𝑤 close to zero, which does not occur in the intended applications. Our simulations include the term 𝑚 𝑣 , but by numerical tests, it was found that neglecting this term just very slightly changes the solution.

Table 1 :

 1

		bentonite	rock	concrete	heater TM	liner
	Initial dry density 𝜌 𝑑 [Mg/m 3 ]	1.60 (a)	2.64	2.2	7.8	7.8
	Initial porosity 𝜙 [-]	0.41	0.016	0.11 (b)	-	0.1
	Initial pressure 𝑝 𝑤 [MPa]	-135	0.7	-10	-	-135
	Saturated permeability 𝑘 𝑠 [m 2 ]	(17), (18)	5 ⋅ 10 -18 5.7 ⋅ 10 -20 (b)	-	4 ⋅ 10 -21
	Relative permeability 𝑘 𝑟 = 𝑆 𝑒 𝑛 ,	n=3	n=1	n=8	-	n=1
	𝑤 𝑚𝑎𝑥 [kg/kgbentonite]	(16)	=𝜙/𝜌 𝑑 (c)	=𝜙/𝜌 𝑑 (c)	-	=𝜙/𝜌 𝑑
	𝑤 𝑟𝑒𝑠 [kg/kgbentonite]	0.001	0.0001	0.001	-	(c) 0.001
	𝜆 1 van Genuchten [-]	0.14	0.595 (d)	0.404 (b)	-	0.17 (e)
	𝑃 0 [MPa]	(16)	1.74 (d)	13.1 (b)	-	4.5 (e)
	𝜆 0 van Genuchten [-]	1.5	-	-	-	-
	𝑃 00 [MPa]	1000	-	-	-	-
	Elastic modulus E [GPa]	(24), (25)	60	30	210	210
	Poisson ratio 𝜈 [-]	0.3	0.25	0.3	0.3	0.3
	Thermal conduct. 𝜆 [Wm -1 K -1 ]	(23)	3.6	1.7	46	46
	Linear thermal expansion 𝛼 𝑣𝑠 [-]	10 -5	7.8 ⋅ 10 -6	10 -5	10 -5	10 -5

Material properties. Identifiers (a)-(e) refer to notes given below the table.

Acknowledgement:

DECOVALEX is an international research project comprising participants from industry, government and academia, focusing on development of understanding, models and codes in complex coupled problems in sub-surface geological and engineering applications; DECOVALEX-2019 is the current phase of the project.

The authors appreciate and thank the DECOVALEX-2019 Funding Organisations Andra, BGR/UFZ, CNSC, US DOE, ENSI, JAEA, IRSN, KAERI, NWMO, RWM, SÚRAO, SSM and Taipower for their financial and technical support of the work described in this paper. The statements made in the paper are, however, solely those of the authors and do not necessarily reflect those of the Funding Organisations.

The introduced THM model was used for the simulation of THM processes occurring in the FEBEX in situ test. We can finally conclude the following.

The performed simulations used high-quality data available for the FEBEX bentonite behaviour obtained from careful laboratory testing. The availability of such data is crucial for the model application.

It was demonstrated that the introduced THM model and its implementation can provide the simulation of the processes occurring in FEBEX in situ test with the results which are in a good agreement with the experimental data or trends visible from the experimental data. In such a way, the FEBEX in situ test provided a validation of the model and its implementation.

The reported results show the significance of the role of some relations and parameters in the model. We can mention the expression of saturation and retention function through water mass content and maximum water content higher than just filling of bulk pores (oversaturation); modelling THM processes with specific parameters, which are evolving in time and depend on the current state of all processes; the significance of finding a proper expression for the vapour diffusion, which is suitable for the regions close to heat sources and simulation over longer time periods.

It was shown that the application of the model could encounter some uncertainties in the description of the situation, e.g. uncertainties induced by the heterogeneity of the geological environment and the barrier itself (see the discussion about the different power of the heaters in the previous section). Such uncertainty should be studied and estimated by simulations under different conditions and heterogeneities of the environment.

The FEBEX in situ test provides unique and high-quality data, which enables the development of proper mathematical models involving coupled THM processes and validation of such models against the experimental data. It also enables comparison of different approaches to construction and implementation of the models, see [START_REF] Gens | DECOVALEX-2019[END_REF] and [START_REF] Gens | HM and THM interactions in bentonite engineered barriers for nuclear waste disposal[END_REF]. This paper would like to contribute to the core idea of the DECOVALEX project which says that the combination of experiments and modelling is the way how to get the knowledge necessary for understanding processes which are crucial for the construction of a safe deep geological repository of the high-level radioactive waste.


R Ravaud

G Lemarquand
email: guy.lemarquand@univ-lemans.fr

MAGNETIC FIELD IN MRI YOKELESS DEVICES: ANALYTICAL APPROACH

This paper presents a three-dimensional analytical expressions for studying the static magnetic field produced by Magnetic Resonance Imaging structures. This medical imaging technique uses a very high and uniform magnetic field produced by ring permanent magnets with rotating polarizations. However, the manufacturing of such ring permanent magnets is difficult to realize. Consequently, such ring permanent magnets are replaced by assemblies of tile permanent magnets uniformly magnetized. Unfortunately, the magnetic field produced by these tile permanent magnets uniformly magnetized is both less important and less uniform than the one produced by an idealized ring permanent magnet. We propose in this paper to study the influence of the number of tile permanent magnets used on the magnetic field properties.

INTRODUCTION

Magnetic Resonance Imaging is a medical imaging technique widely used in radiology to visualize the internal structure of a body. This technique generally uses a very high static magnetic field produced by permanent magnets [START_REF] Marble | Strong, stray static magnetic fields[END_REF]- [START_REF] Chen | An improved discrete configuration of a cylinder magnet for portable nuclear magnetic resonance instruments[END_REF]. It is well known that the magnetic field homogeneity of imaging region is the most important criterion for optimizing MRI structures [START_REF] Ren | Study on shimming method for open permanent magnet of mri[END_REF] [START_REF] Shirai | Development of a compact mouse mri using a yokeless permanent magnet[END_REF]. On the other hand, the magnetic field must be the greatest [START_REF] Abele | A general method for flux confinement in permanent magnet structure[END_REF]- [START_REF] Li | Efficiency of permanent magnet assemblies for mri devices[END_REF]. Several approaches have been carried out for optimizing MRI structures so as to improve the magnetic field homogeneity as well as its magnitude [START_REF] Miyata | 3-d magnetic field analysis of permanent-magnet type of mri taking account of minor loop[END_REF] [START_REF] Podol'skii | Development of permanent magnet assembly for mri devices[END_REF]. Strictly speaking, these approaches are based on either analytical [START_REF] Marble | Designing static fields for unilateral magnetic resonance by a scalar potential approach[END_REF] or numerical approaches [START_REF] Ryu | Development of permanent magnet assembly for mri devices[END_REF]. For example, the finite-element method has been used for studying MRI devices made of permanent magnets [START_REF] Takahashi | 3-d fem analysis of residual magnetism produced by x-gradient coil of permanent magnet type of mri[END_REF]. Furthermore, some original methods have been proposed for designing MRI structures made of permanent magnets [START_REF] Sinha | Design concepts of optimized mri magnet[END_REF] by employing inverse methods [START_REF] Thompson | An inverse approach to the design of mri main magnets[END_REF]. All these studies deal with the optimization of the static field produced by permanent magnets or currents in coils. In the case of MRI made of permanent magnets, this static magnetic field is generated by tile permanent magnets with uniform polarizations [START_REF] Halbach | Strong rare earth cobalt quadrupoles[END_REF] [START_REF] Ravaud | Discussion about the magnetic field produced by cylindrical halbach structures[END_REF]. However, most of the studies dealing with the magnetic field produced by MRI structures assume that the magnetic field is produced by a ring permanent magnet with a rotating polarization. However, this is not the case in practice because it is very difficult to manufacturate such ring permanent magnets. Consequently, the ring permanent magnet generating the magnetic field is replaced by an assembly of tile permanent magnets uniformly magnetized. As a consequence, the magnetic field magnitude and its homogeneity turn out to be overestimated theoretically because the uniform polarizations of the tile permanent magnets generate a decrease in the homogeneity of the magnetic field. The analytical method we use in this paper is based on the coulombian model of a magnet. This model implies the existence of fictitious magnetic charges located on the faces or in the volume of a permanent magnet [START_REF] Babic | Calculation improvement of 3d linear magnetostatic field based on fictitious magnetic surface charge[END_REF]. Such an analytical method has been widely used by many authors for calculating the magnetic field produced by ring permanent magnets radially and axially magnetized [START_REF] Ravaud | Analytical calculation of the magnetic field created by permanent-magnet rings[END_REF][21], tile permanent magnets radially magnetized [START_REF] Ravaud | The three exact components of the magnetic field created by a radially magnetized tile permanent magnet[END_REF], tile permanent magnets uniformly magnetized [START_REF] Ravaud | Analytical expression of the magnetic field created by tile permanent magnets tangentially magnetized and radial currents in massive disks[END_REF] [24]. First, this paper presents the exact three-dimensional analytical expressions that allow us to study the homogeneity and the magnitude of the magnetic field produced by an assembly of tile permanent magnets uniformly magnetized. Then, we discuss the influence of the number of tile permanent magnets used on the magnetic field properties.

MODELING THE MAGNETIC FIELD INSIDE THE MRI STRUCTURE

We propose in this section to determine the expressions of the magnetic field components in cartesian coordinates for studying the homogeneity of the field inside the MRI structure. For this purpose, let us first consider the configuration shown in Fig. 1. The eight tile permanent magnets are uniformly magnetized. The parameters used in this paper are shown in Fig 2 in which we have represented only one tile permanent magnet uniformly magnetized. We have to consider only one tile permanent magnet whose polarization is both uniform and arbitrary. Indeed, the total magnetic field produced by an assembly As the tile permanent magnets are uniformly magnetized, each tile permanent magnet can be represented by fictitious magnetic pole surface densities located on its faces in the coulombian approach. By denoting σ * j,i , the fictitious magnetic pole surface density located on the i th face of the magnet j whose angular width equals θ j+1 -θ j , the fictitious manetic pole surface densities are given by the scalar product J. n j where j = 1..4 (as shown in Fig 2). After mathematical manipulations, we find σ * j,1 , σ * j,2 , σ * j,3 and σ * j,4 in the following reduced forms:

σ * j,1 = -J cos θ -a - θ j + θ j+1 2 σ * j,2 = -J sin θ j+1 -θ j 2 -a σ * j,3 = J cos θ -a - θ j + θ j+1 2 σ * j,4 = J sin θ j -θ j+1 2 -a (1)
It is useful to define the three-dimensional Green's function as follows:

G 0 ( r, r ) = 1 4π r -r ( 2
)
where r is the observation point and r a point located on the fictitious charge distribution. For the rest of this paper, we use the four following definitions of the three-dimensional Green's function:

G 0 ( r, r ) (1) = (4π) -1 r 2 1 + r2 -2r 1 r cos(θ -θ) + (z -z) 2 -1 2 G 0 ( r, r ) (2) = (4π) -1 r 2 + r2 -2rr cos(θ -θ j+1 ) + (z -z) 2 -1 2 G 0 ( r, r ) (3) = (4π) -1 r 2 2 + r2 -2r 2 r cos(θ -θ) + (z -z) 2 -1 2 G 0 ( r, r ) (4) = (4π) -1 r 2 + r2 -2rr cos(θ -θ j ) + (z -z) 2 -1 2 (3)
It is emphasized here that the use of the three-dimensional Green's function is absolutely equivalent to the coulombian model of a magnet because the structure considered in this paper is ironless. For a tile permanent magnet having an angular width that equals θ j+1 -θ j , the x-component of the magnetc field is given by:

H x (x, y, z) = cos(θ)∂ r 4 k=1 S k G 0 ( r, r ) (k) σ * j,k ( r )d 2 r k -sin(θ)r -1 ∂ θ 4 k=1 S k G 0 ( r, r ) (k) σ * j,k ( r )d 2 r k (4)
The y-component of the magnetc field is given by:

H y (x, y, z) = sin(θ)∂ r 4 k=1 S k G 0 ( r, r ) (k) σ * j,k ( r )d 2 r k + cos(θ)r -1 ∂ θ 4 k=1 S k G 0 ( r, r ) (k) σ * j,k ( r )d 2 r k (5)
The z-component of the magnetic field is given by:

H z (x, y, z) = -∂ z 4 k=1 S k G 0 ( r, r ) (k) σ * j,k ( r )d 2 r k (6)
The previous relations can be developed as follows:

H x (x, y, z) = cos(θ)∂ r θ j+1 θ j z 2 z 1 G 0 ( r, r ) (1) σ * j,1 r 1 d θdz + cos(θ)∂ r r 2 r 1 z 2 z 1 G 0 ( r, r ) (2) σ * j,2 drdz + cos(θ)∂ r θ j+1 θ j z 2 z 1 G 0 ( r, r ) (3) σ * j,3 r 2 d θdz + cos(θ)∂ r r 2 r 1 z 2 z 1 G 0 ( r, r ) (4) σ * j,4 drdz -sin(θ)r -1 ∂ θ θ j+1 θ j z 2 z 1 G 0 ( r, r ) (1) σ * j,1 r 1 d θdz -sin(θ)r -1 ∂ θ r 2 r 1 z 2 z 1 G 0 ( r, r ) (2) σ * j,2 drdz -sin(θ)r -1 ∂ θ θ j+1 θ j z 2 z 1 G 0 ( r, r ) (3) σ * j,3 r 2 d θdz -sin(θ)r -1 ∂ θ r 2 r 1 z 2 z 1 G 0 ( r, r ) (4) σ * j,4 drdz (7) H y (x, y, z) = sin(θ)∂ r θ j+1 θ j z 2 z 1 G 0 ( r, r ) (1) σ * j,1 r 1 d θdz + sin(θ)∂ r r 2 r 1 z 2 z 1 G 0 ( r, r ) (2) σ * j,2 drdz + sin(θ)∂ r θ j+1 θ j z 2 z 1 G 0 ( r, r ) (3) σ * j,3 r 2 d θdz + sin(θ)∂ r r 2 r 1 z 2 z 1 G 0 ( r, r ) (4) σ * j,4 drdz + cos(θ)r -1 ∂ θ θ j+1 θ j z 2 z 1 G 0 ( r, r ) (1) σ * j,1 r 1 d θdz + cos(θ)r -1 ∂ θ r 2 r 1 z 2 z 1 G 0 ( r, r ) (2) σ * j,2 drdz + cos(θ)r -1 ∂ θ θ j+1 θ j z 2 z 1 G 0 ( r, r ) (3) σ * j,3 r 2 d θdz + cos(θ)r -1 ∂ θ r 2 r 1 z 2 z 1 G 0 ( r, r ) (4) σ * j,4 drdz (8) H z (x, y, z) = -∂ z θ j+1 θ j z 2 z 1 G 0 ( r, r ) (1) σ * j,1 r 1 d θdz -∂ z r 2 r 1 z 2 z 1 G 0 ( r, r ) (2) σ * j,2 drdz -∂ z θ j+1 θ j z 2 z 1 G 0 ( r, r ) (3) σ * j,3 r 2 d θdz -∂ z r 2 r 1 z 2 z 1 G 0 ( r, r ) (4) σ * j,4 drdz (9)
After mathematical manipulations, H x (x, y, z), H y (x, y, z) and H z (x, y, z) can be written as follows:

H x (x, y, z) = J cos(θ) 4πµ 0 sin θ j -θ j+1 2 + α f (θ j ) + sin θ j -θ j+1 2 -α f (θ j+1 ) + J cos(θ) 4πµ 0 2 i=1 2 k=1 (-1) (i+k) r i (-z + z k )N[ θ] - J sin(θ) 4πµ 0 sin θ j -θ j+1 2 -α h(θ j ) + sin θ j -θ j+1 2 + α h(θ j+1 ) - J sin(θ) 4πµ 0 θ j+1 θ j 2 i=1 2 k=1 r 2 i (-z + z k ) cos α + θ j -θ j+1 2 -θ sin(θ -θ) (r 2 + r 2 i -2rr i cos(θ -θ)) ( θ) d θ ( 10
)
H y (x, y, z) = J sin(θ) 4πµ 0 sin θ j -θ j+1 2 + α f (θ j ) + sin θ j -θ j+1 2 -α f (θ j+1 ) + J sin(θ) 4πµ 0 2 i=1 2 k=1 (-1) (i+k) r i (-z + z k )N[ θ] + J cos(θ) 4πµ 0 sin θ j -θ j+1 2 -α h(θ j ) + sin θ j -θ j+1 2 + α h(θ j+1 ) + J cos(θ) 4πµ 0 θ j+1 θ j 2 i=1 2 k=1 r 2 i (-z + z k ) cos α + θ j -θ j+1 2 -θ sin(θ -θ) (r 2 + r 2 i -2rr i cos(θ -θ)) ( θ) d θ (11) H z (x, y, z) = J 4πµ 0 sin θ 1 -θ 2 2 -α g(θ 2 ) + sin θ 1 -θ 2 2 + α g(θ 1 ) + J 4πµ 0 θ 2 θ 1 2 i=1 2 k=1 r i cos(α + θ 1 -θ 2 2 -θ) (θ j ) d θ (12)
with .

g(θ j ) = 2 i=1 2 k=1 log [r i -r cos(θ -θ j ) + (θ j )] f (θ j ) = 2 i=1 2 k=1 (-1) (i+k) (-cos(θ -θ j ) log[X j ] + sin(θ -θ j ) arctan[Y j ]) h(θ j ) = 2 i=1 2 k=1 sin(θ -θ j ) log [X j ] -cos(θ -θ j ) arctan [Y j ] (13)
N[ θ] = θ j+1 θ j r -r i cos(θ -θ) cos(α + θ j -θ j+1 2 -θ) r 2 + r 2 i -2rr i cos(θ -θ) ( θ) d θ (14)
X j = z -z k + (θ j ) Y j = (z -z k ) (r i -r cos(θ -θ j )) r sin(θ -θ j ) (θ j ) (β) = r 2 + r 2 i + (z -z k ) 2 -2rr i cos(θ -β) (15)
We represent in 4 show that there is a undulation of the magnetic field that comes from the uniform polarizations of the tile permanent magnets.

The end effects on the magnetic field value

The homogeneity of the magnetic field depends also on the length of the Magnetic Resonance Imaging structure. For example, the finite dimensions along the z-direction (axial direction) generates a curvature in the field, as shown in Fig 5 . Indeed, the magnetic lines are less and less packed together near the edges (z = 0 m and z = 1.6 m). Consequently, this implies a decrease in the y-field magnitude when the observation point is determined far from the middle (z = 0.8 m) of the MRI structure.

Influence of the outer radius on the magnetic field magnitude

Another important parameter that can be used for enhancing the magnetic field value is the outer radius. As shown in Fig 6, the more important the outer radius is, the greater the magnetic field value is. Consequently, despite the expensive cost of permanent magnets, this implies that for a given inner radius r in , the value of the outer radius r out must be the greatest.

INFLUENCE OF THE NUMBER OF TILE PERMANENT MAGNETS ON THE FIELD PROPERTIES

Influence of the number of tile permanent magnets on the field value

This section presents the influence of the number of tile permanent magnets on the field value. Strictly speaking, we take the mean value of the magnetic field on a radial length that equals r = r 1 2 . We take the following dimensions: r 1 = 0.3 m, r 2 = 0.6 m, z 1 = 0 m, z 2 = 1.6 m, z = 0.8 m, J = 1.5 T. Furthermore, the angular width ∆θ of each tile permanent magnet is the same and is given by the number of tile permanent magnets used.

∆θ = θ j+1 -θ j = N 2π ( 16
)
where N is the number of tile permanent magnets used in the structure. We represent in Fig 7 the mean value of the magnetic field < B >= B y = µ 0 H y versus the number of tile permanent magnets used. Figure 7 clearly shows that the greater the number of tile permanent magnets is, the higer the magnetic field is. However, according to this criterion, the evolution of the increase in magnetic field value is not constant but decreases with the number of tile permanent magnets used.

Influence of the number of tile permanent magnets on the field homogeneity

This section presents the influence of the number of tile permanent magnets on the field homogeneity. For this purpose, we define the concept of wave rate W R as follows:

W R = 1 2 B max -B min B moy ( 17
)
where B max and B min are the maximum minimum values of the magnetic field, B moy is the mean value of the magnetic field. According to this criterion, the weaker the wave rate is, the greater the magnetic field homogeneity is. Let us study the influence of the number N of tile permanent magnets on the wave rate. For this purpose, we represent the wave rate versus N in Fig 8 . Figure 8 shows that the wave rate is constant when we use at least 16 tile permanent magnets. Consequently, it is not necessary to use further tile permanent magnets for enhancing the magnetic field homogeneity in the MRI structure with the dimensions taken in this paper.

Two-dimensional representation of the magnetic field magnitude for two axial observation points

We illustrate now the use of our three-dimensional analytical expressions with the 2D representation of the magnetic field for eight radial observation points (r = 0.029 m, r = 0.027 m, r = 0.025 m, Figure 9 clearly shows that the magnetic field undulation is the greatest near the structure, that is, when r = 0.029 m. This undulation wave decreases when the radial observation point goes away from r 1 . Such illustrations give elements of information about the homogeneity of the magnetic field in the MRI structure. Consequently, it gives indications about the region in space in which the magnetic field undulation is weak.

CONCLUSION

This paper has presented exact three-dimensional analytical expressions for studying the magnetic field produced by a MRI structure. These expressions are useful because all the effects due to the uniform polarizations of the tile permanent magnets are taken into account. We have demonstrated that the number of tile permanent magnets used has a great importance on the field value and its homogeneity. For example, with the dimensions taken in this paper, we have seen that the use of at least 20 tile permanent magnets allow one to obtain a high magnetic field. On the other hand, we have shown that the wave rate is not improved when we use at least 16 tile permanent magnets. Fur- thermore, we have discussed the interest of using permanent magnets with great outer radius. Indeed, this allows us to improve significantly the magnetic field magnitude inside the MRI structure.

Figure 1 .Figure 2 .

 12 Figure 1. Illustration of a MRI structure with 8 tile permanent magnets uniformly magnetized

Figure 3 .

 3 Figure 3. Representation of the y-component of the magnetic field versus the angle θ for 16 tile permanent magnets: we take the following r 1 = 0.3 m, r 2 = 0.6 m, z 1 = 0 m, z 2 = 1.6 m, z = 0.8 m, J = 1.5 T, r = r 1 2 .

  Fig 3 the y-component of the magnetic field versus the angle θ and in Fig 4 the x-component of the magnetic field versus θ. Figures 3 and

Figure 4 .Figure 5 .

 45 Figure 4. Representation of the x-component of the magnetic field versus the angle θ for 16 tile permanent magnets: we take the following dimensions: r 1 = 0.3 m, r 2 = 0.6 m, z 1 = 0 m, z 2 = 1.6 m, z = 0.8 m, J = 1.5 T, r = r 1 2 .

Figure 6 .

 6 Figure 6. Representation of By versus r out [m] with 16 tile permanent magnets: we take the following dimensions: r 1 = 0.3 m, z 1 = 0 m, z 2 = 1.6 m, z = 0.8 m, J = 1.5 T, r = 0 m.

Figure 7 .

 7 Figure 7. Representation of the magnetic field By versus the number of tile permanent magnets used: we take the following dimensions: r 1 = 0.3 m, r 2 = 0.6 m, z 1 = 0 m, z 2 = 1.6 m, z = 0.8 m, J = 1.5 T, r = r 1 2 .

Figure 8 .

 8 Figure 8. Representation of the wave rate of the magnetic field versus the number of tile permanent magnets used: we take the following dimensions: r 1 = 0.3 m, r 2 = 0.6 m, z 1 = 0 m, z 2 = 1.6 m, z = 0.8 m, J = 1.5 T, r = r 1 2 .

Figure 9 .

 9 Figure 9. Representation of the magnetic field for 8 radial observation points (r = 0.029 m, r = 0.027 m, r = 0.025 m, r = 0.023 m, r = 0.021 m, r = 0.019 m, r = 0.017 m, r = 0.015 m)at a given axial observation point (z = 0.08 m): we take the following dimensions: r 1 = 0.3 m, r 2 = 0.6 m, z 1 = 0 m, z 2 = 1.6 m, J = 1.5 T.

"Magnetic field produced by a tile permanent magnet whose polarization is both uniform and tangential," Progress in Electromagnetics Research B, vol. 13, pp. 1-20, 2009.