Zhenyu Yin

Pierre Yves Hicher

Micromechanics-based model for cement-treated clays

Keywords: microstructure, cementation, constitutive model, clay, anisotropy, debonding

Cementation is produced by mixing a certain amount of cement with the saturated clay. The purpose of this paper is to model the cementation effect on the mechanical behavior of cement-treated clay. A micromechanical stress-strain model is developed considering explicitly the cementation at inter-cluster contacts. The inter-cluster bonding and debonding during mechanical loading are introduced in two ways: an additional cohesion in the shear sliding and a higher yield stress in normal compression. The model is used to simulate isotropic compression and undrained triaxial tests under various confining stresses on cement-treated Singapore clay with various cement contents. The applicability of the present model is evaluated through comparisons between numerical and experimental results. The evolution of local stresses and local strains in inter-cluster planes is discussed in order to explain the induced anisotropy due to debonding at contact level under the applied loads.

Stabilization of soft ground by deep cement mixing and jet grouting methods, which both introduce cement into the ground, has been widely used in geotechnical projects. The proper understanding of the cementation influence on the mechanical behavior of cement-treated clay is of great importance. Experimental results show an increase of shear strength and shear modulus with cement content, accompanied by a large post-peak reduction in the strength of the treated soil. [START_REF] Mitchell | Proc Int. Conf. Soil Mech. Found. Eng[END_REF][START_REF] Tatsuoka | [END_REF][START_REF] Balasubramaniam | Proc. of 11th Asian Regional Conf. Soil Mech. Found. Eng. Balkema[END_REF][START_REF] Horpibulsuk | [END_REF][5] Cementation also makes the yield stress and swelling index higher than those of untreated soil. Bond degradation is often observed if the compression load exceeds the yield stress, as it has also been observed in natural structured soils. [6][7][8][9][10][11][12] In this study, the recently developed micromechanical model of remolded clay is first briefly introduced. The model is then extended by considering the cementation effect on the mechanical properties of cementtreated clays. The model is used to simulate isotropic compression and undrained triaxial tests under various confining stresses on cement-treated clays with different cement contents. The evolution of local stresses and strains in inter-cluster planes due to externally applied loads is also discussed in order to describe the mechanism of induced anisotropy due to debonding at the micro contact level.

Yin et al. [13][14][15][16][17][18] assumed that clay material can be considered as an assembly of clusters. The deformation of a representative volume of clay is generated by mobilizing and compressing of all clusters. Thus, the stressstrain relationship can be derived as an average of the deformation of all local contact planes. The model includes the inter-cluster behavior, the influence of density state, and the overall stress-strain relationship. 16 For the α-th contact plane, the local forces f α j and the local movements δ α i can be denoted as follows: 1 for local coordinate system). In order to obtain a more direct comparison between the local behavior and the overall stress-strain behavior, we define a local normal stress σ α = f α n Nl/(3V ) and a local shear stress τ α = f α s Nl/(3V ), where l is branch length and N/V is the total number of contact per unit volume. The corresponding local normal strain is defined as ε α = δ α n /l and a local shear strain is defined as γ α = δ α s /l. Isotropic compression and oedometer tests by Horpibulsuk et al. [START_REF] Horpibulsuk | [END_REF] and Kamruzzaman et al. 5 on two cement-treated clays show the significant influence of the cement content on the swelling index (see Fig. 2(a)) and on the compression yield stress (see Fig. 2(b)). This influence can be expressed by Eq. ( 1) for the swelling index κ and by Eq. ( 2) for the bonding pressure p b defined as the difference of the compression yield stress between cement-treated clay p yc and untreated clay p yu : p b = p yc -p yu , as follows

f α j = {f α n , f α s , f α t } and δ α i = {δ α n , δ α s , δ α t } (see Fig.
κ = κ 0 exp (-β κ c) , ( 1
)
p b = β p1 [exp (β p2 c) -1] , ( 2
)
where κ 0 is the swelling index for untreated clay, c is the cement content, β κ , β p1 and β p2 are material constants (see Figs. By comparing the post-yield stress-strain curves of cement-treated samples to those of untreated samples, 4,5 a mechanical bond degradation process can be observed. This bond degradation is an irreversible phenomenon that, experimentally, appears to be controlled by plastic strain accumulation. 13,19 The bonding and debonding can be linked directly to the inter-cluster bonds due to the cementation formed when adding cement into clay slurry.

The friction angle at failure was measured based on undrained triaxial tests by Horpibulsuk et al. [START_REF] Horpibulsuk | [END_REF] and

p b =β p1 [exp(β p2 c)-1] c b =β c1 [exp(β p2 c)-1] κ 0 =0.092, β κ =4.2 φ μ0 =40Ο, φ μmax =40Ο, β f =0 φ μ0 =23Ο, φ μmax =43Ο, βf=3.4 φ μ =φ μ0 +(φ μmax-φ μ0 )[1-exp(β f c)] β p1 =10, β p2 =30 β c1 =110, β c2 =10 β c1 =9 000, β c2 =0.1 β p1 =850, β p2 =2 κ 0 =0.086, β κ =11 φ μ /(Ο)
Fig. 2. Effect of cement content on properties of cement treated clays.

Kamruzzaman et al. 5 on both cement-treated clays by extending the curves to an axial strain up to 15 % according to the previous trend. As shown in Fig. 2(c), for Ariake clay, having an initially high friction angle, the cement content influences slightly the friction angle at failure (we obtained φ μ = 40 • measured from undrained stress paths, which is slightly different from φ μ = 38 • provided by Horpibulsuk et al. [START_REF] Horpibulsuk | [END_REF] ). However, for Singapore clay with an initially small friction angle, the cement content influences significantly the friction angle at failure. This influence can be expressed by Eq. ( 3) as follows

φ μ = φ μ0 + (φ μ max -φ μ0 ) [1 -exp (-β f c)] , (3)
where φ μ0 is the friction angle at failure for the untreated clay, φ μ max is the maximum friction angle at failure for the cement-treated clay, and β f is a material constant.

The apparent cohesion due to cement bonds is measured from the position of the failure line in the p -q plane. Figure 2(d) shows the significant influence of the cement content on the apparent cohesion c b for both cement-treated clays. This influence can be expressed by Eq. ( 4), having a form similar to Eq. ( 2) for the bonding pressure, as follows

c b = β c1 [exp (β c2 c) -1] , (4)
where β c1 and β c2 are material constants. The stress-strain curves of the two cement-treated clays show also the influence of the cement content on the initial shear modulus, which is similar to its influence on the swelling index (or bulk modulus). Since the shear modulus can not be accurately measured based on digitized data, we will not discuss it in this paper. Upon increase of the applied shear stresses, the intercluster bonds are progressively damaged. As a result, a decrease of the deviatoric stress takes place after the peak for cement-treated samples.

The critical state lines (CSL) are obtained from the undrained triaxial tests for untreated clays, as shown in Figs. 3(a) and 3(b). The void ratio corresponding to this state is e c . The CSL can be written as follows for clay

e c = e cr -λ lg p p ref , ( 5
)
where p is the mean effective stress, λ is the slope of the CSL in the e-lgp plane, e cr is the reference void ratio at critical state corresponding to p ref = 100 kPa.

For cement-treated clays, if the failure state in the p -q plane (corresponding to 15 % of axial strain) is considered as the critical state, the CSL in the e-lgp plane can be obtained for all cement-treated samples. For cement-treated Ariake clay (see Fig. 3(a)), the critical state line can be assumed to be parallel to that of the untreated clay. If this assumption is also adopted for cement-treated Singapore clay, the critical state lines can be obtained as shown in Fig. 3(b). Therefore, the influence of cement content on the location of the critical state line can be obtained by making the reference void ratio function of the cement content, as shown in Fig. 3(c), expressed as follows e cr = e cr0 + (e cr max -e cr0 ) [1exp (-β e c)] , (6) where e cr0 is the reference critical state void ratio for the untreated clay, e cr max is the maximum reference critical state void ratio for cement-treated clay, β e is a material constant. Horpibulsuk et al. [START_REF] Horpibulsuk | [END_REF] have shown that the clay fabric consists of clay cluster (aggregation of clay particles) and that the role of the induced cementation is to weld clay clusters. Kamruzzaman et al. 5 have made SEM photos untreated and cement treated by Singapore clay with various cement contents under various consolidation pressures, which demonstrate that the size of the clay clusters can be considered the same for all cases. Therefore, the cement-treated clay can be seen as an assembly of clusters (aggregates of clay-cement particles). The deformation of an assembly can be obtained by integrating the movement of the inter-cluster contacts in all orientations. Thus, the effect of the intercluster bonding can be explicitly taken into account. As the orientation dependent properties are explicitly introduced, the induced anisotropy due to inter-cluster debonding can be modeled in a direct way. Therefore, the micromechanics-based model developed by Yin et al. 13,14 for natural clays can be extended for cementtreated clays. In addition to the micromechanical model of remolded clay, three modifications were made to account for the effect of cement content on the intercluster bonding and debonding, listed below.

In order to take into account the effect of the intercluster bonding, the plastic law of the shear sliding at inter-cluster contacts proposed by Yin et al. 14,15 is adopted. The yield function for shear sliding is assumed to be of Mohr-Coulomb type, written in a contact plane (e.g. σ, τ s , τ t , see coordinate system in Fig. 1) as follows

F 1 (σ, τ, c b , H 1 ) = τ σ + c b /tan φ μ -H 1 (γ p ) , (7)
where c b is the inter-cluster cohesion due to the intercluster bonding, H 1 (γ p ) is a hardening/softening function. The values of c b and φ μ depend upon the cement content (see Eqs. ( 4) and ( 3)). The influence of the cement content on the shear behavior is thus introduced. For untreated clay, c b = 0, and the equation can be reduced to that by Yin et al. [15][16][17][18] The local dilatancy equation is derived from that proposed by Yin et al. [16][17][18][19] as follows

d ε p d γ p = D (tan φ μ -tan φ m ) tan φ m tan φ μ • e e c -1 , ( 8
)
where tan φ m = τ /(σ +c b / tan φ μ ), D is a material constant which controls the dilatancy rate during shearing. The void ratio at critical state e c is calculated using Eq. ( 5) where the influence of the cement content is introduced by e cr (see Eq. ( 6)).

The degradation of the inter-cluster bonds can be modeled as a damage of the bonded contacts. Therefore, the damage law proposed by Yin et al. 13,14 is adopted in the expression of the inter-cluster cohesion, as follows

c b = c b0 exp (-ξ d γ p ) , ( 9
)
where c b0 is the initial inter-cluster cohesion, ξ d is the factor of damage representing the influence of the tangential displacement in the damage law. Therefore, during loading each contact sliding produces a progressive damage of the cohesion.

In order to account for the inter-cluster bonding effect on the compressive behavior of cement-treated clay, we adopt the plastic law for the normal compression of inter-cluster contacts proposed by Yin et al. 13,14

F 2 (σ, σ p ) = σ -σ p , ( 10
)
where σ p = σ pi (1+χ) is the compression yield stress, σ pi is the intrinsic compression yield stress corresponding to the yield stress measured from an isotropic compression test on a reconstituted sample of untreated clay, χ is the bonding ratio defined by χ = σ p /σ pi -1. The value of σ p depends on the cement content (σ p =σ pi +p b , see Eq. ( 2) for p b ). Thus, the influence of the cement content on the compression behavior is introduced. The hardening function controlling the evolution of σ pi is defined as follows

σ pi = σ pi0 exp ε p c p , ( 11
)
where c p is the compression index determined from the compression curve plotted in the ε p -lgσ plane. The damage law proposed by Yin et al. 13,14 is also adopted in the expression of the bonding ratio at intercluster level

χ = χ 0 exp (-ξ n ε p ) , ( 12
)
where χ 0 is the initial bonding ratio depending of the cement content, ξ n is the factor of damage representing the influence of the normal displacement in the damage law. Therefore, during loading each contact produces a progressive damage of the bonding ratio. For untreated clay χ 0 = 0 (or σ p = σ pi ), the equation can be reduced to that by Yin et al. [15][16][17][18] Based on experimental investigation, the dependency of the CSL location on the cement-content is implemented into the model, i.e. the Eq. ( 12) is replaced by the Eqs. ( 5) and (6).

Combined all above equations with others by Yin et al., 13,14 the model can be used to reproduce mechanical behavior of cement-treated clays.

The determination of parameters for cement-treated Ariake clay was based on the isotropic compression and undrained triaxial tests under a consolidated stress of 200 kPa on the untreated Ariake clay. All determined parameters, summarized in Table 1, were used to predict all the other tests. Figure 4 shows the comparison between experimental results and numerical simulations of isotropic compression tests on untreated and cement-treated Ariake clay with different cement content. A satisfactory agreement is achieved by using the corresponding set of parameters in Table 1.

ε a /(%) ε a /(%) ε a /(%) ε a /(%) (a
Figure 5 shows the comparisons between experiments and simulations for undrained triaxial tests on samples of Ariake clay with different cement contents consolidated up to 200 kPa. The influence of the cement content on the behavior of cement-treated Ariake clay was well reproduced by the proposed model: (1) The increase of stiffness and shear strength with cement content, accompanied by a large post-peak reduction in the strength of the treated soil; (2) Samples with low cement content (c = 0 %, 6 %, 9 %, 12 %) exhibit a contractive behavior and samples with high cement content (c = 18 %) a dilative one.

The distributions of the local normal stiffness and bond in contact planes of various orientations are plotted in Fig. 6 (Fig. 6(1a) for untreated clay, Fig. 6(2a) for cement content of 9 % and Fig. 6(3a) for cement content of 18 %):

(a) The normal stiffness is a function of the normal stress and of the amount of cement. The normal stiffness distribution at the end of the isotropic consolidation has a circular shape (see the bold line in Fig. 6(a)) which corresponds to an isotropic distribution of the normal stress for all contact planes. During the undrained shearing, different cement contents give different evolutions of the stiffness distribution: For 0 % cement content, the distribution becomes anisotropic with an elliptical shape having the long axis in the vertical direction; For 9 % cement content, the distribution becomes also anisotropic with higher stiffness in the vertical direction and lower stiffness in the horizontal direction, and then shrinks during shearing due to bond damage; The same evolution is found for 18 % cement content with higher stiffness values due to stronger bonds.

(b) During the undrained shearing, the shear stress distribution expands while maintaining a similar shape for all cases (Fig. 6(b)). This distribution shows a concentration of the maximum local shear stresses between the 45 • and 55 • orientations. Step 3

Step 2

Step 1

Step 3

Step 2

Step 1

Step 3 Step 2 cb=0 (2a) Normal stress σ n

Step 3

Step 2

Step 1 (2b) Shear stress τ

Step 3

Step 2

Step 1 (2c) Shear strain γ

Step 3 (2d) Shear bonding cohesion cb

Step 3

Step 2

Step 1 (3a) Normal stress σ n

Step 3

Step 2

Step 1 (3b) Shear stress τ

Step 3

Step 2 Step 1 (3c) Shear strain γ

Step 3 (3d) Shear bonding cohesion cb (c) During the undrained shearing, the shear strain distribution for all cases (Fig. 6(c)) shows that very large strains occurred at the end of the loading in contact planes around the 70 • orientation.

(d) The distribution of shear bonding cohesion changes slightly from step 1 to 2 due to small mobilized shear plastic strains for both 9% and 18% cement contents. Then, from step 2 to 3, the bonding cohesion reduces significantly in the 55 • -70 • contact planes, but does not change very much in the other planes. This agrees with the shear strain distribution in Figs. 6(2d) and 6(3d), because the damage of the shear bonding cohesion is controlled by the amount of plastic shear strain.

In the present model, stresses and bonds in each plane are considered as internal state variables and their different evolution in each individual plane introduces, in a natural manner, the stress-induced anisotropy. Since many soil properties are stress-dependent (e.g. the plastic hardening parameter k p and the shear modulus k r are related to k n which is stress-dependent), the induced anisotropic behavior during undrained shearing is linked to the stress distribution in each plane.

Based on experimental results, a micromechanical stress-strain model was developed considering explicitly the cementation at inter-cluster contacts. The inter-cluster bonding is introduced by considering two aspects: an additional cohesion in shear sliding and a higher yield stress in normal compression. Then, isotropic compression and undrained triaxial tests under various confining stresses on cement-treated Ariake clay with different cement contents were simulated to evaluate the capability of the present model to take into account the changes in mechanical properties linked to the cement content. Finally, the predicted behavior in contact planes was examined in the case of undrained triaxial tests with different cement content, which clearly indicates the development of anisotropy induced by the externally applied loads.

  2(a) and 2(b)).

Fig. 1 . 3 κ=κ 0

 130 Fig. 1. Local coordinate at inter-particle contact.

9 e 4 Fig. 3 .

 943 Fig. 3. Effect of cement content on the position of critical state line.

Fig. 4 .

 4 Fig. 4. Comparisons between experiments and simulations for isotropic compression tests on cement treated Ariake clay with cement content of 0%, 6%, 9%, 12%, 18%.

Fig. 5 .

 5 Fig. 5. Effect of cement content on undrained shearing behaviour of cement treated Ariake clay for tests under confining pressure of 200 kPa.

( 1a )

 1a Normal stress σ n (1b) Shear stress τ (1c) Shear strain γ (1d) Shear bonding cohesion cb

Fig. 6 .

 6 Fig. 6. Local stresses, strains and bonds distributions in rose diagram, (1) is for untreated clay, (2) is for cement content of 9 % and (3) is for cement content of 18 %.

  This work was supported by Research Fund for the Doctoral Program of Higher Education of China (20110073120012) and Shanghai Pujiang Talent plan (11PJ1405700).

Table 1 .

 1 Local coordinate at inter-particle contact.

	Parameters		e0		κ0		λ	σ pi0		φμ	kp		D	e cr0	ξn	ξ d
	Ariake clay		3.44		0.086	0.446	10		40 •	0.15		6	1.911	10	3
		1 600						1 600						400
															9%
	p/kPa	0 400 800 1 200	9% 6% 0%	12%		18%	p/kPa	0 400 800 1 200	9%	18%	0% 6%	-100 0 300 100 200 Du/kPa	6% 0%	18%
		0	200	400	600		0	5	10	15	20		0	5	10	15	20
				p'/kPa
		1 600						1 600						400	12%
	p/kPa	400 800 1 200	9% 6%	12%		18%	p/kPa	400 800 1 200			18%	6% 0% 12% 9%	Du/kPa	0 300 100 200	9% 6% 0%	18%
		0	0%					0						-100
		0	200	400	600		0	5	10	15	20		0	5	10	15	20
				p'/kPa


Qi-Yin Zhu

Yin-Fu Jin

Zhen-Yu Yin

Pierre-Yves Hicher

Influence of natural deposition plane orientation on oedometric consolidation behavior of three typical clays from southeast coast of China *

Keywords: Clay, Compressibility, Consolidation, Creep, Permeability, Natural deposition

The parameters obtained from oedometric consolidation tests are commonly used in the development of constitutive modeling and for engineering practice. This paper focuses on the influence of the natural deposition plane orientation on oedometric consolidation behavior of three natural clays from the southeast coast of China. Oedometer tests were conducted on intact specimens prepared by sampling at a series of angles relative to the natural deposition plane. For each specimen, yield stress, compressibility indexes, secondary compression, and permeability coefficients were determined. The influence of the sampling angle on these properties was investigated, revealing that yield stress, compression index, swelling index, creep index, ratio of secondary compression coefficient to compression index (C e /C c ) and permeability coefficient were all dependent to some extent on the sampling angle. These findings indicate the role of the anisotropy due to the natural deposition on the oedometric consolidation behavior.

Introduction

Natural clayey deposits are widely distributed in the coastal regions of southeast China. In general, these deposits show diverse ground characteristics due to variations in depositional environment. Moreover, their post-depositional geotechnical characteristics, such as stress history, leaching processes, and variations in excess pore water pressure, are very complex. These natural deposits are usually extremely soft and present considerable problems in civil engineering construction. Thus, the design and construction of projects built on these natural deposits require complete geotechnical knowledge of the natural clays therein to prevent excessive settlementdeformation and associated structural damage [START_REF] Zhang | Achievements and problems of geotechnical engineering investigation in China[END_REF]. Presently, the parameters used for design and construction are usually determined from laboratory tests on clay samples that have been cut vertically. However, the mechanical behavior of natural clays cut at different angles from the deposition plane can strongly vary [START_REF] Lade | Effects of stress rotation and changes of b-values on cross-anisotropic behaviour of natural, K 0 -consolidated soft clay[END_REF][START_REF] Shen | Study on Strength Criterion of Intact Soft Clay after Monotonic Principal Stress Rotation[END_REF]Yin and Chang, 2009a;2009b;Yin et al., 2009;[START_REF] Huang | Simulation of yield characteristics and principal stress rotation effects of natural soft clay[END_REF]. Therefore, the present design methods using the parameters measured from vertical clay samples may lack reliability and are likely to generate inaccurate predictions, resulting in engineering problems.

Parameters from oedometric consolidation tests are the most basic parameters commonly used in the development of elastoplastic or viscoplastic models and for engineering practice [START_REF] Chen | Consolidation of sensitive clay with vertical drain[END_REF][START_REF] Yin | Identifying parameters controlling soil delayed behaviour from laboratory and in situ pressuremeter testing[END_REF][START_REF] Li | Environmental effects induced by excavation[END_REF][START_REF] Wang | Undrained shear strength of K 0 consolidated soft soils[END_REF][START_REF] Wang | Construction behavior simulation of a hydraulic tunnel during standpipe lifting[END_REF]2012;[START_REF] Li | One-dimensional viscous behavior of clay and its constitutive modeling[END_REF][START_REF] Karstunen | Modelling time-dependent behaviour of Murro test embankment[END_REF][START_REF] Feng | Some observations on the oedometric consolidation strain rate behaviors of saturated clay[END_REF][START_REF] Yin | An anisotropic elastic-viscoplastic model for soft clays[END_REF][START_REF] Zhang | Achievements and problems of geotechnical engineering investigation in China[END_REF][START_REF] Duan | Theoretical elastoplastic analysis for foundations with geosyntheticencased columns[END_REF][START_REF] Wang | Stress-dilatancy of natural soft clay under undrained creep condition[END_REF][START_REF] Yin | Elastic viscoplastic modeling for natural soft clays considering nonlinear creep[END_REF]. For instance, the compression and swelling indexes, the secondary compression coefficient and the permeability are all oedometric properties that may potentially influence decisions prior to initiating civil engineering projects. For this reason, oedometric properties of clays have been widely investigated (e.g., [START_REF] Miao | Timedependent deformation behavior of Jiangsu marine clay[END_REF][START_REF] Zeng | Change of hydraulic conductivity during compression of undisturbed and remolded clays[END_REF][START_REF] Li | Low secondary compressibility and shear strength of Shanghai Clay[END_REF]. However, these experimental studies have primarily been conducted on clay samples cut vertically, therefore having an angle with the natural deposition plane depending on its in situ orientation, and few have investigated the influence of the sampling angle relative to the natural deposition plane on oedometric parameters of clays.

In this study, we performed oedometer tests carried out on samples of three natural clays from Southeast China coastal regions sampled at various angles relative to the natural deposition plane. The effect of the sampling angle was experimentally quantified, and then used to study the oedometric consolidation behavior of clays, including compression, swelling, creep, and permeability properties. The implications of these findings for construction projects on clay deposits were discussed.

Clay samples and test program

Sampling and physical properties of clays

The clay samples, hereafter designated based on their origins, were taken from the Minhang District of Shanghai, Liuheng Island at Zhoushan, and Puzhou at Wenzhou (Fig. 1). Some basic physical properties of the tested clays are summarized in Table 1, indicating significant silt content. Considering the percentages of silt and sand, we can describe Shanghai clay as silt clay and Zhoushan clay as slightly sandy. The Casagrande (1936)'s plasticity chart was used to classify the three clays; all are above the A-line and can therefore be considered as inorganic clays (Fig. 2). The natural deposition has a significant influence on the anisotropy of clay structures. Fig. 3 shows the photos from scanning electronic microscopy (SEM) of three natural depositions in the vertical direction. These photos show clearly that the particles arrangement of natural clays is anisotropic. Using the software 'ImageJ', we analyzed the orientation of the clay particles on the vertical plane. An angular distribution plotted in a rose diagram gives the percentage of the particle surface area as a function of the particle orientation. The rose diagrams indicate a structural anisotropy with a preferential orientation of the particles towards the horizontal direction. A certain number of particles non-uniformly oriented in the 30°-150° quadrants can also be observed (Fig. 4).

Test program

Fig. 5 illustrates a block of natural soil sampled from the field with the drive sampling techniques described by [START_REF] Clayton | Site Investigation[END_REF]. The sedimentation, or natural deposition plane of the clay is assumed to mostly coincide with the in situ horizontal plane. To characterize the influence of the sampling angle on the consolidation behavior of clays, specimens prepared for oedometer testing were cut from the soil block at several different angles. The specimens of Zhoushan clay were sampled at angle of 0° (i.e., perpendicular to the deposition plane), 30°, 60°, and 90° (i.e., parallel to the deposition plane), while specimens of Shanghai and Wenzhou clays were sampled at =0°, 22°, 45°, 68°, and 90°.

Specimens were subjected to conventional oedometer tests according to Chinese Industry Standard (1999). All laboratory tests were conducted in a standard fixed-ring oedometer with 30 cm 2 brass rings (diameter d=61.8 mm) and 20 mm in thickness. Each specimen was subjected to successive vertical stresses of 12.5, 25, 50, 100, 200, 400, 800, 200, 50, 12.5, 50, 200, 800, and 1600 kPa, corresponding to a loading stage (12.5 to 800 kPa), an unloading stage (800 to 12.5 kPa) and a reloading stage (12.5 to 1600 kPa). Each stress step lasted 24 h.

Experimental results and analysis

Yield stress

Fig. 6 shows the compression curves (e-log v ) of specimens of the three clays with different sampling angles relative to the natural deposition plane. Differences among the four compression curves for Zhoushan clay were observed, with the compression curve of the 0° specimen exhibiting the smallest change in void ratio and the compression curve of the 60° specimen exhibiting the biggest change (Fig. 6a). In contrast, the five compression curves for Shanghai and Wenzhou clay specimens were highly similar, regardless of the sampling angle (Figs. 6b and6c).

Based on these findings, the yield stress ( p0 ) of each specimen was determined using the Casagrande method [START_REF] Casagrande | The Determination of the Preconsolidation Load and Its Practical Significance[END_REF]. Fig. 6d shows the change in yield stresses of Zhoushan, Shanghai and Wenzhou

Compressibility

The compression (C c ) and swelling (C s ) indexes can be measured from the compression curves (Fig. 6). C c is measured from the slope of the compression line between 200 and 800 kPa. Fig. 7a shows the change in C c with the sampling angle for the three clays, along with the correlation coefficients for each linear regression, indicating overall that C c varies very slightly with the sampling angle. Two values of C s were measured: C si related to the slope of the initial compression line between 12.5 and 25 kPa fully induced by the natural clay deposition, and C sr related to the slope of the unloading and recompression lines additionally influenced by the first loading stage. Fig. 7b shows the change in C si with the sampling angle for the three clays, putting into evidence the relatively large variations of these index values for each clay sample. Fig. 7c indicates that C sr decreases slightly with the sampling angle.

Secondary compression

The change in void ratio occurring over one log cycle of time during the secondary compression phase is known as

C e : e , log e C t ( 1
)
where e is the void ratio and t is the length of time since the end of the primary consolidation. C e is often used to assess the long-term deformation of clays, as it can account for the creep rate of soils [START_REF] Yin | An anisotropic elastic-viscoplastic model for soft clays[END_REF]. We therefore determined C e for a variety of sampling angles for the three clays. Fig. 8 shows an example of C e values determined from plots of e vs. logt at each load increment for the 0° specimen of Shanghai clay. The duration of the primary consolidation phase for Shanghai clay was determined to be about 15 min, followed by the secondary consolidation. Based on the evolution of e over time during the secondary consolidation stage, C e was calculated for each load increment using Eq. ( 1). This approach was used to quantify C e for all specimens.

Fig. 9 shows changes in C e for three clays sampled at different angles relative to the deposition plane, based on levels of vertical stress applied. For all three clays, C e increases with the increase in the consolidation stress regardless of the sampling angle and reaches a maximum at a stress about twice the preconsolidation pressure; thereafter, C e decreases irrespective of clay sample or sampling angle. This observation is consistent with previous studies [START_REF] Mesri | Time and stress-compressibility interrelationship[END_REF][START_REF] Suneel | Compressibility characteristics of Korean marine clay[END_REF] on various clays in the vertical direction.

Note that, for each clay, all samples have approximately the same initial void ratio, then the relationship between secondary compression C e and vertical stress can also be considered an expression of the relationship between C e and void ratio upon completion of the test. We analyzed the relationship between C e and the void ratio within the normally consolidated range of the clays, as shown in Fig. 10. C e appears to decrease roughly linearly with the void ratio in the double logarithm plane, which is in agreement with the evolution suggested by [START_REF] Wang | Stress-dilatancy of natural soft clay under undrained creep condition[END_REF]. Meanwhile, the sampling angle relative to the natural deposition plane has a visible influence on the relationship between C e and the void ratio that also depends on the clay. For Zhoushan clay, C e is the lowest for the 30° specimen (Fig. 10a). However, for Shanghai clay, C e generally decreases with an increase in the sampling angle (Fig. 10b), and C e in Wenzhou clay increases with an increase in the sampling angle, apart from the 0° specimen, which has the highest C e for this clay (Fig. 10c).

2 -1.8 -1.6 -1.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0°2 2°4 5°6 8°9 0°( c) -2.3 -2.1 -1.9 -1.7 -1.5 -0.4 -0.3 -0.2 -0.1 0 loge 0°2 2°4 5°6 8°9 0°( b) -2.3 -2.2 -2.1 -2 -1.9 -1.8 -0.5 -0.4 -0.3 -0.2 -0.1 0 0°3 0°6 0°9 0°( a)
For each clay we selected the maximum value of C e and the value of C e at the vertical stress of 1600 kPa to represent the creep properties of the natural deposit from which the samples had been taken. Plotting the relationship between the maximum of C e and the sampling angle for each of the three clays shows that C e varies from 0.024 to 0.033 for the highly plastic Wenzhou clay (Fig. 11a), whereas C e has a much smaller range for the low plasticity Shanghai clay (C e =0.014-0.018) and Zhoushan clay (C e =0.010-0.012). Notably, C e corresponding to a vertical stress of 1600 kPa varies slightly with the deposition plane angle, likely because the soil microstructure was modified towards a similar state during loading up to a high stress level (Fig. 11b). [START_REF] Mesri | Time and stress-compressibility interrelationship[END_REF] summarized previously published C e /C c values for a number of natural soils and found that this ratio ranged overall between 0.025-0.010. [START_REF] Mesri | C /C c concept and K 0 during secondary compression[END_REF] Zhoushan clay with an average value of 0.032 (Fig. 12a), between 0.025 and 0.052 for Shanghai clay with an average value of 0.036 (Fig. 12b), and between 0.018 and 0.047 for Wenzhou Clay with an average value of 0.03 (Fig. 12c). Tables inset within Fig. 12 list the average values of C e /C c for each of the clay sampling angles, and these data were used to plot the change in C e /C c as a function of the sampling angle (Fig. 12d). For the three clays, C e /C c first decreases, then increases as the sampling angle increased, with 0° and 90° samples generating the largest C e /C c .

C e /C c ratio

Permeability

The permeability coefficient k is calculated by

v w v , k C m (2)
where C v is the coefficient of consolidation, w is the unit weight of water and m v is the coefficient of volume change. Since C v and m v varied with the sampling angle relative to the deposition plane, we calculated C v from strain-time curves using the Terzaghi 1D consolidation theory. Fig. 13 shows the e-logk relationships for Zhoushan, Shanghai and Wenzhou clays tested at different sampling angles relative to the natural deposition plane. The results show that the permeability coefficient k decreases linearly with the decrease in the void ratio in a logarithmic plane for all specimens.

The initial permeability k 0 and the permeability index C k are widely used in engineering calculations; hence, their changes in the sampling angle were used to evaluate the effects of the orientation of the natural deposition plane. Eq. (3) proposed by [START_REF] Berry | The consolidation of peat[END_REF], i.e., was used to determine the two permeability parameters k 0 and C k based on k, e, and e 0 . Given the initial void ratio e 0 and the linear regression between e and logk for different sampling angles, the corresponding k 0 could be calculated. Meanwhile, C k is the slope of the corresponding regression line, such as that plotted in Fig. 13a. We thereby determined k 0 and C k for all samples. Fig. 14a shows that the sampling angle in-fluences k 0 in a positive, linear manner, with regression slope values of 0.0089, 0.0074, and 0.0050 for Zhoushan, Shanghai, and Wenzhou clays, respectively. The k 0 value for the vertical samples (0°) is the smallest, and the largest for the horizontal samples (90°) for all three clays. Previous results [START_REF] O'kelly | Compression and consolidation anisotropy of some soft soils[END_REF] showed that the ratio of permeability values for horizontal relative to vertical samples is always larger than 1.0, and the present data are consistent with this finding. On the other hand, plotting C k versus the sampling angle shows no significant influence of the latter on C k (Fig. 14b).

Discussion

We found that the yield stress of Shanghai clay gradually decreases with the sampling angle, with about 17 kPa difference between the maximum and minimum p0 . For Zhoushan and Wenzhou clays, however, p0 does not strongly vary. Furthermore, the highly plastic Wenzhou clay has a large range of C e that varies with the sampling angle, while the low plasticity Zhoushan and Shanghai clays have small ranges of C e values. Similar phenomena can be found for C si and C k . Consequently, we can conclude that the anisotropic behavior is more obvious for clays of highly plasticity.

For isotropic compression, experimental evidence suggests that the initial anisotropy is erased by isotropic loading up to a pressure that is two or three times larger than the preconsolidation pressure [START_REF] Leoni | Anisotropic creep model for soft soils[END_REF]. If adopting this in an oedometer test, by assuming K 0 =0.5 (friction angle =30°), the initial anisotropy will be removed after a vertical stress larger than 1.3 or 2 times preconsolidation pressure and the current conditions can be regarded as yielding. Fig. 11 shows the relationships between the maximum value of C e (before yielding) and the value of C e at the vertical stress of 1600 kPa (after yielding) and the sampling angle. The significant influence of initial anisotropy of clays before yielding on the anisotropic behavior of C e can be observed. As the compression index C c and swelling index C sr was measured after yielding, and C si was measured before yielding of clay, C c and C sr vary slightly with the sampling angle but not C si as shown in Fig. 7, also demonstrating the rationality of the above observation. In addition, the evolution of initial permeability k 0 with the sampling angle for three clays also gives positive supports (Fig. 14).

Conclusions

The effect of the orientation of the natural deposition plane on oedometric consolidation behavior of three natural clays from coastal regions of Southeast China was investigated by conducting a series of oedometer tests on specimens prepared by sampling a single clay block at a series of angles relative to the natural deposition plane. Conclusions can be drawn as follows.

1. The anisotropic behavior is more obvious for clays of highly plasticity. The highly plastic Wenzhou clay has a large range of C e that varies with the sampling angle, while the low plasticity Zhoushan and Shanghai clays have small ranges of C e values. Similar phenomena can be found for C si and C k . The increase in k 0 with the sampling angle is smaller for high plasticity Wenzhou clay than for low plasticity Zhoushan and Shanghai clays.

2. The initial anisotropy of clays before yielding has a significant influence on the anisotropic behavior of clays. The maximum value of C e varies very much and the value of C e corresponding to a vertical stress of 1600 kPa varies slightly with the deposition plane angle for all three clays. The compression index C c and swelling index C sr were found to vary slightly with the sampling angle, but not C si . The initial permeability k 0 depends much on the sampling angle for the three clays which are relevant to the initial anisotropy of clays.

The results of this study clearly show the anisotropic oedometric consolidation behavior of natural clays. It is meaningful to determine the anisotropic parameters associated with oedometric consolidation and use them for engineering designs and projects on clays. The perspectives of this study are now to relate these behaviors at the specimen size to microstructural analyses in order to better understand the physical origins of the material anisotropy.

Fig. 2 Fig. 1

 21 Fig. 2 Classification of clays by liquid limit and plasticity index CL: low plastic inorganic clays, sandy, and silty clays; OL: low plastic inorganic or organic silty clays; CH: high plastic inorganic clays; OH: high plastic fine sandy and silty clays

Fig

  Fig. 4 Orientation of clay particles for natural depositions on vertical plane: (a) Zhoushan clay; (b) Shanghai clay; (c) Wenzhou clay

Fig. 5 Fig. 6

 56 Fig. 5 Sample block of soil and sampling angle relative to natural deposition plane

Fig. 7 Fig. 8

 78 Fig. 7 Influence of sampling angle on (a) compression index, (b) swelling index for initial compression and (c) swelling index for unloading and recompression

Fig. 9

 9 Fig. 9 Relationship between C e and vertical stress for three clays sampled at different angles relative to their deposition planes: (a) Zhoushan clay, (b) Shanghai clay and (c) Wenzhou clay

-

Fig. 10

 10 Fig. 10 Relationship between C e and e in the normally consolidated range for three clays with different sampling angles: (a) Zhoushan clay, (b) Shanghai clay, and (c) Wenzhou clay

Fig. 11

 11 Fig. 11 Relationship between C e and deposition plane angle for three clays: (a) maximum C e ; (b) C e corresponding to a vertical stress of 1600 kPa

  Fig. 13 Evolution of permeability with void ratio for clays from (a) Zhoushan, (b) Shanghai and (c) Wenzhou

  Fig. 14 Relationship between k 0 and sampling angle for three clays (a); Relationship between C k and sampling angle for three clays (b)

Table 1 Physical properties of tested clays

 1 Shanghai clay and Zhoushan clay have low plasticity, while Wenzhou clay is high plastic clay.

	Clay	Depth, h (m)	Liquid limit, w L	Plasticity index, I P	Initial void ratio, e 0	Percentage of clay (%)	Percentage of silt (%)	Percentage of sand (%)
	Shanghai	12	42.5	20	1.08	26	72.4	1.6
	Zhoushan	11	40.7	20	1.13	37	57.8	5.2
	Wenzhou	8	63	35	1.98	44	54.8	1.2

  clays functions relative to the sampling angles. It is clear that, for Shanghai clay, the value of p0 is the biggest for the 0° sample and gradually decreases as the sampling angle increases. The difference between the maximum and minimum values of p0 for Shanghai clay is about 17 kPa. However, for Zhou-

	shan and Wenzhou clays, p0 was not strongly in-
	fluenced by the sampling angle. Furthermore, com-
	paring the average measured yield stress with the in
	situ vertical effective stress indicates that Zhoushan
	clay is underconsolidated, Shanghai clay is slightly
	overconsolidated (OCR= p0 / v0 =1.06) and Wenzhou
	clay is normally consolidated.

* Project supported by the National Natural Science Foundation of China (Nos. 41240024 and 41372285), the Shanghai Pujiang Talent Plan (No. 11PJ1405700), the Research Fund for the Doctoral Program of Higher Education of China (No. 20110073120012), and the European Project CREEP (No. PIAPP-GA-2011-286397)


M.Sc Sébastien Roujol

PhD Mario Ries

PhD Bruno Quesson

PhD Chrit Moonen

PhD Baudouin Denis De Senneville
email: baudouin@imf.u-bordeaux2.fr

Real-time MR-Thermometry and Dosimetry for interventional guidance on abdominal organs

The use of Proton Resonance Frequency shift (PRF) based MR-Thermometry for interventional guidance on abdominal organs is hampered by the constant displacement of the target due to the respiratory cycle and the associated thermometry artifacts.

Ideally, a suitable MR-Thermometry method should for this role achieve a sub-second temporal resolution while maintaining a precision comparable to those achieved on static organs, while not introducing significant processing latencies.

Here, a computationally effective processing pipeline for 2D image registration coupled with a multi-baseline phase correction is proposed in conjunction with high frame-rate MRI as a possible solution.

The proposed MR-thermometry method was evaluated for five minutes at a frame-rate of 10 images/s in the liver and the kidney of 11 healthy volunteers and achieved a precision of less than 2 °C in 70 % of the pixels while delivering temperature and thermal dose maps on the fly. The ability to perform MR-Thermometry and Dosimetry in-vivo during a real intervention was demonstrated on a pig kidney during an HIFU-heating experiment.

INTRODUCTION

MR-thermometry relying on the water proton resonance frequency (PRF) is gaining importance for monitoring and guiding thermal therapies such as radio-frequency [START_REF] Laumonier | Real-time monitoring of hepatocellular carcinoma radiofrequency ablation by quantitative temperature MRI[END_REF], laser [START_REF] Puls | Laser ablation of liver metastases from colorectal cancer with MR thermometry: 5-year survival[END_REF] or focused ultrasound thermal ablation [START_REF] Fry | Precision high intensity focusing ultrasonic machines for surgery[END_REF][START_REF] Hindley | MRI guidance of focused ultrasound therapy of uterine fibroids: early results[END_REF][START_REF] Hynynen | Preclinical testing of a phased array ultrasound system for MRI-guided noninvasive surgery of the brain-A primate study[END_REF]. Typically, PRF-based MR-thermometry relies on the voxel-wise evaluation of phase differences between sequentially acquired gradient echo images.

However, for the use on abdominal organs, this renders the method very sensitive to motion artifacts and magnetic field changes. These motion artifacts can be coarsely classed into the two following types: Intra-scan motion artifacts and Inter-scan motion artifacts.

Intra-scan motion artifacts are caused by displacement during the MR-acquisition process and lead to image blurring and object ghosting. Commonly this type of artifact is addressed by the use of fast MR-acquisition schemes. Alternatively, intra-scan motion artifacts can also be avoided by using respiratory gated sequences which, in turn, reduce the temporal resolution to the respiratory frequency.

Inter-scan motion artifacts are due to organ displacement between the MR-acquisitions.

Although the image quality is not degraded, this leads to a mis-registration between subsequent phase images and thus to artifacts in the subtraction process. Furthermore, since any displacement or plastic deformation of the abdominal organs will in general also lead to a modified demagnetization field and R1.1 R1.17 thus to a change of the local magnetic field [START_REF] Poorter | The proton resonance frequency shift method compared with molecular diffusion for quantitative measurement of two dimensional time dependent temperature distribution in phantom[END_REF][START_REF] Young | An evaluation of the effects of susceptibility changes on the water chemical shift method of temperature measurement in human peripheral muscle[END_REF][START_REF] Peters | Proton-resonance frequency shift MR thermometry is affected by changes in the electrical conductivity of tissue[END_REF], additional phase artifacts are introduced.

To overcome these problems, several correction strategies have been proposed, such as respiratory gating [START_REF] Weidensteiner | Stability of real-time MR temperature mapping in healthy and diseased human liver[END_REF], cardiac triggering [START_REF] Quesson | Dynamic cardiac MR temperature imaging with the Proton Resonance Frequency technique: a study on the feasibility and accuracy[END_REF], navigator echoes [START_REF] De Zwart | Online correction and visualization of motion during MRI-controlled hyperthermia[END_REF], multi-baseline acquisition to sample periodic changes [START_REF] Vigen | Triggered, navigated, multi-baseline method for proton resonance frequency temperature mapping with respiratory motion[END_REF][START_REF] Denis De Senneville | Real time adaptive methods for treatment of mobile organs by MRI controlled High Intensity Focused Ultrasound[END_REF] and reference-less phase corrections [START_REF] Rieke | Referenceless PRF shift thermometry[END_REF]. Furthermore, although the local temperature is a precise indicator for the energy deposition, it does not directly allow to estimate tissue damage and thus to determine the therapy endpoint. For this purpose, the concept of the equivalent thermal dose was introduced to reflect the biological effects of elevated temperatures on tissue by Sapareto et al [START_REF] Sapareto | Thermal dose determination in cancer therapy[END_REF]. The tissue destruction is achieved when the equivalent thermal dose exceeds the lethal dose (which is taken as 43 °C during 240 mn). Due to its exponential dependence on the temperature, artifacts in the temperature measurements have a large impact on the precision of the thermal dose estimate and will accumulate during the intervention due to the integration. Furthermore, since the temporal integration has to be performed on a voxel-by-voxel basis, a precise image registration is paramount for a precise calculation of the thermal dose.

Although the required correction schemes are computationally intensive, MR-image guidance of semi-or fully-automated interventions, such as high intensity focused ultrasound (HIFU) ablations, require the temperature and thermal dose measurements to be available in real-time with preferably short calculation latencies.

Here, a computationally efficient pipeline for 2D motion compensated PRF-thermometry and thermal dose measurements on moving abdominal organs is proposed. The method is designed to address both, inter-scan and intra-scan artifacts by applying high frame-rate MRI coupled with a realtime image registration of all incoming MR-images. Phase artifacts due to the periodical motion of the respiration cycle are addressed either by a multi-baseline strategy [START_REF] Vigen | Triggered, navigated, multi-baseline method for proton resonance frequency temperature mapping with respiratory motion[END_REF][START_REF] Denis De Senneville | Real time adaptive methods for treatment of mobile organs by MRI controlled High Intensity Focused Ultrasound[END_REF], or alternatively by applying a phase correction based on a model of the phase variation in dependence of the principal displacement component which is adjusted from training-data on a voxel-by-voxel basis [START_REF] Hey | Online correction of respiratory-induced field disturbances for continuous MRthermometry in the breast[END_REF].

Furthermore, temporal filtering of the temperature maps allows to choose a balance between temporal resolution and the precision of the MR-thermometry which can be freely adjusted according to the employed interventional modality and the available signal-to-noise ratio of the target area.

All computational intensive calculations are off-loaded to a dedicated graphics processing unit (GPU) to allow sustained MR-thermometry and dosimetry under real-time conditions over durations of several minutes with short image latencies.

The potential of the real-time reconstruction pipeline to remove MR-thermometry artifacts dynamically is evaluated in-vivo on the abdomen of 11 healthy volunteers under free-breathing conditions, and on a pig kidney during an HIFU-heating experiment.

MATERIAL AND METHODS

Figure 1 gives an overview of the employed data processing sequences for the multi-baseline and the linear phase modeling approaches.

Common to both approaches is the separation of the intervention into a pre-treatment step and a hyperthermia step. The pre-treatment step is performed prior to the intervention and records a reference dataset of magnitude and phase images during the motion cycle of the organ. In the multi-baseline approach this data is used to establish a lookup-table which serves as a reference phase during the hyperthermia step. For that purpose, the motion cycle must be sampled with a sufficient density in order to avoid discretization errors. With typical imaging frame-rates of 5-10 images per second and a respiration frequency of 3-5 s this pre-treatment step can be completed in a relatively short duration of 10-20 s. For the linear phase-modeling approach the data collection of the pre-treatment step is subsequently used to fit a model function on a voxel-by-voxel basis in order to relate the recorded phase variations to their principal displacement component (PDC).

During the subsequent hyperthermia step, all incoming images are first co-registered to the common position of the reference image [START_REF] Maintz | A survey of medical image registration[END_REF]. Subsequently, the phase image is corrected for local magnetic field changes based on the data collected in the pre-treatment step and the temperature is calculated. Then, temporal filtering based on an infinite impulse-response (IIR) filter ( 18) is applied to the temperature maps to increase the Signal-To-Noise ratio (SNR) and finally the thermal dose is calculated.

The proposed approach takes benefit of a combined CPU/GPU architecture by offloading computational intensive calculations to the GPU and thus freeing the CPU for pipeline management and data preparation.

Organ displacement estimation

In the presented work a region based PDC estimation is initially performed on the T 2 * weighted images. This method serves also for preconditioning the more complex optical flow based image registration. Note that all corrections are estimated on the magnitude images but applied to the complex MR-images in order to avoid interpolation problems with spatial phase wraps.

PDC estimation: An affine transformation of six free parameters (two translations, two scales, a rotation and a shear) is estimated using a differential Gauss-Newton approach (19). Since the complexity of the estimated deformation is limited, the algorithm is unsuitable to cope with fold-over MR-artifacts, mixtures of static/dynamic parts of the entire field-of-view and complex motion patterns.

Therefore it is restricted to a region of interest (ROI), which is manually set at the beginning of the intervention containing the full path of the targeted organ.

Optical flow based image registration:

The implemented approach was originally proposed by Horn and Schunck [START_REF] Schunck | Determining optical flow[END_REF] for optical image stabilization and computes a displacement field on a pixel-by-pixel basis. It includes a modification proposed by Cornelius and Kanade [START_REF] Cornelius | Adapting optical flow to measure object motion in reflectance and X-Ray image sequences[END_REF] which relaxes the intensity conservation condition from a global to a regional conservation. Finally, in order to optimize the computation time and to stabilize the convergence of the algorithm, a multi-resolution scheme was used [START_REF] Pratikakis | Robust multiscale deformable registration of 3D ultrasound images[END_REF] which iterates the registration algorithm from a four-fold down-sampled image step-by-step to the full image resolution.

Susceptibility and displacement related phase changes with motion

Although spatial displacements in phase images can be corrected by a suitable image registration, unwanted phase shifts generated by modified local magnetic susceptibility perturbations generally remain. This effect induces an apparent temperature modification which can bias or even completely mask the true temperature evolution induced by the energy deposition. Precise modeling of the inhomogeneous magnetic field in-vivo and under real-time conditions is difficult to achieve and thus several alternative simplified strategies have been proposed to allow to correct motion related errors in PRF-based MR thermometry [START_REF] De Zwart | Online correction and visualization of motion during MRI-controlled hyperthermia[END_REF][START_REF] Vigen | Triggered, navigated, multi-baseline method for proton resonance frequency temperature mapping with respiratory motion[END_REF][START_REF] Denis De Senneville | Real time adaptive methods for treatment of mobile organs by MRI controlled High Intensity Focused Ultrasound[END_REF][START_REF] Rieke | Referenceless PRF shift thermometry[END_REF].

In the presented work, two correction strategies have been implemented, a multi-baseline approach and a linear phase model approach, which are both based on a training set composed of the first N dynamic acquisitions in the time series (typically N=100), acquired prior to MR-thermometry.

Multi-baseline approach:

The first N dynamic acquisitions in the time series are stored in a look-up table [START_REF] Vigen | Triggered, navigated, multi-baseline method for proton resonance frequency temperature mapping with respiratory motion[END_REF][START_REF] Denis De Senneville | Real time adaptive methods for treatment of mobile organs by MRI controlled High Intensity Focused Ultrasound[END_REF], which holds magnitude images and the corresponding registered phase images. During the intervention, arriving magnitude images are matched to the image with the closest correspondence in the look-up table by computing and comparing the inter-correlation coefficients. This task is performed on the CPU parallel to the image registration process performed on the GPU. Subsequently, the matching phase image of the GPU collection is subtracted to obtain a corrected temperature map.

Linear phase model approach:

A second correction algorithm was implemented, which extends the multi-baseline approach by applying a model of local magnetic field changes to the phase data stored in the look-up table. The objective is to interpolate intermediate positions in order to reduce discretization errors and the required number of images in the correction data set, and to improve temperature accuracy by using a synthesized reference phase image. For this purpose, the first N acquisitions are used to evaluate a linear relation between motion and phase variations. To prevent errors arising from the 2π-periodicity of the phase, all phase values φ t (x,y) are first reordered over the time function of their corresponding displacement vector amplitude and subsequently phase unwrapped. Then, for each individual voxel, a system of N equations expressing the unwrapped registered phase (φ t ) as a linear combination of the six free transformation parameters (noted T i t , 0 ≤ i < 6 is the index of the estimated transformation parameters, and 0 ≤ t ≤ N-1 is the index of each image of the pre-treatment step) is given:

ϕ t  x,y = ∑ i=0 5 T i t ⋅P i  x,y  +P 6  x,y  , ∀ t, 0≤t≤N-1 [1]
where (x,y) denotes voxel coordinates, P i , (0 ≤ i < 6) are the unknowns and can be intuitively seen as a parameterized magnetic field model. This results in an overdetermined system resolved using a Singular Value Decomposition (SVD) .

For each image acquired after the learning phase, the transformation coefficients T i t (calculated from the magnitude image registration) are used with the P i (calculated from [START_REF] Laumonier | Real-time monitoring of hepatocellular carcinoma radiofrequency ablation by quantitative temperature MRI[END_REF]) to compute a synthetic reference phase map (φ reco ) as follows:

ϕ reco  x,y = ∑ i=0 5 T i t ⋅P i  x,y +P 6  x,y  [2]
This reference phase map is then subtracted from the acquired motion corrected phase image to suppress the background phase prior to the temperature calculation.

Correcting for local temperature aliasing

Since the 2π periodicity of the image phase can lead to aliasing artifacts in the temperature maps, a temporal phase unwrapping on a pixel-by-pixel basis is applied.

This process is valid under the condition that the temperature variation between two successive acquisitions does not create a phase variation greater than π. Note that on mobile targets, each phase image is registered to a common reference position and subsequently phase corrected, before temporal unwrapping is performed. Therefore, background-phase changes do not contribute to this limitation.

Correcting for magnetic field drift

Although recent MRI designs have excellent spatial magnetic field homogeneity, their field is not entirely stable over time when sustained high frame-rate imaging is applied [START_REF] El-Sharkawy | Monitoring and correcting spatiotemporal variations of the MR scanner's static magnetic field[END_REF]. Therefore, the subtraction of phase images acquired at different times can be affected with a different bias for each pixel [START_REF] Peters | Proton-resonance frequency shift MR thermometry is affected by changes in the electrical conductivity of tissue[END_REF][START_REF] Riederer | Recent technical advances in MR imaging of the abdomen[END_REF]. In this study, this perturbation is corrected by subtracting a global temperature offset (a magnitude-weighted average is used to give less importance in poor signal area) obtained from a region of interest which is chosen in the moving organ, adjacent to the ablation area.

Temporal temperature filtering

An infinite impulse-response (IIR) low-pass filter of 5 th order ( 25 On-line in-vivo heating study on a pig kidney: MRI guided HIFU was performed in vivo in the kidney of a pig under general anesthesia. No focal point position adjustment with respect to target displacements was performed. The MR sequence employed the following parameters: 1500 dynamic sagittal images, one slice, TR=100 ms, TE=41 ms, bandwidth in readout direction=2085 Hz, flip angle=35°, FOV=320×140 mm 2 , slice thickness=6 mm, matrix=128×56, using the integrated phased array coil of the HIFU system. The data processing sequence designed for the linear phase modelling approach was used to perform the MR-thermometry and dosimetry in real-time. and liver. Over the 11 human volunteers, the SNR was evaluated to 9.64 ± 2.4 (min=7.2, max=14.3) in the kidney and 7.5 ± 3.1 (min=4.5, max=14.3) in the liver. In average over all volunteers, the temperature stability is improved from an initial value of over 8 °C to 2.12 °C (kidney) and 2.66 °C (liver) using the multi-baseline approach, and to 1.5 °C (kidney) and 2.16 °C (liver), when the linear phase model is used. This precision can be furthermore improved by over 20 % if a drift correction is applied: While a temperature standard deviation of 1.51 °C (kidney) and 2.07 °C (liver) were obtained with the multi-baseline correction, the correction based on a linear phase model achieves a further reduction to 1.26 °C (kidney) and 1.77 °C (liver). Additional temporal filtering results in a final precision of 0.86 °C (multi-baseline), 0.79 °C (linear model) in the kidney and 1.05 °C, 0.98 °C in the liver. None of the data were rejected due to miss-registration or faulty phase correction due to the absence of spontaneous motion. Figure 3 shows as an example the least precise result of the examined volunteer group. Even in this case, both correction strategies ensured 2 °C of temperature stability in 70 % of all pixels of both the kidney and the liver.

RESULTS

Performance comparison between the CPU only and the CPU/GPU implementation

On-line in-vivo heating study on a pig kidney

The overall temperature accuracy over the whole kidney was found to be 0.65 °C ± 0.11 (min=0.4, max=0.99) with a measured SNR of 6 (although a SNR of 6 provides a theoritical temperature standard deviation of 1.47 °C, the temporal filtering reduces the initial temporal resolution of 10 Hz down to 2 Hz improving the temperature precision by a factor of √5). Figure 4.a shows the temperature image after 39 seconds of sonication corrected with the linear phase model approach. The evolution of the temperature at the focal point position is reported on figure 4.b. An hyperthermia of 12 °C was reached, which leads to a final thermal dose of 10 % of the lethal dose as shown in figure 4.c [START_REF] Seror | Quantitative magnetic resonance temperature mapping for real-time monitoring of radiofrequency ablation of the liver: an ex vivo study[END_REF].

DISCUSSION

Organ displacement estimation:

The main limitation of the presented approach is that through-plane motion can not be compensated, which places severe constrains on the type of motion that can be adressed and the possible slice placement. For the case of respiratory induced motion, which predominantly causes a linear displacement of the abdominal organs in head-foot direction and to a lesser extend in anteriorposterior direction, this limits the choice of the slice placement to the coronal or sagital orientation.

R1.11

R1.12 R2.4

Although a full 3D motion compensation would be preferable, this is in practise hard to achieve with image based compensation techniques due to the inherent compromise between acquisition time and volume coverage.

PDC estimation: This approach was found to be well suited for targets which show only little plastic deformation during the respiratory cycle, such as the kidney. However, since such targets are in general in the vicinity of either elastic soft-tissue, such as digestive tubes or liver, or static tissue, such as the quadratus lumborum muscle or the vertebral column, the algorithm needs to be confined to a ROI encompassing only the full path of the selected target and excluding such problematic areas. For targets which display a complex spatially variant deformation such as the liver, this approach was found to be unsuitable since a global adaption of only six free transformation parameters lead to a poor representation of the complex deformation. The main advantages of this approach are the short computation times and its robustness against noise and local intensity variations.

Optical flow based image registration: This technique allows to estimate complex organ deformations

with sub-pixel precision on a voxel-by-voxel basis. Furthermore, the underlaying regularity constraint [START_REF] Schunck | Determining optical flow[END_REF] which assumes a continuous differentiable motion field was found well suited for elastic targets which display local deformations during the respiratory cycle such as the liver or the spleen. Compared to Ross et al. [START_REF] Ross | Real-time liver motion compensation for MrgFUS[END_REF], the implemented approach does not require any landmark identification and thus the estimated motion field is independent of preceding user interactions or preparational segmentation steps.

Although this approach is computationally intensive and thus challenging to implement for realtime imaging, the proposed CPU/GPU implementation renders this method compatible with sustained high frame-rate real-time imaging.

The main disadvantage of the original approach, as suggested by Horn and Schunk (19), is that it relies on the preservation of the global signal amplitude within the image. This renders the method sensible to contrast variations and to the effect of residual through-plane motion since any temporal change in the global intensity of the image is attributed to motion. The stability of the algorithm can be greatly improved by relaxing the global regularity constraint as suggested by Cornelius and Kanade [START_REF] Cornelius | Adapting optical flow to measure object motion in reflectance and X-Ray image sequences[END_REF]. Although this modification limits mis-registration of the area of interest caused by contrast changes or tissue disappearance outside the target region, this condition can still be violated. During hyperthermia several MR relevant tissue properties such as T 1 and T 2 relaxation times can change during heating process and thus lead to local signal intensity variations in the heated region. The impact of this miss-registration effect on large ablation areas, which may be used for advanced ablation R2.5

strategies such as spatial volumetric temperature feedback [START_REF] Mougenot | Automatic spatial and temporal temperature control for MR-guided focused ultrasound using fast 3D MR thermometry and multispiral trajectory of the focal point[END_REF], needs to be carefully investigated in future studies.

Susceptibility related phase changes with motion:

The proposed multi-baseline correction for susceptibility related phase changes is designed to correct MR-thermometry artifacts on abdominal organs which are subject to a periodic displacement due to the respiratory cycle. Since phase variations are "learned" in a preceding learning step and subsequently applied to correct the MR-thermometry during the intervention, the method can intrinsically not correct for MR-thermometry artifacts associated with spontaneous motion. If during hyperthermia new positions are observed, a recalibration of the phase correction data is required, which is a disadvantage compared to other correction approaches such as reference-less correction approaches [START_REF] Rieke | Referenceless PRF shift thermometry[END_REF].

However, in practise this limitation is less severe, since the proposed correction approach is paired with high frame-rate imaging and thus a re-calibration can be completed in the relatively short time of two to five respiratory cycles.

Furthermore, compared to reference-less correction methods ( 14), the proposed multi-baseline correction allows an accurate correction of susceptibility related phase changes even in regions with complex susceptibility distributions or signal discontinuities such as organ boundaries, since it requires neither a uniform susceptibility distribution in the target nor time-consuming 2D spatial phase-unwrapping steps [START_REF] Chavez | Understanding phase maps in MRI: a new cutline phase unwrapping method[END_REF].

In addition, since all images are co-registered to a common reference position before the phase corrections are applied, temporal aliasing of the temperature arising from the 2π periodicity of the phase can simply be removed by applying a temporal unwrapping step on a voxel-by-voxel basis.

The proposed multi-baseline correction can be improved further by applying a linear phase model on a voxel-by-voxel basis with respect to the following three aspects: First, the correction is not constrained to positions present in the collection, but can also be interpolated to intermediate positions, which significantly reduces discretization errors and the required number of images in the correction data set [START_REF] Hey | Online correction of respiratory-induced field disturbances for continuous MRthermometry in the breast[END_REF]. Second, while the multi-baseline based correction can intrinsically not correct for motion amplitudes higher than the ones observed in the learning phase, the phase model allows extrapolation and can still provide an estimate of the reference phase. Finally, since the reduced set of six P i coefficients is derived from the resolution of an overdetermined system of 50-100 reference images, noise will be reduced on the resynthesized phase images. Ideally, the noise contribution on temperature uncertainty can thus be reduced by a factor √2 (as optimally σ(φ reco )=0). Since the lower precision boundary of PRF-based MR-thermometry on static targets is directly proportional to the inverse of the SNR R2.5

of the employed sequence [START_REF] Conturo The | Signal to Noise in Phase Angle Reconstruction : Dynamic Range Extension Using Phase Reference Offsets[END_REF], it is possible to evaluate the effectiveness of the correction strategies for dynamic phase artifacts by a direct comparison (assuming a noise free reference phase): On the 11 human volunteers the SNR study leads to a minimal achievable temperature standard deviation of 1.12 °C ± 0.26 (kidney, min=0.71, max=1.48) and 1.52 °C ± 0.51 (liver, min=0.71, max=2.22). The experimental study approaches those theoretical values when the linear phase model is used: After drift correction, the temperature precision in the kidney and the liver averaged over all volunteers were found to be 12 % and 16 % larger than the boundary imposed by the SNR, respectively (see figure 2).

On the other hand, the correction using the linear phase model introduces several new problems compared to the multi-baseline approach: The magnetic field perturbation is estimated assuming a simple linear magnetic field variation with organ displacement. Although this assumption holds in general for small displacements, the precision of this simple model must be carefully evaluated in regions displaying large motion amplitudes or large susceptibility variations, such as in the vicinity of the digestive tubes (see arrow (1) in figure 3.c), of the quadratus lumborum muscle (see arrow (2) in figure 3.c), of the vertebral column (see arrow (3) in figure 3.c), and in the upper part of the liver (see arrow (4) in figure 3.c). This can in practice be achieved by simply mapping the fitting error ε(x,y) from equation [START_REF] Laumonier | Real-time monitoring of hepatocellular carcinoma radiofrequency ablation by quantitative temperature MRI[END_REF]. This also allows to discard regions where low signal levels in conjunction with large susceptibility variations hampered a successful temporal phase unwrapping which is required to adjust the linear model to the phase data.

Rapid acquisition with subsequent temporal filtering:

The presented MR-Thermometry approach addresses intra-scan motion artifacts by applying high frame-rate MRI and the remaining inter-scan artifacts using image registration and a multibaseline phase correction. For the case of respiratory induced motion of abdominal organs, MR-acquisition times of typically 100 ms per slice or less are required to avoid intra-scan artifacts. This places severe restriction on the available processing time for image reconstruction, image registration and all MR-thermometry related calculations in order to maintain sustained real-time MR-thermometry.

Figure 1 shows that the motion estimation step remains the most time consuming step. For motion estimation based on/combined with optical flow, computation times can easily exceed typical MRacquisition times when the single CPU based implementation is used. The GPU based implementation allows to reduce the bulk computation time for all image processing steps below the used MR-acquisition time. This means that the achievable MR-thermometry frame-rate is below 30 images/s only limited by the constraints of the MR-sequence in order to satisfy the real-time condition. For interventional guidance such as energy deposition control or retroactive adjustment of the beam position, short laten-cies have a great importance. Denis de Senneville [START_REF] Denis De Senneville | Real time adaptive methods for treatment of mobile organs by MRI controlled High Intensity Focused Ultrasound[END_REF] showed that large image latencies have to be compensated with a temporal prediction and that the robustness of such algorithms improves with short latency. The two proposed corrections ensure a latency of 27.3 ms (image matrix of 128×128), which is suitable for interventional MR-guidance and is realizable on affordable commodity hardware.

Commonly, high frame-rate imaging for MR-thermometry is considered to be limited by the precision penalty imposed by the associated low SNRs on the temperature measurements. Since all MR temperature images in the time series are registered to the same reference position, temporal filtering can be applied to improve temperature accuracy prior to thermal dose calculation. This, in turn, allows to choose a balance between temporal resolution and the precision of the MR-thermometry which can be freely adjusted according to the employed interventional modality and the available SNR of the target area. For a direct tissue ablation using high temperature elevations of up to 30 °C, the temporal resolution will be privileged. On the other hand, recent experiments employing hyperthermia to control local drug delivery or gene expression [START_REF] Moonen | Spatio-temporal control of gene expression and cancer treatment using magnetic resonance imaging-guided focused ultrasound[END_REF], which represent a great potential for cancer therapy in kidney and liver, have much higher demands on the precision of the thermometry since it employs low temperature elevations of 4-5 °C. Theoretically, a perfect co-registration allows to double the SNR, and thus the precision of the temperature measurements, by reducing the bandwidth by a factor of 4. In the volunteer study, a filter with a transition band between 0.15 and 0.45 resulted in an improvement of the precision of MR-thermometry by a factor 1.85 which corresponds well with the expected theoretical result.

CONCLUSION

The proposed approach for 2D motion compensated MR-thermometry and dosimetry addresses both, inter-scan and intra-scan motion artifacts, by applying high-framerate MRI coupled with real-time image registration and multi-baseline phase correction of all incoming MR-images. This, in conjunction with the use of parallel processing on affordable commodity graphics hardware, allows to achieve a sub-second temporal resolution with very short image latencies over sustained imaging periods of several minutes which is purely limited by the constraints of the MR-acquisition.

The temperature artifacts related to the periodic respiratory motion of the abdominal organs were found to be reduced to the boundary imposed by the SNR of the employed sequence in both, kidney and liver. Additional temporal filtering of the temperature maps allows to freely readjust the balance between temporal resolution and additional precision of MR-Thermometry.

The method was found robust and artifact free in all examined cases and well able to follow the temperature evolution of an in-vivo HIFU ablation. This renders the method well suitable for the MR-

  ), using an elliptic approximation in the denominator of the transfer function (also referred to as a Cauer filter), was used for temporal low-pass filtering of the calculated temperature values. The filter coefficients were designed using the signal processing toolbox of MATLAB with the following characteristics: f pass- band =1.5 Hz to 8 Hz, f stop-band =4.5 Hz to 8.5 Hz, ripple pass-band: 3 dB, ripple stop-band: 50 dB. Passband and stop-band were adjusted depending on the dynamic scan-time of the employed sequence to result in an overall temporal resolution of 2 Hz.Experimental setupMRI imaging: Dynamic MR temperature imaging was performed on a Philips Achieva 1.5 T (Philips Healthcare, Best, The Netherlands) with a single-shot gradient recalled echo-planar (EPI) sequence[START_REF] Bruder | Image reconstruction for EPI with nonequidistant k-space sampling[END_REF]. MR images were obtained using an in house developed reconstructor.HIFU System: A MR compatible Philips HIFU prototype (Philips Healthcare, Helskinki, Finland) composed of a table top containing a 256 elements HIFU transducer (Imasonic, Besançon), integrated in the 1.5 T Achieva-Intera MRI was used to perform a temperature elevation. The transducer radius and aperture were 120 mm and 126 mm, respectively, creating an ellipsoid focal point (1×1×7 mm 3 ).Reconstructor description:The reconstructor was a dual processor (INTEL 3.1 GHz Penryn, four with 8 GB of RAM and dual 1 GB/s network interface cards. The GPU was a NVIDIA GTX280 card with 1 GB of DRAM connected over a PCIe x16 link. The GPU implementation was realized using CUDA[START_REF]NVIDIA CUDA: Compute Unified Device Architecture[END_REF].On-line in-vivo study: Dynamic MRI was performed under free-breathing conditions on the abdomen of 11 healthy volunteers under real-time conditions. An imaging frame-rate of 10 images/s was maintained for 300 seconds of MR-imaging while MR-Thermometry was performed in real-time. The MR sequence employed the following parameters: 3000 dynamic sagittal images, one slice, TR=100 ms, TE=26 ms, bandwidth in readout direction=2085 Hz, flip angle=35°, FOV=256×168 mm 2 , slice thick-ness=6 mm, matrix=128×84, using a four element phased array body coil. Statistical evaluation of the temperature stability was performed on the kidney and the liver of each volunteer individually by averaging the temporal temperature standard deviation over a ROI which was manually set in an area with maximal SNR and avoiding areas showing a complex susceptibility distribution (such as organ boundaries or major vessels).

Figure 1

 1 Figure 1 details the computation time of each processing step of the two proposed pipelines for an image resolution of 128×128 pixels between the CPU only and the CPU/GPU implementation. The most time consuming task, the image registration, was accelerated by a factor of 3 for the PDC estimation and a factor of 10 for the optical flow computation. The resulting overall latency for the entire pipeline (including 13 ms for data transport and 1.2 ms for image reconstruction) was reduced from 95 ms (CPU only) to 27.3 ms (CPU/GPU).

Figure 2

 2 Figure 2 details the precision improvement of each separate phase correction step for kidney

EC-FP6-project LSHB-CT-2005-512146 "Agence National de Recherche (project MRgHIFU-ALKT)", Fondation InNaBioSanté (project ULTRAFIT)", Philips Medical System.

guidance of hyperthermia ablations in abdominal organs under free-breathing conditions and as the basis for more advanced automatic spatial and temporal temperature control algorithms used in conjunction with dynamic ultrasound beam-steering.

LIST OF FIGURES